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}}$.



Square root of a large sparse symmetric positive definite matrix
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> 


