Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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.

share|improve this question
My first attempt would be trying rational Krylov methods; see for instance – Federico Poloni Apr 15 '14 at 10:18
up vote 6 down vote accepted

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

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, rational approximations to square root, Krylov subspace based approaches and Extended Krylov subspace method

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

share|improve this answer
Thank you, perfect. This is the best answer. I also ran across the Hale Higham and Trefethen contour integral paper yesterday and think that either that or rational approximation will be the best option for me. – Nick Alger Apr 16 '14 at 20:56

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,

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.

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.