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:**55979**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:**55979**Joined:****Location:**Amsterdam-
**Contact:**

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

There are a couple of loose ends I'd like to tie up. Nothing important you understand.

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:**55979**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:**222**Joined:**

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

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.

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

- FaridMoussaoui
**Posts:**222**Joined:**

ISayMoo wrote: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.

Your COMPLEX function should satisfy the Cauchy-Riemann equations to be differentiable.

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

FaridMoussaoui wrote:ISayMoo wrote: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.

Your COMPLEX function should satisfy the Cauchy-Riemann equations to be differentiable.

AFAIR this is the same as saying that the function is holomorphic. A restrictive assumption. How would you compute the Hessian?

- FaridMoussaoui
**Posts:**222**Joined:**

The Hessian of a complex function?

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

FaridMoussaoui wrote:The Hessian of a complex function?

Why not?

- FaridMoussaoui
**Posts:**222**Joined:**

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

ISayMoo wrote:

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}[$]