# How to obtain the closed form expression of least-square sphere fitting?

I see here and there that there is a closed form expression of the least-square fit of a sphere of radius $$r$$ and center $$\mathbf{c}$$ to $$N$$ data points $$\{\mathbf{x_i}\}_{i\in(1,\cdots,N)}$$.

How is this expression obtained?

I considered this energy to be minimized: $$E(r,\mathbf{c}) = \sum_{i=1}^N (\|\mathbf{x}_i-\mathbf{c}\|^2 - r^2)^2,$$ then I tried to find radius $$r^*$$ and center $$c^*$$ that minimize $$E$$, so that it makes its gradient null. I obtained: $$\frac{\partial E}{\partial r}(r^*,\mathbf{c}^*) = 0 \ \Rightarrow \ -4\sum_{i=1}^N (\|\mathbf{x}_i - \mathbf{c}^*\|^2 - {r^*}^2) r^* = 0 \ \Rightarrow\ {r^*}^2 = \frac{1}{N}\sum_{i=1}^N \|\mathbf{x}_i - \mathbf{c}^*\|^2 ,$$ which match what I found on the cited references above.

But I couldn't find any correct way to obtain the system of equation $$A\mathbf{c} = \mathbf{b}$$ in order to find the center $$\mathbf{c}^* = (A^T A)^{-1}A^T\mathbf{b}$$.

Is it by expressing the partial deivative of $$E$$ with respect to $$\mathbf{c}$$? $$\frac{\partial E}{\partial \mathbf{c}} (r^*,\mathbf{c}^*) = 0 \ \Rightarrow\ -4\sum_{i=1}^N(\|\mathbf{x}_i-\mathbf{c}^*\|^2 - {r^*}^2)(\mathbf{x}_i - \mathbf{c}^*) = 0.$$ I tried to develop further this equation by replacing $$r^*$$ by its expression above but I couldn't put it in the form $$A\mathbf{c}^* = \mathbf{b}$$

$$\def\vec{\boldsymbol}\def\peq{\mathrel{\phantom{=}}{}}$$Because\begin{align*} \frac{\partial E}{\partial \vec{c}}(r, \vec{c}) &= -4\sum_{k = 1}^n (\|\vec{x}_k - \vec{c}\|^2 - r^2)(\vec{x}_k - \vec{c})\\ &= 4\left( \sum_{k = 1}^n (\|\vec{x}_k - \vec{c}\|^2 - r^2) \right) \vec{c} - 4\sum_{k = 1}^n (\|\vec{x}_k - \vec{c}\|^2 - r^2) \vec{x}_k\\ &= 4\left( \color{blue}{\sum_{k = 1}^n \|\vec{x}_k - \vec{c}\|^2 - nr^2} \right) \vec{c} + 4r^2 \sum_{k = 1}^n \vec{x}_k - 4\sum_{k = 1}^n \|\vec{x}_k - \vec{c}\|^2 \vec{x}_k, \end{align*} Plugging $${r^*}^2 = \dfrac{1}{n} \sum\limits_{k = 1}^n \|\vec{x}_k - \vec{c}^*\|^2$$ into $$\dfrac{\partial E}{\partial \vec{c}}(r^*, \vec{c}^*) = \vec{0}$$ yields$$\begin{gather*} \sum_{k = 1}^n \|\vec{x}_k - \vec{c}^*\|^2 \vec{x}_k - \frac{1}{n} \left( \sum_{k = 1}^n \|\vec{x}_k - \vec{c}^*\|^2 \right) \left( \sum_{k = 1}^n \vec{x}_k \right) = \vec{0}.\tag{1} \end{gather*}$$ Since $$\|\vec{x}_k - \vec{c}^*\|^2 = \|\vec{x}_k\|^2 + \|\vec{c}^*\|^2 - 2\langle\vec{x}_k, \vec{c}^*\rangle$$ for all $$1 \leqslant k \leqslant n$$, then$$\begin{gather*} \sum_{k = 1}^n \|\vec{x}_k - \vec{c}^*\|^2 \vec{x}_k = \sum_{k = 1}^n (\|\vec{x}_k\|^2 + \|\vec{c}^*\|^2 - 2\langle\vec{x}_k, \vec{c}^*\rangle) \vec{x}_k\\ = \sum_{k = 1}^n \|\vec{x}_k\|^2 \vec{x}_k + \|\vec{c}\|^2 \sum_{k = 1}^n \vec{x}_k - 2\sum_{k = 1}^n \langle\vec{x}_k, \vec{c}\rangle \vec{x}_k, \end{gather*}$$$$\begin{gather*} \left( \sum_{k = 1}^n \|\vec{x}_k - \vec{c}^*\|^2 \right) \left( \sum_{k = 1}^n \vec{x}_k \right) = \left( \sum_{k = 1}^n (\|\vec{x}_k\|^2 + \|\vec{c}^*\|^2 - 2\langle\vec{x}_k, \vec{c}^*\rangle) \right) \left( \sum_{k = 1}^n \vec{x}_k \right)\\ = \left( \sum_{k = 1}^n \|\vec{x}_k\|^2\right) \left( \sum_{k = 1}^n \vec{x}_k \right) + n\|\vec{c}\|^2 \sum_{k = 1}^n \vec{x}_k - 2\left( \sum_{k = 1}^n \langle\vec{x}_k, \vec{c}\rangle \right) \left( \sum_{k = 1}^n \vec{x}_k \right), \end{gather*}$$ and (1) is equivalent to$$\begin{gather*} 2\sum_{k = 1}^n \left( \langle\vec{x}_k, \vec{c}\rangle - \frac{1}{n} \sum_{j = 1}^n \langle\vec{x}_j, \vec{c}\rangle \right) \vec{x}_k = \sum_{k = 1}^n \|\vec{x}_k\|^2 \vec{x}_k - \frac{1}{n} \left( \sum_{k = 1}^n \|\vec{x}_k\|^2 \right) \left( \sum_{k = 1}^n \vec{x}_k \right).\tag{2} \end{gather*}$$ Now,\begin{align*} &\peq \sum_{k = 1}^n \left( \langle\vec{x}_k, \vec{c}\rangle - \frac{1}{n} \sum_{j = 1}^n \langle\vec{x}_j, \vec{c}\rangle \right) \vec{x}_k = \sum_{k = 1}^n \langle \vec{x}_k - \vec{\overline{x}}, \vec{c} \rangle \vec{x}_k\\ &= \sum_{k = 1}^n ((x_k - \overline{x}) x_c + (y_k - \overline{y}) y_c + (z_k - \overline{z}) z_c) \begin{bmatrix}x_k \\ y_k \\ z_k\end{bmatrix}\\ &= \begin{bmatrix} x_c \sum\limits_{k = 1}^n x_k (x_k - \overline{x}) + y_c \sum\limits_{k = 1}^n x_k (y_k - \overline{y}) + z_c \sum\limits_{k = 1}^n x_k (z_k - \overline{z})\\ x_c \sum\limits_{k = 1}^n y_k (x_k - \overline{x}) + y_c \sum\limits_{k = 1}^n y_k (y_k - \overline{y}) + z_c \sum\limits_{k = 1}^n y_k (z_k - \overline{z})\\ x_c \sum\limits_{k = 1}^n z_k (x_k - \overline{x}) + y_c \sum\limits_{k = 1}^n z_k (y_k - \overline{y}) + z_c \sum\limits_{k = 1}^n z_k (z_k - \overline{z}) \end{bmatrix}\\ &= \begin{bmatrix} \sum\limits_{k = 1}^n x_k (x_k - \overline{x}) & \sum\limits_{k = 1}^n x_k (y_k - \overline{y}) & \sum\limits_{k = 1}^n x_k (z_k - \overline{z})\\ \sum\limits_{k = 1}^n y_k (x_k - \overline{x}) & \sum\limits_{k = 1}^n y_k (y_k - \overline{y}) & \sum\limits_{k = 1}^n y_k (z_k - \overline{z})\\ \sum\limits_{k = 1}^n z_k (x_k - \overline{x}) & \sum\limits_{k = 1}^n z_k (y_k - \overline{y}) & \sum\limits_{k = 1}^n z_k (z_k - \overline{z}) \end{bmatrix} \begin{bmatrix}x_c \\ y_c \\ z_c\end{bmatrix}, \end{align*}$$\begin{gather*} \sum_{k = 1}^n \|\vec{x}_k\|^2 \vec{x}_k = \sum_{k = 1}^n (x_k^2 + y_k^2 + z_k^2) \begin{bmatrix}x_k \\ y_k \\ z_k\end{bmatrix} = \begin{bmatrix} \sum\limits_{k = 1}^n (x_k^2 + y_k^2 + z_k^2) x_k\\ \sum\limits_{k = 1}^n (x_k^2 + y_k^2 + z_k^2) y_k\\ \sum\limits_{k = 1}^n (x_k^2 + y_k^2 + z_k^2) z_k \end{bmatrix},\\ \frac{1}{n} \left( \sum_{k = 1}^n \|\vec{x}_k\|^2 \right) \left( \sum_{k = 1}^n \vec{x}_k \right) = \sum_{k = 1}^n \|\vec{x}_k\|^2 \vec{\overline{x}}\\ = \sum_{k = 1}^n (x_k^2 + y_k^2 + z_k^2) \begin{bmatrix}\overline{x} \\ \overline{y} \\ \overline{z}\end{bmatrix} = \begin{bmatrix} \sum\limits_{k = 1}^n (x_k^2 + y_k^2 + z_k^2) \overline{x}\\ \sum\limits_{k = 1}^n (x_k^2 + y_k^2 + z_k^2) \overline{y}\\ \sum\limits_{k = 1}^n (x_k^2 + y_k^2 + z_k^2) \overline{z} \end{bmatrix}, \end{gather*}$$ then (2) is equivalent to$$2\begin{bmatrix} \sum\limits_{k = 1}^n x_k (x_k - \overline{x}) & \sum\limits_{k = 1}^n x_k (y_k - \overline{y}) & \sum\limits_{k = 1}^n x_k (z_k - \overline{z})\\ \sum\limits_{k = 1}^n y_k (x_k - \overline{x}) & \sum\limits_{k = 1}^n y_k (y_k - \overline{y}) & \sum\limits_{k = 1}^n y_k (z_k - \overline{z})\\ \sum\limits_{k = 1}^n z_k (x_k - \overline{x}) & \sum\limits_{k = 1}^n z_k (y_k - \overline{y}) & \sum\limits_{k = 1}^n z_k (z_k - \overline{z}) \end{bmatrix} \begin{bmatrix}x_c \\ y_c \\ z_c\end{bmatrix} = \begin{bmatrix} \sum\limits_{k = 1}^n (x_k^2 + y_k^2 + z_k^2) (x_k - \overline{x})\\ \sum\limits_{k = 1}^n (x_k^2 + y_k^2 + z_k^2) (y_k - \overline{y})\\ \sum\limits_{k = 1}^n (x_k^2 + y_k^2 + z_k^2) (z_k - \overline{z}) \end{bmatrix}.$$