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}$…
 A: $\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}.
$$
