Using matrix exponential to solve linear system

Question

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?

Answer 1

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.

Answer 2

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()

Answer 3

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