What methods exist for finding the squareroot of a nondiagonalizabe positive complex matrix?

5$\begingroup$ If, by positive, you mean that $\langle Ax,x\rangle$ is nonnegative for all vectors $x$, then the matrix is diagonalizable. So your question might need rethinking $\endgroup$ – Yemon Choi Feb 4 '10 at 3:07

1$\begingroup$ For a Jordan cell $A$ with eigenvalue $t\neq 0$ write $A=tB$. $B$ has 1's on the main diagonal and $1/t$'s immediately above it. $N=BI$ is nilpotent, so a square root $C$ of $B=I+N$ can be found using the binomial formula (which gives a finite sum). Then $\sqrt{t}C$ will be a square root of $A$. If $A$ is a Jordan $n$ by $n$ cell with eigenvalue 0, then $A$ has no square roots for $n>1$ (for rank reasons). The case when $A$ is arbitrary follows from the above. $\endgroup$ – algori Feb 4 '10 at 3:54

1$\begingroup$ Wikipedia has a little section on non symmetric/hermitian positive matrices at en.wikipedia.org/wiki/… $\endgroup$ – Mariano SuárezÁlvarez Feb 4 '10 at 3:55

2$\begingroup$ It should be noted that, in general, algori's procedure gives many different square roots.. $\endgroup$ – Mariano SuárezÁlvarez Feb 4 '10 at 3:59

4$\begingroup$ Would someone please step up and give their comments as an answer? :) $\endgroup$ – Pete L. Clark Feb 4 '10 at 4:14
Following Pete's advice, here is my comment with some more details added:
For a Jordan block $A$ with eigenvalue $t\neq 0$ write $A=tB$. $B$ has 1's on the main diagonal and $1/t$'s immediately above it. $N=BI$ is nilpotent, so a square root $C$ of $B=I+N$ can be found using the binomial formula (which gives a finite sum). Then $\sqrt{t}C$ will be a square root of $A$.
If $A$ is arbitrary, then find the Jordan form $B$ of $A$ so that $A=C^{1}BA$. If there are no zero eigenvalues, then we can find a square root of each block and then conjugate back.
If there are Jordan blocks with eigenvalue 0, the problem gets a bit trickier. The square of a Jordan $m$ by $m$ block with zero eigenvalue is conjugate to the union of two $m/2$ by $m/2$ blocks if $m$ is even and to the union of an $(m1)/2$ by $(m1)/2$ block and an $(m+1)/2$ by $(m+1)/2$ block if $m$ is odd. This allows one to compute a square root of a union of two Jordan blocks of equal sizes or of a union of an $n$ by $n$ block and an $(n+1)$ by $(n+1)$ block.
Let $a_1\leq \ldots\leq a_k$ be the sizes of the zero eigenvalue Jordan blocks (including 1 by 1 ones, so $\sum a_i$ is the dimension of the generalized eigenspace with eigenvalue 0). $A$ has a square root, iff $a_1\ldots,a_k$ can be obtained from a sequence $b_1\leq \ldots\leq b_l$ of positive integers by replacing an even $m$ with $m/2,m/2$, an odd $m$ with $(m1)/2,(m+1)/2$ and leaving 1's untouched. (I know this looks messy but can't think of anything better.)
Of course, a square root of a matrix is not unique (if it exists).
Note that if all eigenvalues of $A$ are positive, then numerically it's probably easier to use the binomial formula straight away:
$$\sqrt{A}=\sqrt{t}(I+\frac{1}{2t}X\frac{1}{8t^2}X^2 +\cdots).$$
Here $X=AtI$ and the formula is valid for $t$ greater then the maximum eigenvalue of $A$.

1$\begingroup$ After some googling, I've found this paper: maths.manchester.ac.uk/~nareports/narep89.pdf. The paper should contain an efficient algorithm for real matrices, and references for algorithms for complex matrices. It seems to me that using a polynomial function for the square root might be inefficient (a lot of matrix multiplications). $\endgroup$ – user2734 Feb 4 '10 at 7:45

1$\begingroup$ unknown  here we have to find powers of the same matrix, not to multiply arbitrary matrices; computing powers is cheap e.g. one can find a Jordan form first (which the general method involves anyway). $\endgroup$ – algori Feb 4 '10 at 15:03
I'm not sure if you're looking for an efficient algorithm or just a method that works. In the latter case, we can ignore the positivity condition and just look for square roots of general complex matrices. algori essentially gave the answer in the comments: if you diagonalize the matrix, the Jordan blocks with nonzero eigenvalues automatically have (nonunique) square roots.
The nilpotent Jordan blocks are a little more delicate: First, arrange them in decreasing order by size. Then in order for a square root to exist, it is necessary and sufficient that the (2n1)st block is at most one larger than the (2n)th for all positive n. Approximate proof: If you square a nilpotent block it becomes two blocks of approximately equal size.

$\begingroup$ I presume sbseminar.wordpress.com/2008/12/02/… is relevant? $\endgroup$ – Yemon Choi Feb 4 '10 at 4:49

One does not actually need the Jordan form to construct a square root in the case of nonzero eigenvalues, although simply knowing that it exists is quite helpful. Let $a_1$, ..., $a_n$ denote the nonzero complex eigenvalues. Then a square root $S$ is given by $$S=P(A),$$ where $P$ is any polynomial such that $$(d/dz)^m P(z) = (d/dz)^m \sqrt z$$ for all $z=a_k$ and $0 \le m\le d1$, where $d$ is the dimension of the space. (Here we have used complex derivatives, and you can choose any branch of the complex square root you like. One may find higherorder roots or compute functions such as $\exp (A)$ similarly.)
This follows from the Jordan form, which tells us that for any polynomial $Q$ the value of $Q(A)$ depends only on the values of $Q$ and its derivatives up to order $b_k1\le d1$ at each of the eigenvalues of $A$, where $b_k$ is the size of the maximal Jordan block for the corresponding eigenvalue $a_k$. For example $$Q\left( \left[ \begin{array} [c]{cccc}% \lambda & 1 & & \\ & \lambda & 1 & \\ & & \lambda & 1\\ & & & \lambda \end{array} \right] \right) =\left[ \begin{array} [c]{cccc}% Q\left( \lambda\right) & Q^{\prime}\left( \lambda\right) & \frac{1} {2!}Q^{\prime\prime}\left( \lambda\right) & \frac{1}{3!}Q^{\prime \prime\prime}\left( \lambda\right) \\ & Q\left( \lambda\right) & Q^{\prime}\left( \lambda\right) & \frac {1}{2!}Q^{\prime\prime}\left( \lambda\right) \\ & & Q\left( \lambda\right) & Q^{\prime}\left( \lambda\right) \\ & & & Q\left( \lambda\right) \end{array} \right]$$ In particular, $$S^2=(P(A))^2=(P^2)(A)=A,$$ by choosing $$Q=P^2.$$ [By construction, $Q(z)=z$ to order $d1$ at the eigenvalues of A.]
In the infinitedimensional case one has the holomorphic calculus (also called the Dunford calculus), although in the finitedimensional case one needs only polynomials. Indeed, in the finitedimensional case the Dunford calculus simply extends the identity above for applying a polynomial $Q$ to a $4 \times 4$ Jordan block to the case that $Q$ is a holomorphic function.
This is not an exhaustive solution even in the case of nonzero eigenvalues, since it only finds square roots which are polynomial functions of the original matrix. For example, all reflections are square roots of the identity.

$\begingroup$ Just FYI: each repeated edit bumps a post to the top of the stack, and this sometimes causes an annoyance for others. With more than 10 edits since posting, perhaps the post has reached a stable condition? $\endgroup$ – Todd Trimble♦ May 27 '14 at 20:34

$\begingroup$ Oops, sorry! It would be good if an "edit without bumping" button were included. $\endgroup$ – J Tyson May 27 '14 at 23:49

$\begingroup$ That's alright. Yes, I wish there were something like that. $\endgroup$ – Todd Trimble♦ May 27 '14 at 23:53
I asked a question on similar topic a while ago. The answer I got is using newton's method.
My original question was ""Approximate the square root of (1X) efficiently through (nested) products However, I think the method applied to your problem.
Here is an reference: Newton's Method for the Matrix Square Root
There are also papers for pth root and inverse pth root: A SchurNewton method for the matrix pth root
The general idea is that 1) we need to scale your matrix, so that its eigenvalues fall into a specific circle. 2)then, apply Newton's method with initial value $y_0=I$.
The update rule for pth root is
$$ y_{k+1}=\frac{1}{p}[(p1)y_k+y_k^{1p}A] $$
Then we will have $y_k^{p}A\to I$
For the record, since this question was recently resurrected: computing a Jordan form is not recommendable in numerical (i.e., not symbolical) computation. The standard algorithm for implementation on a computer in floatingpoint arithmetic is the Schur method (described very quickly: compute a Schur form, then get the entries one superdiagonal at a time via a sort of backsubstitution). You will find everything you need in Chapter 6 of N. Higham's book Functions of matrices.
Another caveat, already mentioned in the comments, is that there are in general multiple square roots (think about the diagonal case with $2^n$ choices of signs, or the $2\times 2$ solutions of $X^2=I$ which are infinite). This is discussed at large in Chapter 1 of the same book.
As usual, you will most often get better results by sticking with the default implementations for your programming language of choice.
There are some theoretically applicable methods here. In particular, you can pick a large enough circle that contains all the eigenvalues of your matrix and apply Cauchy Integral Formula as described in the above article. This method does not require the matrix spectrum, only a bound on its eigenvalues.