# Most accurate way to calculate matrix powers and matrix exponential for a positive semidefinite matrix

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}$$ is a PSD matrix, where $$k$$ can get quite 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$$.
• 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 and then exponentiate the singular values of $$M$$
• 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. Not knowing if here or stackoverflow being a good place to share this question at, I will be cross-posting it on their with the same title and content.

• Scicomp stack exchange is a more appropriate place. For what concerns 2, computing the eigendecomposition is a bad idea generally. If the size of the matrix is large, this is going to be even much more problematic. Also, expm computes the full exponential, which is not what you need, because you need the action of such matrix on a vector $x$. So, other techniques should be used (for a survey, take a look at Cleve Moler's famous paper "Nineteen dubious way ..." ).
– VoB
Feb 25, 2021 at 19:44
• There are some robust routines to accomplish this goal now, some of them don't use the positive semidefinitness of the matrix, but they work. There's a method based on rational approximations of $e^{x}$ for $x<=0$ that in the end requires to solve about 10 shifted linear systems. Those techiniques, for sure, are not present in Python
– VoB
Feb 25, 2021 at 19:48
• Also, forget naive Taylor expansions. In the mentioned paper you can see what happens if you try to use them
– VoB
Feb 25, 2021 at 19:49
• Don't to a Taylor expansion! The monomial basis $1,x,x^2, \ldots$ is terribly conditioned. If you want to compute it by a polynomial series, you'd be better of with Chebyshev polynomials. Alternately, you might look into using Lanczos. There is a fair amount of analysis on computing the matrix exponential via Lanczos. It should be reasonably straightforward and inexpensive to compute to machine precision this way. Feb 25, 2021 at 20:23
• There are even some python packages to compute $\exp(M)x$: weinbe58.github.io/QuSpin/generated/…. You can take that output and then take the inner product with $x$ to get the quadratic form you want. Feb 25, 2021 at 20:24

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" $$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$$. There are software packages to compute $$Q_k f(T_k) Q_k^T x$$. This can easily be turned into the Lanczos approximation to $$x^Tf(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, but without reorthogonalization $$k$$ may need to be larger. However, reorthogonalizaton is more computationally expensive, so there is some tradeoff.

• thanks so much! The question is also posted in scicomp and you might want to add your answer there too: scicomp.stackexchange.com/questions/36922. I think your answer covers the $e^M$ case. Do you have any idea of how to handle $M^k$? I assume by taking $f$ to be quadratic? Also the wiki link does not clarify how to replace $M$ with $f(M)$ could you link to something for that. Feb 25, 2021 at 21:35
• btw, the Golub Meurant et. al. reference is so useful for theoretical understanding. I wish I could have double-upvoted the post. Thank you so much! Feb 25, 2021 at 21:47
• You can do it the same way with $f(x) = x^i$. If $i<k$, then Lanczos will be exact. However, again computing $A^i$ is not a good idea numerically. Feb 25, 2021 at 23:19
• again I don't get how I can aply $f$? I don't see it in the wiki article either. How does one apply the $f$, considering that $T_k$ is a tridiagonal matrix? Feb 26, 2021 at 0:16
• Yeah, but this can be done in $O(k^2)$ becaues $T_k$ is tridiagonal. Most numerical libraries have a function for the eigendecomposition of a tridiagonal matrix. Feb 26, 2021 at 4:13