# Using matrix exponential to solve linear system

Consider the system of linear equations:

$$Ax=b \tag{1} \label{eq1}$$

where

• $$A\in\mathbb F^{n\times n}$$, diagonalizable dense matrix, over the field $$\mathbb F$$ of real or complex numbers,
• $$x\in \mathbb F^{n\times 1}$$ is a vector of unknowns
• $$b\in \mathbb F^{n\times 1}$$ is a known right-hand side vector
• $$n$$ is in the order of 1000–10000

Unlike the usual system of linear equations, I do not have access to $$A$$ itself; however, I have access to the matrix exponential $$e^A$$. The matrix exponential is accessible both as an explicit matrix and, consequently, as a function acting on a vector.

What are my options for finding the solution of $$\eqref{eq1}$$ knowing only the matrix exponential? Computing matrix logarithm does not seem the best option as I doubt being able to get anything numerically stable and reasonably efficient.

Am I missing something simple?

• 1. Have you tried some ready-to-use methods for the matrix logarithm such as Matlab's logm(E) \ b, before writing it off as unstable / inefficient? 2. Are $A$ and $E=\exp(A)$ sparse or dense? – Federico Poloni Apr 11 at 7:32
• Is $B=e^A$ close to unity? Then you can use Taylor expansion $A = log(B) = (B-I) + (B-I)^2/2 - (B-I)^3/3 + ...$ – Maxim Umansky Apr 11 at 14:55
• @FedericoPoloni $A$ and $E$ are dense. I would have to export the matrix to Matlab and try it there, which I should do. That is a "bruteforce" method in my initial assessment – and I certainly have to try it at least on a couple of examples. – Anton Menshov Apr 11 at 19:02
• @MaximUmansky nope, $e^A$ is not close to unity. – Anton Menshov Apr 11 at 19:02
• (What is $\|A\|$, just to have an idea?) – Federico Poloni Apr 11 at 19:05

You are effectively asking how to compute $$y=(\log M )^{-1}b$$, where $$M=e^A$$ is the given matrix. There are several methods for computing $$f(M)b$$ without forming $$f(M)$$, and they are reviewed here. One general method is to use Cauchy’s theorem, $$y=\dfrac{1}{2\pi i}\int_\Gamma f(z)(zI - M)^{-1}b\,dz,$$ with $$f(x) = 1/\log(x)$$. $$\Gamma$$ is a contour that encloses all the eigenvalues of $$M$$, so you need to first estimate the magnitude of the largest eigenvalue, say, with a power method. Then you approximate the integral with the trapezoidal rule. You need to solve several shifted systems of the form $$(zI-M)x =b$$, for which a preliminary reduction to Hessenberg form is useful.

• great reference. It seems like exactly what I hoped for. I will have to analyze it against computing $\log M$, but now I have an alternative approach. – Anton Menshov Apr 11 at 19:05

I am expanding my comment into an answer. I don't think that the method is efficient, but I think that it can be used to obtain the matrix $$A$$ from $$e^{A}$$.

We know that

$$\frac{d e^{tA}}{dt} = e^{tA} A\, ,$$

so, we could use

$$\left.\frac{d e^{tA}}{dt}\right|_{t=0} = A\, ,$$

if we can approximate the derivative

$$\frac{d e^{tA}}{dt} \approx D(A)\, .$$

For example, we could use a forward finite difference

$$\left.\frac{d e^{tA}}{dt}\right|_{t=0} \approx \frac{e^{hA} - I}{h}\, ,$$

but the problem is that we need to compute the fractional power of the matrix $$e^{A}$$. Maybe we could take advantage of higher-order approximations and just use integer powers of the matrix, but a couple I tried didn't work properly.

It seems to work, but I doubt it is efficient.

import numpy as np
from scipy.linalg import logm, fractional_matrix_power as powm
import matplotlib.pyplot as plt

eA = np.array([
[1, -1],
[1, 2]])
A = logm(eA)
rel_error = []
steps = [1, 1e-1, 1e-2, 1e-3, 1e-4, 1e-5]
for h in steps:
A1 = np.real((powm(eA, h) - np.eye(2))/h)
rel_error.append(np.linalg.norm(A - A1)/np.linalg.norm(A))

plt.loglog(steps, rel_error)
plt.xlabel("Relative error")
plt.xlabel("$$h$$")
plt.savefig("matexp.png", dpi=300, bbox_inches="tight")
plt.show()


• I think you used a pretty bad example to test your method here. Based on OP's question, his matrix is in the order of $1000 \times 1000$ to $10000 \times 10000$, and I believe your method won't converge at all for a $1000 \times 1000$ matrix. Replace your $2 \times 2$ matrix by this n = 1000 eA = np.random.rand(n,n) and you would see by decreasing $h$ you won't get a convergence. – Alone Programmer Apr 12 at 3:07
• @AloneProgrammer, I tested it with the example provided by Federico Poloni. – nicoguaro Apr 12 at 3:31
• I guess your finite-difference expression can be also interpreted as a consequence of the property (or definition) of the exponential function, $exp(x) = lim(1+x/n)^n$ for large n, or $x = lim (n(e^{x/n}-1))$. If we had a good way to find the $n^{th}$ root of our given matrix $e^A$ then this would be indeed a possible way to approximate matrix A. – Maxim Umansky Apr 12 at 6:14
• Here is a good survey of calculation methods for finding a nth root of a matrix, cis.upenn.edu/~cis610/bini_higham_pth_root_na_2005.pdf. – Maxim Umansky Apr 12 at 7:05
• Thanks for the explanation, now I understand what you mean! – Federico Poloni Apr 12 at 7:50

If we diagonalize matrix $$A$$ by finding the transformation S such that $$A = S D S^{-1}$$ where D is a diagonal matrix and the diagonal elements of $$D$$ are the eigenvalues $$\lambda_k$$ then the same transformation makes $$e^A$$ diagonal, and the eigenvalues are $$e^{\lambda_k}$$. So diagonalizing $$e^A$$ and taking log of eigenvalues we find matrix $$D$$, which is sufficient to solve the linear system; and using the transformation $$S$$ we can find the original matrix $$A.$$

Alternatively, if $$B=e^A$$ is close enough to unity, you can use Taylor expansion to find $$A$$, $$A = log(B) = (B-I) + (B-I)^2/2 - (B-I)^3/3 + ...$$

• I think that's just how matrix logarithms work – Thijs Steel Apr 11 at 5:53
• That’s how matrix logarithm’s work in exact arithmetic and for normal matrices in floating point. For general diagonalizable matrices, $S$ mat be ill-conditioned so that significant errors may result from this approach. – Amit Hochman Apr 11 at 8:21
• You are slowly re-inventing the standard methods to compute logarithms of dense matrices. :) The next step is using the identity $\log(A) = 2^s \log(A^{1/2^s})$ to get a matrix closer to the identity, and replace Taylor with Padé which has a better convergence radius. – Federico Poloni Apr 11 at 20:33
• @Federico Poloni How about other expansions for $ln(x)$, in particular the second one in math.com/tables/expansion/log.htm? I guess it can be applied for matrices too, if we first calculate the inverse of our given $e^A$; and convergence of this series is not restricted to norm $|| e^A-I ||$ < 1. – Maxim Umansky Apr 12 at 5:51
• It's possible that it could be made to work, but I have no direct experience. Convergence of that series is restricted to matrices with no eigenvalues $\lambda < 1/2$, if I understand correctly (not clear how this generalizes to complex numbers). The problem is, it is easy to ensure that a matrix has no large eigenvalues (small norm implies that), but it is more difficult to ensure that it has no small eigenvalues). – Federico Poloni Apr 12 at 7:48