# 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. Commented May 27, 2018 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. Commented May 27, 2018 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...). Commented May 27, 2018 at 19:56
• @BillBradley I mean the unique PSD root. Commented May 28, 2018 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... Commented May 28, 2018 at 15:51

## 1 Answer

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.

The method was actually considered by many: https://link.springer.com/article/10.1007/bf02083211 , https://epubs.siam.org/doi/abs/10.1137/S0895479895292400 . Still better approximation may be achieved using rational functions (Optimal Finite Difference Grids and Rational Approximations of the Square Root I. Elliptic Problems. Ingerman, Druskin, Knizhnerman, 2000.)

There are also some alternatives, like Pade approximation https://epubs.siam.org/doi/abs/10.1137/S089547989731631X or Schur factorization https://www.sciencedirect.com/science/article/pii/002437958380010X

• There's no need to factor the polynomial, just evaluate $P(A)X$ by left-multiplying $X$ by $A$ repeatedly. Commented May 28, 2018 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$". Commented May 29, 2018 at 10:58