# Square root of a large sparse symmetric positive definite matrix

I am trying to calculate

$$Y = A^{\frac 12} X$$

where $A$ is a very large and sparse positive definite matrix, say, $10^4 \times 10^4$. Matrix $X$ is known and, say, $10^4 \times 100$. Is there any method that can compute or approximate $A^{\frac{1}{2}}$ efficiently?

I noticed that matrix $A$ can be written as $D + B$, where $D=\mbox{diag}(A)$ and $B=A-D$, which is very sparse. But I still do not know how to compute the square root of the sum of matrices, $(D + B)^{\frac{1}{2}}$.

• If you have an a priori information about the condition number $K$ of $A$, you can just approximate the square root by a polynomial on $[K^{-1},1]$ with the desired precision. You can do it even on the entire interval $[0,1]$, but the convergence is then too slow to be of practical use. – fedja May 27 '18 at 16:59
• @Mahdi I just mean that if all eigenvalues are between $K^{-1}$ and $1$ and $|P_n(x)-\sqrt x|<\delta$ there, then $\|P(A)-A^{1/2}\|\le\delta$ as well Since the computation of $P(A)x$ requires just $deg(P)$ multiplications of a vector by $A$, this can be of use if we can make $\delta$ small enough without making $deg(P)$ too large. – fedja May 27 '18 at 19:48
• Small sanity check: how do you know that $A$ is positive definite? (I mean I have once been asked to prove that a certain crazy matrix with complicated expressions for entries was positive definite, spent three days representing it as a Gram matrix, and then was told that "but yes, that is how we computed it in the first place!". Perhaps you also have some background objects that can be used directly and not read from the matrix in complicated ways...). – fedja May 27 '18 at 19:56
• @BillBradley I mean the unique PSD root. – messcode May 28 '18 at 5:26
• @BrendanMcKay Are you sure it'll work for non-commuting $B,D$? I don't even see why your expression gives a symmetric matrix... – fedja May 28 '18 at 15:51

I completely agree with fedja: there is a nice method here (which, unfortunately, does not always work well). If you know bounds for the spectrum of $A$, say $0<a<\lambda<b$, then you (sometimes) can compute the square root very efficiently by approximating the square root by a polynomial on this interval $$\sqrt{x}\approx P(x),\,a<x<b$$ and then simply computing $P(A)$. By the way, it is not necessary (and not even reasonable) to keep the intermediate matrix $P(A)$ in the computer memory; the simplest way to avoid this is to factor the polynomial.
• There's no need to factor the polynomial, just evaluate $P(A)X$ by left-multiplying $X$ by $A$ repeatedly. – Federico Poloni May 28 '18 at 17:34
• @FedericoPoloni Exactly. I guess you meant "Do it in $\deg(P)$ multiplications by the standard scheme $aA^3x+bA^2x+cAx+dx=A(A(Ax+bx)+cx)+dx$". – fedja May 29 '18 at 10:58