8
$\begingroup$

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

Improve this question
$\endgroup$
14
  • 1
    $\begingroup$ 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. $\endgroup$
    – fedja
    Commented May 27, 2018 at 16:59
  • 1
    $\begingroup$ @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. $\endgroup$
    – fedja
    Commented May 27, 2018 at 19:48
  • 1
    $\begingroup$ 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...). $\endgroup$
    – fedja
    Commented May 27, 2018 at 19:56
  • 1
    $\begingroup$ @BillBradley I mean the unique PSD root. $\endgroup$
    – messcode
    Commented May 28, 2018 at 5:26
  • 1
    $\begingroup$ @BrendanMcKay Are you sure it'll work for non-commuting $B,D$? I don't even see why your expression gives a symmetric matrix... $\endgroup$
    – fedja
    Commented May 28, 2018 at 15:51

1 Answer 1

Reset to default
4
$\begingroup$

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

Improve this answer
$\endgroup$
2
  • 2
    $\begingroup$ There's no need to factor the polynomial, just evaluate $P(A)X$ by left-multiplying $X$ by $A$ repeatedly. $\endgroup$
    – Federico Poloni
    Commented May 28, 2018 at 17:34
  • $\begingroup$ @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$". $\endgroup$
    – fedja
    Commented May 29, 2018 at 10:58

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .