Explicit formula for square root of a $3\times3$ positive definite matrix I have an algorithm for calibrating a vector magnetometer. The input is $N$ readings of the $x$, $y$, $z$ axes: $(x_1, x_2, \dots, x_N)$, $(y_1, y_2, \dots, y_n)$, and $(z_1, z_2, \dots, z_N)$. The algorithm fits an ellipsoid to the data by estimating a symmetric $3\times3$ matrix $A$. In order to calibrate the system, it needs to calculate $\sqrt{A}$. I am adapting the algorithm for a microcontroller with very little memory, so cannot load standard matrix manipulation libraries. Is there an explicit formula for calculating the square root of $3\times3$ positive definite matrix?
 A: There may not be an explicit formula, but there may be some methods that give a fair approximation to the square root.
For instance, if you know that the eigenvalues of $A$ are in $[0,1]$, you can take a polynomial approximation $P(x)$ to $\sqrt{x}$ on $[0,1]$, and then consider $P(A)$ that approximates $\sqrt{A}$.
We can always reduce to "eigenvalues of $A$ in $[0,1]$" by multiplying by some positive constant.
Another way is to use the Newton approximation to the solution of the equation $x^2 = A$. Start with an initial term say $x_0 = I$, $x_{n+1} = \frac{1}{2}( x_n+ A \cdot x_n^{-1})$, it should converge quickly to $\sqrt{A}$.
A: I would rather use an iterative method (e.g. QR iteration or Jacobi's algorithm) to diagonalise $A$ and find its square root, but there is a semi-closed-form formula for $B=\sqrt{A}$. By Cayley-Hamilton theorem,
$$
B^3-\operatorname{tr}(B)B^2+\operatorname{tr}(B^{-1})\det(B)B-\det(B)I=0.\tag{1}
$$
Multiply both sides of $(1)$ by $B$, we obtain
$$
B^4-\operatorname{tr}(B)B^3+\operatorname{tr}(B^{-1})\det(B)B^2-\det(B)B=0.\tag{2}
$$
Substitute $(1)$ into $(2)$ and using the fact that $B^2=A$, we obtain
$$
A^2+\operatorname{tr}(B)\left[-\operatorname{tr}(B)A+\operatorname{tr}(B^{-1})\det(B)B-\det(B)I\right]+\operatorname{tr}(B^{-1})\det(B)A-\det(B)B=0.
$$
Therefore
\begin{align}
B&=\frac{-A^2+\left(\operatorname{tr}(B)^2-\operatorname{tr}(B^{-1})\det(B)\right)A+\operatorname{tr}(B)\det(B)I}
{\left(\operatorname{tr}(B)\operatorname{tr}(B^{-1})-1\right)\det(B)}\\
&=\frac{-A^2+(\alpha+\beta)A+\gamma\delta I}
{\alpha\gamma-\delta}\tag{3}\\
\end{align}
where
\begin{align}
\alpha&=\sqrt{\lambda_1\lambda_2}+\sqrt{\lambda_2\lambda_3}+\sqrt{\lambda_3\lambda_1},\\
\beta&=\operatorname{tr}(A),\\
\gamma&=\sqrt{\lambda_1}+\sqrt{\lambda_2}+\sqrt{\lambda_3},\\
\delta&=\sqrt{\det(A)}
\end{align}
and $\lambda_1,\lambda_2,\lambda_3$ are the eigenvalues of $A$.
So, to find $B$ using $(3)$, we need the eigenvalues $A$. Let $s=\operatorname{tr}(A^2)$. The characteristic equation of $A$ is then
$$
x^3-\beta x^2+\frac12(\beta^2-s)x-\delta^2=0\tag{4}
$$
and $\lambda_1,\lambda_2,\lambda_3$ are its roots.
Now $(4)$ can be solved by radicals.
