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

Question

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 = [email protected](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)
    results.append((xexpMx - lanczos_quad_form(np.exp,a_,b_,np.linalg.norm(x)))/xexpMx)

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!

Answer 1

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)

def lanczos_quad_form(f,a_,b_,norm_b=1):
    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)
lanczos_quad_form(np.exp,a_,b_,np.linalg.norm(x))

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 = [email protected](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)
(xexpMx - lanczos_quad_form(np.exp,a_,b_,np.linalg.norm(x)))/xexpMx

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

Answer 2

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

Answer 3

tch already gave you a good answer, but I have a suggestion for a simpler problem that you could address pretty easily and which should help you sanity check your results. Let

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