# Accurate Way to Calculate Matrix Powers and Matrix Exponential for Sparse Positive Semidefinite Matrices

I do need to numerically calculate the following forms for any $$x\in\mathbb{R}^n$$, possibly in python:

1. $$x^T M^k x$$, where $$M\in\mathbb{R^{n\times n}}$$ is a PSD sparse matrix, $$n$$ can be quite large in the orders of $$10^3$$ or higher, and where $$k$$ can get large values $$-$$possibly up to order of hundreds. I prefer $$k$$ to be a real number, but it is ok if it can only be an integer, if that makes a considerable accuracy difference.
2. Similarly I am interested in $$x^Te^{-tM}x$$ where $$t$$ is a non-positive real value.

For case 1 I can either:

• Use the scipy.linalg.fractional_matrix_power to calculate $$M^k$$ and then derive $$x^TM^Kx$$, or
• Use scipy.linalg.svd to find SVD of $$M$$ as $$U\Lambda U^T$$ and then finally evaluate the desired value using $$x^TU\Lambda^k U^Tx$$. According to Cleve Moler et. al's paper "Nineteen dubious..." (method 6) SVD shouldn't be a bad idea, but I am still not sure for my overall calculation which is $$x^T MX$$ this is the best way, and also people here discouraged me.
• Finally if $$k$$ is integer, and again based on SVD I can calculate $$\|x^TM^{k/2}\|^2$$.

For case 2

• Again I can use off the shelf scipy.linalg.expm
• I can do SVD for $$M$$ and then go with $$x^TU\exp(\Lambda) U^Tx$$.
• Finally since I am only interested in $$x^T \exp(M) x$$, and not exactly $$\exp(M)$$ it self, I can consider the Taylor expansion of $$x^T{\rm expm}(M)x\approx \sum_{i=0}^{l} \frac{1}{i!}x^TM^ix$$ for some $$l$$ that controls the precision, and $$x^TMx$$ can be calculated based on case 1.

Can anybody guide me about what is the most precise way to calculate either of these expressions, up to hopefully machine precision? Any of these methods, or they're better solutions out there? I would be happy also with references.

P.S. I have cross-posted this question on math.stackexchange.com, mathoverflow, and stackoverflow and here is suggested as the write place to post it.

UPDATE I initially accepted the answer by tch but running their code --depending on the input-- does produce unstable results. So I wanted to share the experience and see if I am doing something wrong, or there is problem in the solution they proposed that I am getting unstable results:

When I try the following

n = 3000
U,_ = np.linalg.qr(np.random.randn(n,n))
U = U.astype(np.float64)
lam = np.exp(np.abs(np.linspace(-5,5,n).astype(np.float64) + np.random.randn(n)*0.4))
print (lam)
M = U@np.diag(lam)@U.T.astype(np.float64)
x = np.random.randn(n).astype(np.float64)


And then loop over the answers from k=10 to k=49:

for i in range(10, 50):
_,(a_,b_) = lanczos(M,x,i)


I sometimes get meaningless results:

Can somebody shed a light on this please?

UPDATE 2

The problem might have been from how we import scipy, from what I can see. I am not able to replicate the issue anymore after explicitly importing scipy.linalg instead of all of scipy. For now I will accept tch's answer. Thanks a lot!

• Directed here from mathoverflow.net/q/384941/1898 – Federico Poloni Feb 25 at 20:32
• Thank you. Also from math.stackexchange.com/questions/4039692 – Cupitor Feb 25 at 20:37
• In my opinion the problem size (N=1000) doesn't seem quite big enough to pursue any exploitation of sparsity, especially since the best algorithms in this space seem to all be based on some sort of spectral decompositions, not LU-like ones. I see the mention of (semi) positivity .. is it also symmetric? – rchilton1980 Feb 25 at 21:00
• Any reason why you use svd and not eigh? – Federico Poloni Feb 25 at 21:10
• Yes, an eigenvalue algorithm can exploit the symmetry, you should use it instead. – rchilton1980 Feb 25 at 22:05

This is a slightly modified version of my response on math.stackexchange. One standard approach to computing matrix functions times a vector $$f(M)x$$ or quadratic forms $$x^Tf(M)x$$ when $$M$$ is symmetric is via the Lanczos algorithm. Lanczos computes an orthonormal basis $$Q_k = [q_1, \ldots, q_k]$$ for Krylov subspace $$\operatorname{span}(x,Ax,\ldots, A^{k-1}x)$$ by a Gram-Schmidt like procedure. This results in a factorization $$AQ_k =Q_kT_k + \beta_k q_{k+1}e_k^T$$, where $$e_k$$ is the $$k$$-th canonical basis vector and $$T_k$$ is a $$k\times k$$ symmetric tridiagonal matrix.

The "Lanczos approximation" to $$f(M)x$$ is defined as $$Q_k f(T_k) Q_k^T x$$, and the approximation to $$x^Tf(M)x$$ is defined as $$x^TQ_k f(T_k) Q_k^T x$$. Note that $$f$$ can be any scalar function, for instance $$f(x) = \exp(x)$$ or $$f(x) = x^i$$. There are software packages to compute $$Q_k \exp(T_k) Q_k^T x$$. This can easily be turned into the Lanczos approximation to $$x^T\exp(M)x$$ by taking the inner product with $$x$$. The Lanczos approximation to the quadratic form can be viewed as a certain Gaussian quadrature approximation to $$x^T f(M) x$$ and converges extremely quickly in $$k$$ if $$f$$ is analytic (see Golub Meurant "Matrices Moments and Quadrature").

This approach has runtime $$O(T_{mv}k + nk)$$ (or $$O(T_{mv}k + nk^2)$$ with full reorthogonalization), where $$T_{mv}$$ is the cost of evaluating $$v\mapsto Mv$$. This will almost certainly be cheaper than an SVD or any $$n^3$$ algorithm since $$k$$ can typically be taken to be fairly small even for a machine precision accurate output.

Note that if $$x$$ is unit length $$Q_k^Tx = e_1$$. This means the Lanczos approximation to $$x^Tf(M)x$$ is given by $$e_1^T f(T_k) e_1$$ and that $$Q_k$$ doesn't need to be stored. There are also some subtleties about whether full reorthogonalization should be used or not. Whether or not it is used, the output will still be accurate for the matrix exponential, but without reorthogonalization $$k$$ may need to be larger. However, reorthogonalizaton is more computationally expensive, so there is some tradeoff.

Lanczos and a function to compute the approximation of the quadratic form are easy to implement. Note that if you do not want to use reorthogonalization, you could rewrite the code to avoid storing Q.

def lanczos(A,q0,k,reorth=True):

n = len(q0)

Q = np.zeros((n,k),dtype=A.dtype)
a = np.zeros(k,dtype=A.dtype)
b = np.zeros(k-1,dtype=A.dtype)

Q[:,0] = q0 / np.sqrt(q0.T@q0)

for i in range(1,k+1):
# expand Krylov space
qi = A@Q[:,i-1] - b[i-2]*Q[:,i-2] if i>1 else A@Q[:,i-1]

a[i-1] = qi.T@Q[:,i-1]
qi -= a[i-1]*Q[:,i-1]

if reorth:
qi -= Q@(Q.T@qi) # regular GS

if i < k:
b[i-1] = np.sqrt(qi.T@qi)
Q[:,i] = qi / b[i-1]

return Q,(a,b)

theta,S = sp.linalg.eigh_tridiagonal(a_,b_,tol=1e-30)

return norm_b**2*S[0]@(f(theta)*S[0])


To use the code you simply do something like:

_,(a_,b_) = lanczos(M,x,100)


Even for a matrix of $$n=3000$$, its pretty easy to see that this approach is much faster than SVD or eigenvalue decomposition. For instance, we can generate a matrix with the specificed eigenvalues:

n = 3000
U,_ = np.linalg.qr(np.random.randn(n,n))
lam = np.linspace(0,1,n)
M = U@np.diag(lam)@U.T
x = np.random.randn(n)


For $$k=7$$, the Lanczos approximation obtains full machine precision for the matrix exponential in 50ms on my machine:

xexpMx = (x.T@U)@np.diag(np.exp(lam))@(U.T@x)
_,(a_,b_) = lanczos(M,x,7)


On the other hand, even computing the (symmetric) eigendecomposition of M takes 2-3 seconds.

• Sorry @tch to take my accept back. I just tried to play around with it and couldn't get stability. Could you please point out what I might be doing wrong that this happens? – Cupitor Feb 28 at 6:13
• what do you mean "couldn't get stability"? Hard to say without more details – tch Feb 28 at 17:11
• I meant my attempt of running your code. Did you see my update section in the original post? I have the little snippet that you can run. In the plot you can see the relative error goest to -3e72 – Cupitor Mar 1 at 7:37
• When I run your code it works fine. You should post a full MWE for the code, otherwise it's impossible to determine the cause. Perhaps the RHS you used for lanczos was different than for computing the true solution? Also, as a note, the first entries of a_, and b_ are the same regardless of $k$, so you can run Lanczos once outside the loop and then just take the first i and i-1 entries of a_ and b_ respectively. – tch Mar 1 at 16:35
• Also thanks a lot for the headsup about running Lanczos once. I guess I was able to understand I don't need to calculate Lanczos everytime. I was just trying to do a sanity check on the code, before trying to change it in any way. – Cupitor Mar 1 at 21:02

For symmetric, more generally, normal, matrices there is no danger in diagonalizing a matrix to evaluate a matrix function, because all eigenvectors are perfectly conditioned. You should get answers with errors close to the machine epsilon times the condition number of the function evaluation problem. This is as accurate as you can get in floating point. Here's a matlab program that shows that the methods based on diagonalization yield answers very close to those of the standard libraries. I assume you would get the same in python, but I think I remember reading of some problem with one of these functions in python.

n = 1000;
% create a random PSD matrix and random x, k, and t
U = orth(randn(n));
S = diag(rand(n, 1));
M = U*S*U';
x = rand(n, 1);
k = randn(1);
t = -rand(1);
% try different methods for problem matrix power
v1 = x'*M^k*x
v2 = (x'*U)*diag(diag(S).^k)*(U'*x)
[U1, S1] = eig(M);
v3 = (x'*U1)*diag(diag(S1).^k)*(U1'*x)

% try different methods for expm
v4 = x'*expm(-t*M)*x
v5 = (x'*U1)*diag(exp(-t*diag(S1)))*(U1'*x)


Here's what I got for a sample run:

v1 = 1294.749371499737
v2 = 1294.749371499742
v3 = 1294.749371499842
v4 = 464.4679977927723
v5 = 464.4679977925199


If you want to work which much larger matrices, you do need to make use of the fact that the matrix is sparse and you only need $$f(M)x$$. Several suitable methods are described here: https://www.maths.manchester.ac.uk/~higham/narep/narep436.pdf

and a more recent paper, also by Nick Higham, is: http://eprints.ma.man.ac.uk/1591/1/covered/MIMS_ep2010_30.pdf

Matlab code for this paper is here: https://github.com/higham/expmv/blob/master/expmv.m

and if you want to evaluate $$x'e^{-tM}x$$ for many $$t$$ values, there is code here: https://github.com/higham/expmv/blob/master/expmv_tspan.m

• Thanks so much Amit for your thorough reply! I wish I could accept more than one answer. I will consider finding an equivalent of the matlab codes, or just go ahead with your suggestion of diagonalizing! – Cupitor Feb 27 at 19:15

$$M = \sum_i\lambda_iv_iv_i^\top$$

be the eigenvalue decomposition of $$M$$ and furthermore take the eigenvalues to be sorted in decreasing order: $$\lambda_1 \ge \ldots \ge \lambda_n$$. Then the value of $$x^\top Mx$$ is likely going to be dominated by the projection of $$x$$ onto $$v_1$$:

\begin{align} x^\top M^kx & = \sum_i\lambda_i^k|v_i^\top x|^2 \\ & = \lambda_1^k\left(|v_1^\top x|^2 + \sum_{i = 2}\left(\frac{\lambda_i}{\lambda_1}\right)^k|v_i^\top x|^2\right) \end{align}

The important part to note is that the sum in parentheses is negligible asymptotically in the limit of large $$k$$ assuming that the largest eigenvalue has a multiplicity of 1. You can bound the relative error like so:

\begin{align} \frac{x^\top M x}{\lambda_1^k|v_1^\top x|^2} - 1 & = \sum_{i = 2}\left(\frac{\lambda_i}{\lambda_1}\right)^2|v_i^\top x|^2 \\ & \le \left(\frac{\lambda_2}{\lambda_1}\right)^k\sum_{i = 2}|v_i^\top x|^2 \\ & \le \left(\frac{\lambda_2}{\lambda_1}\right)^k|x|^2 \end{align}

where going from the first to the second line I used the fact that the eigenvalues are in decreasing order. Calculating the largest two eigenvalues and corresponding eigenvectors of $$M$$ is a whole lot easier than calculating the entire eigenvalue decomposition. All you do is use the power method on two vectors that you keep orthogonal at every iteration.

There are a few ways that this analysis could be wrong. The projection of $$x$$ onto $$v_1$$ could be very close to 0. You could also have a matrix where $$\lambda_1$$ and $$\lambda_2$$ are very close or nearly equal. In those instances, you might need to go to higher terms in the asymptotic expansion. But you still have a way of bounding the error after the fact while doing a lot less work than calculating a full eigenvalue decomposition.

• This idea should also carry over tot he matrix exponential. In practice, to use such an approach, you could use a black box algorithm to compute the top $j$ eigenvalues/vectors of $M$. In many cases those algorithms are Lanczos based, and really the Lanczos approximation in my response is implicitly doing something like this. But this approach sort of sheds a bit more light on what's actually happening. – tch Feb 26 at 15:55
• The dirty secret of Lanczos is that you need reorthogonalization and not every package does that well. I learned this the hard way after getting bad results from scipy and then coding this myself. (scipy uses ARPACK under the hood.) SLEPc on the other hand appears to be totally unbreakable. Lanczos is a great trick, but I really don't recommend people use black box eigenvalue solvers without doing a bit of a posteriori sanity checking. – Daniel Shapero Feb 26 at 17:59
• @DanielShapero thanks so much for this intuitive explanation. I am worried that in fact I might have a set of eigenvalues that are quite close to each other. On a different note can you please let me know why I cannot use SLEPc for this problem? – Cupitor Feb 27 at 18:23
• Sorry if I was unclear. SLEPc is a fantastic package and you totally could use it for this problem. But it is a little harder to use than scipy because the Python interface is more or less a straight translation of the underlying C code. – Daniel Shapero Mar 1 at 2:19