Applying matrix square root inverse in matrix-free regime

Question

Let $A$ be a large symmetric positive definite matrix, and suppose that we can efficiently apply $A$ and have a fast solver to apply $A^{-1}$, but we do not have access to the matrix entries for either $A$ or $A^{-1}$.

How can one use this to solve the system, $$\sqrt{A}x=b?$$

---

As suggested in this previous question, the problem can be solved by computing an eigenvalue decomposition of $A$. In my case a full eigenvalue decomposition is infeasible for speed and memory reasons, as are factorizations that focus on elements of the matrix such as Cholesky.

Answer 1

Nick, your answer is certainly valid. However, the additional solve using $Ax = \tilde{b}$ where $\tilde{b} = A^{1/2}b$ can be avoided if you directly solve $f(A)x = b$. This is the subject of a book chapter by Henk van der Vorst http://link.springer.com/chapter/10.1007%2F978-3-642-58333-9_2.

Several other options for computing matrix-square root (or its inverse) are available. As you mentioned solving an ODE is one such approach. Broadly speaking, there are a few other approaches - polynomial approximation to square root http://link.springer.com/article/10.1007%2FBF02083211#page-1, rational approximations to square root http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.131.1310, Krylov subspace based approaches http://epubs.siam.org/doi/abs/10.1137/130920587 and Extended Krylov subspace method http://epubs.siam.org/doi/pdf/10.1137/S0895479895292400.

I have only provided a single example for each case but the literature on computing matrix square root is rich. The particular choice of the method you use will ultimately depend on the matrix and the distribution of eigenvalues. However, I tend to favor the rational approximation based on contour integral http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.131.1310 because even for ill-conditioned matrix, computing the f(A)b boils down to computing about a dozen solves on a shifted system of equations. It may not work if A is only available in matrix-free form, in which case the Extended Krylov method may be your best bet.

Answer 2

After looking through some literature I'm answering my own question. I'll wait on accepting for a few days in case anyone else has a better answer.

First of all, the problem is equivalent to solving, $$A x = \sqrt{A}b,$$ so the only square root part needed is computing $\sqrt{A}b$, after which it is a normal solve.

The following paper discusses two ways to do that,

Numerical approximation of the product of the square root of a matrix with a vector, Allen, EJ and Baglama, J and Boyd, SK , Linear Algebra and its Applications 2000, http://www.sciencedirect.com/science/article/pii/S0024379500000689

The second method transforms the problem into solving an ODE whose solution at time $t=1$ is $A^{1/2}b$. This ODE is, $$\begin{cases} \frac{dx}{dt} = -\frac{1}{2}\left(At + (1-t)I\right)^{-1}(I-A)x(t) \\ x(0) = b, \end{cases}$$ which can be solved through standard techniques such as forward Euler. At each timestep $t_k$, $(I-A)$ must be applied to $x(t_k)$, and then a regularized version of $A^{-1}$: $$\left(At_k + (1-t_k)I\right)^{-1}$$ must be applied to the result. In my particular situation the method for applying $A^{-1}$ can be extended to also apply $(A + \alpha I)^{-1}$. In other situations it should be possible to use $A^{-1}$ to build a preconditioner for $\left(At_k + (1-t_k)I\right)$.

Edit: A major downside of this method is that the ODE becomes stiff as the condition number of $A$ increases.