Children in primary school learn that every positive number has a real square root. The number x is a square root of s, if x2 = s.
Did you know that matrices can also have square roots? For certain matrices S, you can find another matrix X such that X*X = S. To give a very simple example, if S = a*I is a multiple of the identity matrix and a > 0, then X = ±sqrt(a)*I is a square root matrix.
Did you know that some matrices have square roots? Click To TweetA matrix square root
I'm going to restrict this article to real numbers, so from now on when I say "a number" I mean a real number and when I say "a matrix" I mean a real matrix. Negative numbers do not have square roots, so it is not surprising that not every matrix has a square root. For example, the negative identity matrix (-I) does not have a square root matrix.
All positive numbers have square roots, and mathematicians, who love to generalize everything, have defined a class of matrices with properties that are reminiscent of positive numbers. They are called positive definite matrices, and they arise often in statistics because every covariance and correlation matrix is symmetric and positive definite (SPD).
It turns out that if S is a symmetric positive definite matrix, then there exists a unique SPD matrix X, (called the square root matrix) such that X2 = A. For a proof, see Golub and van Loan, 3rd edition, 1996, p. 149. Furthermore, the following iterative algorithm converges quadratically to the square root matrix (ibid, p. 571):
- Let X0=I be the first approximation.
- Apply the formula X1 = (X0 + S*X0-1) / 2 where X-1 is the inverse matrix of X. The matrix X1 is a better approximation to sqrt(S) than X0.
- Apply the formula Xn+1 = (Xn + S*Xn-1) / 2 iteratively until the process converges. Convergence is achieved when a matrix norm ∥ Xn+1 – Xn ∥ is as small as you like.
The astute reader will recognize this algorithm as the matrix version of the Babylonian method for computing a square root. As I explained last week, this iterative method implements Newton's method to find the roots of the (matrix) function f(X) = X*X - S.
Compute a matrix square root in SAS
In SAS, the SAS/IML matrix language is used to carry out matrix computations. To illustrate the square root algorithm, let S be the 7x7 Toeplitz matrix that is generated by the vector {4 3 2 1 0 -1 -2}. I have previously shown that this Toeplitz matrix (and others of this general form) are SPD. The following SAS/IML program implements the iterative procedure for computing the square root of an SPD matrix:
proc iml;
/* Given an SPD matrix S, this function to compute the square root matrix
X such that X*X = S */
start sqrtm(S, maxIters=100, epsilon=1e-6);
X = I(nrow(S)); /* initial starting matrix */
do iter = 1 to maxIters while( norm(X*X - S) > epsilon );
X = 0.5*(X + S*inv(X)); /* Newton's method converges to square root of S */
end;
if norm(X*X - S) <= epsilon then return( X );
else return(.);
finish;
S = toeplitz(4:-2); /* 7x7 SPD example */
X = sqrtm(S); /* square root matrix */
print X[L="sqrtm(S)" format=7.4]; |
The output shows the square root matrix. If you multiply this matrix with itself, you get the original Toeplitz matrix. Notice that the original matrix and the square root matrix can contain negative elements, which shows that "positive definite" is different from "has all positive entries."
Functions of matrices
The square root algorithm can be thought of as a mapping that takes an SPD matrix and produces the square root matrix. Therefore it is an example of a function of a matrix. Nick Higham (2008) has written a book Functions of Matrices: Theory and Computation, which contains many other functions of matrices (exp(), log(), cos(), sin(),....) and how to compute the functions accurately. SAS/IML supports the EXPMATTRIX function, which computes the exponential function of a matrix.
The square root algorithm in this post is a simplified version of a more robust algorithm that has better numerical properties. Higham (1986), "Newton's Method for the Matrix Square Root" cautions that Newton's iteration can be numerically unstable for certain matrices. (For details, see p. 543, Eqns 3.11 and 3.12.) Higham suggests an alternate (but similar) routine (p. 544) that is only slightly more expensive but has improved stability properties.
I think it is very cool that the ancient Babylonian algorithm for computing square roots of numbers can be generalized to compute the square root of a matrix. However, notice that there is an interesting difference. In the Babylonian algorithm, you are permitted to choose any positive number to begin the square root algorithm. For matrices, the initial guess must commute with the matrix S (Higham, 1986, p. 540). Therefore a multiple of the identity matrix is the safest choice for an initial guess.
13 Comments
Very cool indeed! I love how your blog posts are individual gems with the bonus of some of them being a chapter to an encompassing blog post. A reward for being "A Do Loop" follower! Thanks...
Sir, I really thankful to you. Its very much helpful for me.I follow this site more than one and half years.this site help me to grow my knowledge related to SAS and SAS SQL procedure
/* Given an SPD matrix S, this function to compute the square root matrix
33 X such that X*X = S */
34 start sqrtm(S, maxIters=100, epsilon=1e-6);
_
22
200
ERROR 22-322: Syntax error, expecting one of the following: ), ','.
ERROR 200-322: The symbol is not recognized and will be ignored.
35 X = I(nrow(S));
35 ! /* initial starting matrix */
36 do iter = 1 to maxIters while( norm(X*X - S) > epsilon );
37 X = 0.5*(X + S*inv(X));
37 ! /* Newton's method converges to square root of S */
38 end;
39 if norm(X*X - S) <= epsilon then return X;
_
22
202
ERROR 22-322: Syntax error, expecting one of the following: ;, (.
ERROR 202-322: The option or parameter is not recognized and will be ignored.
The error is on a line before the START statement. You've probably forgotten a semicolon or a set of parentheses.
He simply runs your code, which fails, and you tell him that he forgot something? Sounds more like you forgot something.
Instead of
return X;
it must be
return(X);
at least this worked for me. Might be due to SAS Version?
Thanks for writing. My comment was about the FIRST error in his program. You are correct to point out that the log shows a second error near the end of the program.
The RETURN function has supported optional parentheses since SAS/IML 13.1, which was released in 2013 as part of SAS 9.4M1. If you are using an older version of SAS, you need to use parentheses around the argument. I will update the post so that no one else is confused.
Theoretically yes. In practice the numerical errors propagate like crazy because once you have a slight non-commutativity between $A$ and $X$ you are ending with complete nonsense in 20 steps or fewer instead of getting the result with machine precision. Try the method on A={{1,2},{2,4.01}} and you'll see what I mean...
Yes, I mentioned this fact in the second-to-last paragraph and included a reference to Higham (1986), which has better stability.
How can we do square root of a matrix in SAS , without using proc iml , Sir?
You can perform the same computation in another language that supports matrix computations. For example, you could use MATLAB, R, or the numPy package in Python.
Thank you sir for your reply , however i am restricted to use SAS without an IML option.
Running the IML code I obtained a different output instead of the matrix you show in this blog post.
Thank you! I have updated the post to display the correct matrix.
The matrix you saw was the square-root matrix for the 6x6 Toeplitz matrix generated by {6,5,4,3,2,1}. I later updated the program to use the 7x7 Toeplitz matrix generated by {4,3,2,1,0,-1,2}. However, I forgot to update the image. The square-root matrix for the first matrix has all positive entries. The square-root matrix for the second matrix has some negative entries, which is why I decided to switch examples.