I was wondering about matrix square roots. What is the procedure to evaluate $(X^{T}X)^{-1/2}$? Is it by a spectral decomposition of $(X^{T}X)^{-1}$ as $U\lambda U^{T}$ followed by the square root $S$ obtained as $U\lambda^{1/2}$? I would like to hear more.
$\begingroup$
$\endgroup$
Add a comment
|
6 Answers
$\begingroup$
$\endgroup$
1
The best way I can think of to calculate this particular matrix square root is to calculate a singular value decomposition of $X$. You should be able to use the LAPACK routine xgesvd, where 'x' is the appropriate choice of s (single), d (double), or z (complex).
Suppose that $X = U\Sigma V^{H}$ is a singular value decomposition of $X$, where the superscript $H$ stands for Hermitian (complex conjugate) transpose, $U$ and $V$ are unitary, and $\Sigma$ is diagonal with nonnegative elements.
Then $X^{H}X$ = $V\Sigma^{2}V^{H}$ is an eigenvalue decomposition. It follows that
$$(X^{H}X)^{-1/2} = V\Sigma^{-1}V^{H}$$
exists iff all of the singular values of $X$ are strictly positive.
-
$\begingroup$ "all of the singular values of $X$ are strictly positive." - or equivalently, $X$ has full column rank. I think this is the best approach, since the mere formation of the cross-product matrix is a rather ill-conditioned route to take... $\endgroup$ Feb 9, 2012 at 0:10
$\begingroup$
$\endgroup$
I'm not sure whether there's a canonical "the" procedure, but diagonalization certainly works, assuming that $X$ has real entries. Since $X^TX$ is real symmetic, it is diagonalizable. Diagonalize it and raise every eigenvalue to the $-\frac{1}{2}$ power (which is always possible because there can be no negative eigenvalues).
There's no need to invert the matrix before you diagonalize it, since it is so much easier just to invert the diagonal elements afterwards.
answered Feb 7, 2012 at 20:30
hmakholm left over Monicahmakholm left over Monica
282k23 gold badges418 silver badges678 bronze badges
$\begingroup$
$\endgroup$
3
Another method is to use the polar decomposition $X = UH$ with $U$ unitary and $H$ symmetric positive definite. Then $A^{1/2} = H$ because $H^2 = H U^T U H = X^T X = A$, so you just have to invert $H$. This can be done using the Cholesky decomposition: if $H = R^T R$ then $H^{-1} = R^{-1} R^{-T}$, and triangular matrices can be easily inverted.
I'd guess that this method is more stable than diagonalization (because you don't form $A$) but probably a bit slower. The SVD is probably even slower.
As noted in the comments, one method for computing the polar decomposition is to first compute the SVD, and if you do that you can just as well use the SVD directly to compute the inverse square root. However, there are other methods. I have no experience with this, but there is quite some discussion in Chapter 8 of Higham, Functions of Matrices, SIAM, 2008. One method is to use the Newton iteration for the equation $U^TU=I$, which is $X_{k+1} = \frac12 (X_k + X_k^{-T})$ with initial condition $X_0 = A$; this converges quadratically to the polar factor $U$. Better is to add scaling: $X_{k+1} = \frac12 (\mu_k X_k + \mu_k^{-1} X_k^{-T})$. Higham concludes that the resulting algorithm (Algorithm 8.20 in his book) is "remarkably stable, quick to converge, and robust […] And its flop count is less than that for computation of the polar factors for the SVD."
-
$\begingroup$ I would however note that a number of algorithms for the polar decomposition rely on SVD, so one might not be able to avoid it entirely... $\endgroup$ Feb 9, 2012 at 0:08
-
$\begingroup$ I thought about the polar decomposition method, too, when I looked at the problem, but then I thought about calculating it, looked it up, and came to a similar conclusion as J.M.. $\endgroup$ Feb 9, 2012 at 7:09
-
$\begingroup$ I added a bit on computing the polar decomposition. $\endgroup$ Feb 9, 2012 at 9:48
$\begingroup$
$\endgroup$
Hmm, a little less thinking (but more computerpower) requires the Newton-iteration. If $\small Z=X^t X $ then with $\small A_0=I $ and then iterations $\small A_{k+1} = (Z \cdot A_k^{-1}+A_k)/2 $ and $\displaystyle\small \lim_{k \to \infty} A_k^{-1} = Z^{ - \frac12} $ .
Maybe we need some more restrictions on the eigenvalues of Z to get it converging (at least we have numerical issues due to inversions of the iterations of A), but some random examples worked well, so it is at least one of the list of possible methods...
answered Feb 15, 2012 at 16:37
$\begingroup$
$\endgroup$
4
One way to deal with this problem is to approximate the solution by iterative methods. For example, by using Newton's method to solve the matrix equation $X^2-Y=0$. It leads to the following iteration: $$X_{k+1}=\frac{1}{2}(X_k+X_k^{-1}Y)$$ where $x_0=I$. Note that there is no need to explicitly compute $X_k^{-1}$. Rather, it's possible to solve the least-squares problem $\min_z\|Xy-z\|_2$, for every column $z$ of the matrix $X_k^{-1}Y$.
-
$\begingroup$ "...there is no need to explicitly compute $X_k^{-1}Y$..." - your wording's a tad confusing. Did you mean that you can form $X_k^{-1}Y$ a column at a time, add to the corresponding column of $X_k$, and then halve? $\endgroup$ Feb 15, 2012 at 15:56
-
$\begingroup$ I meant computing $X_k^{-1}$. The reason I'd like to get around this operation is that it is both computationally expensive and numerically unstable. $\endgroup$ Feb 15, 2012 at 17:58
-
$\begingroup$ Ah. Then I certainly agree. Though, why wouldn't (pivoted) Gaussian elimination be used here? $\endgroup$ Feb 15, 2012 at 18:17
-
$\begingroup$ Gaussian elimination costs $O(n^3)$ time while an iteration for an iterative LS solution method such as GMRES costs only $O(n^2)$. $\endgroup$ Feb 21, 2012 at 10:01
$\begingroup$
$\endgroup$
2
Will any matrix r with r*r' = inv(X'*X) do? If so I think an easy one to calculate is the inverse of the upper triangle from a QR decomposition of X.
-
1$\begingroup$ You are in fact computing a Cholesky decomposition, not a square root... $\endgroup$ Feb 9, 2012 at 11:50
-
$\begingroup$ Yes, but then the S in the original post was't a square root either. $\endgroup$– dmuirFeb 9, 2012 at 12:35