# Algorithm for taking square root of a matrix

Suppose $$A = A^T$$ and suppose the entries of $$A$$ are in $$\mathbb{Z}^+$$. I want to find all matrices $$M$$ with natural entries so that: $$M^2 = A$$ How can one do this? I know techniques that will get a square-root of an arbitrary matrix, but I want the full set. I want to be able to do this efficiently for large matrices ~$$100 \times 100$$.

Of course, the set must be finite because we are working over positive integers and the matrix $$A$$ gives upper bounds.

• Let's just be clear here: if $A = I_n$, the $n \times n$ identity, then there are $2^n$ possible diagonal square roots, each $1$ in a bit-string of length $n$ corresponding to a $-1$ entries on the diagonal, and each $0$ entry corresponding to a $+1$ on the diagonal. So enumeration of "SOME of the roots" in the simplest possible case takes time $2^{100} \approx 10^{31}$. I think you're probably not going to get a quick answer. – John Hughes Sep 19 '20 at 21:59

## 2 Answers

Find an invertible matrix $$B$$ and diagonal matrix $$D$$ such that $$D=B^{-1}AB$$. Then take all square roots $$D_1,...,D_m$$ (all of them diagonal, there are $$m\le 2^n$$ of these where $$n$$ is the size of $$A$$ because every non-negative number has at most 2 square roots) of $$D$$ and all matrices $$BD_iB^{-1}$$. Look which ones of these matrices are over natural numbers.

• Oh, I just added that to the question. I'm looking for efficient algorithms for this. – mtheorylord Sep 19 '20 at 21:48
• I have changed the answer. Now you have fewer choices. Don't think it can be improved much. – Mark Sapir Sep 19 '20 at 21:59

In a way, if you can generate one, the rest follows, since if you have one solution any other matrix $$M'$$ will have to adhere to the following set of equations: $$(M-M')(M+M')=M^2-M'^2 =0$$ Now, by the definition of $$M,\;M'$$ they are simultaneously orthogonally diagonalizable, thus you can assume WLOG they are diagonal. even more, you know the eigen values of $$M,\;M'$$ sutisfy $$\mu_i^2=\mu'^2_i = \lambda_i\to \mu_i=\mu'_i$$ (These are the eigenvalues respectively of $$M,\;M'$$ and $$\lambda$$ is of $$A$$). Hence to get another matrix that answers your question is equivalent to the condition $$(\mu_i+\mu_i')\cdot(\mu_i-\mu_i')=0$$ For all of the eigenvalues. This gives a complete classification of solutions.

• Checking the integrality of the requested matrix can be done efficiently with tools of field theory. – Alon Yariv Sep 19 '20 at 22:01
• This answer can be generalized for not necessarily square matrices using SVD decomposition – Alon Yariv Sep 19 '20 at 22:04
• Well...it's a complete classification, but not a complete enumeration; you'd need to write down every orthogonal matrix by which you could conjugate them to get a complete answer. But you'd have to write down only those that conjugate them to integral matrices. – John Hughes Sep 19 '20 at 22:47
• @JohnHughes Well, I did what I normally do, reducing the question to an integral point on an algebraic variety. – Alon Yariv Sep 19 '20 at 23:08