Question

What methods exist for finding the square-root of a non-diagonalizabe positive complex matrix?

Answer 1

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=B-I$ 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 $(m-1)/2$ by $(m-1)/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 $(m-1)/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=A-tI$ and the formula is valid for $t$ greater then the maximum eigenvalue of $A$.

Answer 2

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 (2n-1)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.