Find the square root of a 3x3 matrix

Question

Let ${\bf A}$ be a symmetric matrix, and I would like to find $\sqrt{\bf A}$ of a $3\times 3$ matrix. I am interested in finding the $\sqrt{\rm A}$ which is equivalent to the "sqrtm" built-in function in python (which is based on the Schur method: LINK).

I found this formulation, $$ \sqrt{\bf A} = \frac{{\bf A} + \sqrt{\operatorname{det}{\bf A}}\ {\tt I}} {\sqrt{\,{\operatorname{tr}\left(\bf A\right) + 2 \sqrt{\,{\operatorname{det}\left(\bf A\right)}\,}}}} $$ in one of the previous posts LINK. This formulation gives the exact results as the sqrtm in python. However, it is valid for ONLY a $2\times2$ matrix.

My question: is there a similar formulation that can be valid for a $3\times3$ matrix?

Many thanks in advance!

Answer 1

Your link to the sqrtm method has a link to that function's source.

The algorithm does a Schuur decomposition of $A$, finding $Z$ and $T$ such that $A=ZTZ^*$ with $Z$ unitary and $T$ is upper triangular.

Then it tries to find $\sqrt{T}$ using _sqrtm_triu. This method is for finding the square root of an upper triangular. Once it find that (if it succeeds) it returns $\sqrt{A}$ as $Z\sqrt{T}Z^*$.

So now the question is about how _sqrtm_triu works. Consider that if you know $T$ is upper triangular, then (as long as you allow non-real entries) it is easy to find the diagonal entries of an upper triangular $\sqrt{T}$: they have to be square roots of the diagonal entries of $T$. After these are established, you can solve for the entries of the first off-diagonal. And repeat. I can't think of a nice way to write this in a succinct formula at present.

Answer 2

I'm not sure about a formulation for all $3\times 3$ matrices, but we can pretty easily find the square root of diagonalizable matrices. Suppose $$A = \begin{bmatrix} \\ v_1 & v_2 & v_3 \\ \\\end{bmatrix} \begin{bmatrix} \lambda_1 & 0 & 0\\ 0 & \lambda_2 & 0\\ 0 & 0 & \lambda_3 \end{bmatrix} \begin{bmatrix} \\ v_1 & v_2 & v_3 \\ \\\end{bmatrix}^{-1}$$ For eigenvalues $\lambda_i$ and eigenvectors $v_i$. Then $$\sqrt{A}= \begin{bmatrix} \\ v_1 & v_2 & v_3 \\ \\\end{bmatrix} \begin{bmatrix} \sqrt{\lambda_1} & 0 & 0\\ 0 & \sqrt{\lambda_2} & 0\\ 0 & 0 & \sqrt{\lambda_3} \end{bmatrix} \begin{bmatrix} \\ v_1 & v_2 & v_3 \\ \\\end{bmatrix}^{-1}$$ since $$\sqrt{A}^2 = \begin{bmatrix} \\ v_1 & v_2 & v_3 \\ \\\end{bmatrix} \begin{bmatrix} \sqrt{\lambda_1} & 0 & 0\\ 0 & \sqrt{\lambda_2} & 0\\ 0 & 0 & \sqrt{\lambda_3} \end{bmatrix} \begin{bmatrix} \sqrt{\lambda_1} & 0 & 0\\ 0 & \sqrt{\lambda_2} & 0\\ 0 & 0 & \sqrt{\lambda_3} \end{bmatrix} \begin{bmatrix} \\ v_1 & v_2 & v_3 \\ \\\end{bmatrix}^{-1} = A.$$

Answer 3

I found myself in the same situation, had to find a real symmetric square root $B$ of a real 3x3 symmetric positive definite matrix $A$, and that, on a small board with basically no memory, javascript, and very few libraries.

Note that if the matrix has negative eigenvalues, then there is at least no symmetric and real square root.

First, find the eigen values, say $\lambda_1$, $\lambda_2$, $\lambda_3$. Since the characteristic polynomial is of degree $3$, we can find the roots numerically, using either roots or trigonometric functions, see the wikipedia page for the cubic equation. So, set $\mu_i = \sqrt{\lambda_i}$, and define an interpolating polynomial $$ q(t) = \frac{ (t-\lambda_2)(t-\lambda_3)} { (\lambda_1-\lambda_2)(\lambda_1-\lambda_3) } \mu_1 + \frac{ (t-\lambda_3)(t-\lambda_1)} { (\lambda_1-\lambda_2)(\lambda_1-\lambda_3) } \mu_2 + \frac{ (t-\lambda_1)(t-\lambda_2)} { (\lambda_1-\lambda_2)(\lambda_1-\lambda_3) } \mu_3 $$ and set $$ B = q(A). $$ We now note that if $v$ is an eigenvector for $A$ with eigenvalue $\lambda_i$, then $$ Bv = q(A)v = q(\lambda_i)v = \mu_iv $$ and so $v$ is an eigenvector for $B$ with eigenvalue $\mu_i$. It follows from this, since $B$ is diaonalizable, and so $\mathbb{R}^3$ is generated by eigenvectors, that $B^2 = A$.

Note 1:

In this proof, you can really replace $\mu_i$ with any number, and get a new matrix using the polynomial $q$. In particular, you can set $\mu_i = \lambda^{a/b}$ to get a rational power of the matrix $A$.

Note 2:

In the formula for $q(t)$, we divide by differences of eigenvalues. As a result, the $\lambda_i$ must be distict. This method must be modified (simplified really) if there are only two or one eigenvalues. However, if you are implementing this using floating numbers, and you end up with eigenvalues that are very close, then you end up dividing and multiplying with very small numbers, which can be a problem with floating point arithmetic, I am not an expert on this, but be careful :)

Note 3:

This method works for matrices of any size, even nonsymmetric ones, as long as they are diagonalizable. The only problem is finding the roots of the characteristic polynomial, but there should exist simple numerical methods for this (again, not an expert on that topic).

Note 4:

This result can be generalized to basically any matrix, nonsymmetric, nondiagonalizable, etc. This is described in one of the answers here:

Find the square root of a matrix