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}. $$