Finding the Square-Root of a Non-diagonalizable Positive Matrix

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.

Answer 3

One does not actually need the Jordan form to construct a square root in the case of non-zero eigenvalues, although simply knowing that it exists is quite helpful. Let $a_1$, ..., $a_n$ denote the non-zero 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 d-1$, 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 higher-order 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_k-1\le d-1$ 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 $d-1$ at the eigenvalues of A.]

In the infinite-dimensional case one has the holomorphic calculus (also called the Dunford calculus), although in the finite-dimensional case one needs only polynomials. Indeed, in the finite-dimensional 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.

Answer 4

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 (1-X) 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 p-th root and inverse p-th root: A Schur-Newton 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 p-th root is

$$ y_{k+1}=\frac{1}{p}[(p-1)y_k+y_k^{1-p}A] $$

Then we will have $y_k^{-p}A\to I$

Answer 5

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 floating-point 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 back-substitution). 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.

Answer 6

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.