Simple question which stumps me: I have a complex square matrix H. There are some nice methods for calculating exp(H). What about calculating the derivative of exp(H) over elements of H? To be precise: let M = exp(H). I want to calculate dM_{jk} / dH_{mn} numerically, accurately and (relatively) quickly.

- Cuchulainn
**Posts:**58401**Joined:****Location:**Amsterdam-
**Contact:**

I assume you want the Frechet derivative.

Offhand

1. standard finite differences

2. Complex step method

3. Other (Pade)

1 is ill-posed and 2 was discussed here in the scalar case.

http://eprints.ma.man.ac.uk/1260/1/cove ... 009_31.pdf

Offhand

1. standard finite differences

2. Complex step method

3. Other (Pade)

1 is ill-posed and 2 was discussed here in the scalar case.

http://eprints.ma.man.ac.uk/1260/1/cove ... 009_31.pdf

Thanks, this looks like what I was looking for.

One more question, as Lt Colombo says: can you recommend a good algorithm to compute the exponential of a complex matrix?

- Cuchulainn
**Posts:**58401**Joined:****Location:**Amsterdam-
**Contact:**

There are a couple of loose ends I'd like to tie up. Nothing important you understand.One more question, as Lt Colombo says: can you recommend a good algorithm to compute the exponential of a complex matrix?

Maybe

1. Compute the real and imaginary parts of exp(A + iB)

2. Pade approximation

3. Solve dV/dt = HV as an ODE solver I've done it in Boost odeint in the real case (odeint supports complex-valued matrices).

Matlab/Mathematica might document how they do it?

For 2, you could test on a scalar problem and a 2X2 matrix to prove/disprove. Use Pade as a kind of surrogate.

[$]r(z) = (2-z)/(2+z)[$]

- Cuchulainn
**Posts:**58401**Joined:****Location:**Amsterdam-
**Contact:**

If you have a continuous-time dynamic system it might be easier to solve a matrix ODE in matrix [$]M(t):[$]

[$]dM/dt = MA[$] with [$]M(0) = C[$].

Solution [$]M(t) = C e^{At}[$]

It's worked out in forthcoming C++ 2nd ed. book

[$]dM/dt = MA[$] with [$]M(0) = C[$].

Solution [$]M(t) = C e^{At}[$]

It's worked out in forthcoming C++ 2nd ed. book

Code: Select all

```
class MatrixOde
{
private:
// dB/dt = A*B, B(0) = C;
boost::numeric::ublas::matrix<value_type> A_;
boost::numeric::ublas::matrix<value_type> C_;
public:
MatrixOdePde(const boost::numeric::ublas::matrix<value_type>& A, const boost::numeric::ublas::matrix<value_type>& IC)
: A_(A), C_(IC) {}
void operator()( const state_type &x , state_type &dxdt , double t ) const
{
for( std::size_t i=0 ; i < x.size1();++i )
{
for( std::size_t j=0 ; j < x.size2(); ++j )
{
dxdt(i, j) = 0.0;
for (std::size_t k = 0; k < x.size2(); ++k)
{
dxdt(i, j) += A_(i.k)*x(k.j);
}
}
}
}
};
```

- FaridMoussaoui
**Posts:**356**Joined:**

Al-Mohy, A. H. and N. J. Higham, “A new scaling and squaring algorithm for the matrix exponential,” SIAM J. Matrix Anal. Appl., 31(3) (2009), pp. 970–989.One more question, as Lt Colombo says: can you recommend a good algorithm to compute the exponential of a complex matrix?

use your favourite search engine to find this article. You can also use sci-hub to download ANY published paper (for free).

- FaridMoussaoui
**Posts:**356**Joined:**

Your COMPLEX function should satisfy the Cauchy-Riemann equations to be differentiable.Simple question which stumps me: I have a complex square matrix H. There are some nice methods for calculating exp(H). What about calculating the derivative of exp(H) over elements of H? To be precise: let M = exp(H). I want to calculate dM_{jk} / dH_{mn} numerically, accurately and (relatively) quickly.

- Cuchulainn
**Posts:**58401**Joined:****Location:**Amsterdam-
**Contact:**

AFAIR this is the same as saying that the function is holomorphic. A restrictive assumption. How would you compute the Hessian?Your COMPLEX function should satisfy the Cauchy-Riemann equations to be differentiable.Simple question which stumps me: I have a complex square matrix H. There are some nice methods for calculating exp(H). What about calculating the derivative of exp(H) over elements of H? To be precise: let M = exp(H). I want to calculate dM_{jk} / dH_{mn} numerically, accurately and (relatively) quickly.

- FaridMoussaoui
**Posts:**356**Joined:**

The Hessian of a complex function?

- Cuchulainn
**Posts:**58401**Joined:****Location:**Amsterdam-
**Contact:**

Why not?The Hessian of a complex function?

- FaridMoussaoui
**Posts:**356**Joined:**

something like [$] \frac{d^2 f}{dz_i d\overline{z}_j} [$] for a function [$] f(z_1, z_2, ..., z_n) [$]

if it's diagonalizable,

[$]A =PDP^{-1}[$]

[$]\exp(A)=P \exp(D)P^{-1}[$]

The function is exp, so it's everything you need.

if [$]D[$] is diagonal, with diagonal entries [$]\lambda_{1},\lambda_{2},\lambda_{3},\cdots[$],

[$]\exp(D)[$] is also diagonal, with diagonal entries [$]\exp(\lambda_{1}),\exp(\lambda_{2}),\exp(\lambda_{3}),\cdots[$]

[$]A=PDP^{-1}[$] where [$]P[$] is the matrix whose columns are the eigenvectors of [$]A[$] and [$]D[$] is a diagonal matrix whose elements are the eigenvalues of [$]A[$]

notice that

[$]A^{2}=PDP^{-1}PDP^{-1}=PD^{2}P^{-1}[$] and similarly [$]A^{3}=PD^{3}P^{-1}[$], and [$]A^{n}=PD^{n}P^{-1}[$]

so [$]\exp(A)=\sum_{n=0}^{\infty}\frac{1}{n!}A^{n}=P\left[\sum_{n=0}^{\infty}\frac{1}{n!}D^{n}\right]P^{-1}=P\exp(D)P^{-1}[$]

[$]\exp(D)[$] is also diagonal, with diagonal entries [$]\exp(\lambda_{1}),\exp(\lambda_{2}),\exp(\lambda_{3}),\cdots[$]

[$]A=PDP^{-1}[$] where [$]P[$] is the matrix whose columns are the eigenvectors of [$]A[$] and [$]D[$] is a diagonal matrix whose elements are the eigenvalues of [$]A[$]

notice that

[$]A^{2}=PDP^{-1}PDP^{-1}=PD^{2}P^{-1}[$] and similarly [$]A^{3}=PD^{3}P^{-1}[$], and [$]A^{n}=PD^{n}P^{-1}[$]

so [$]\exp(A)=\sum_{n=0}^{\infty}\frac{1}{n!}A^{n}=P\left[\sum_{n=0}^{\infty}\frac{1}{n!}D^{n}\right]P^{-1}=P\exp(D)P^{-1}[$]

GZIP: On