Matrix-Square Root 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.    
 A: 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.
A: 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.
A: 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."
A: 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$.
A: 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... 
A: 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.
