3
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – Federico Poloni Apr 11 at 7:32
  • $\begingroup$ 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 + ...$ $\endgroup$ – Maxim Umansky Apr 11 at 14:55
  • 1
    $\begingroup$ @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. $\endgroup$ – Anton Menshov Apr 11 at 19:02
  • $\begingroup$ @MaximUmansky nope, $e^A$ is not close to unity. $\endgroup$ – Anton Menshov Apr 11 at 19:02
  • $\begingroup$ (What is $\|A\|$, just to have an idea?) $\endgroup$ – Federico Poloni Apr 11 at 19:05
7
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Anton Menshov Apr 11 at 19:05
2
$\begingroup$

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

enter image description here

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Alone Programmer Apr 12 at 3:07
  • $\begingroup$ @AloneProgrammer, I tested it with the example provided by Federico Poloni. $\endgroup$ – nicoguaro Apr 12 at 3:31
  • 1
    $\begingroup$ 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. $\endgroup$ – Maxim Umansky Apr 12 at 6:14
  • $\begingroup$ 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. $\endgroup$ – Maxim Umansky Apr 12 at 7:05
  • $\begingroup$ Thanks for the explanation, now I understand what you mean! $\endgroup$ – Federico Poloni Apr 12 at 7:50
0
$\begingroup$

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

| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ I think that's just how matrix logarithms work $\endgroup$ – Thijs Steel Apr 11 at 5:53
  • 2
    $\begingroup$ 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. $\endgroup$ – Amit Hochman Apr 11 at 8:21
  • 1
    $\begingroup$ 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. $\endgroup$ – Federico Poloni Apr 11 at 20:33
  • $\begingroup$ @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. $\endgroup$ – Maxim Umansky Apr 12 at 5:51
  • $\begingroup$ 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). $\endgroup$ – Federico Poloni Apr 12 at 7:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.