# How to properly do matrix exponents?

sleepingawake86

Wed, 11/02/2022 - 12:12 am

Hello,

I have a matrix A that I want to take to the y power. Doing A^y or matrixop Out=powC(A,y) are not giving me the right answer.

However, if I manually do the matrix multiplication with:

MatrixOp Out=A x A x A x A.....

I get the right answer. However this method is infeasible for large y (I am trying to do up to y=50000). What command do I need to use to get me the correct answer without manually typing this.

You probably have picked the wrong command. Use powR instead of powC; the latter is for complex values. Does it work then? Otherwise, can you explain in a bit more detail what is 'wrong'?

November 2, 2022 at 06:14 am - Permalink

`displayhelptopic "MatrixEigenV"`

Have you tried using MatrixEigenV? Once you have the diagonal eigenvalues, raise each to the desired power, then re-transform with the left and right eigenvector matrices.

November 2, 2022 at 06:19 am - Permalink

50000? Is that even feasible in double precision? I'm curious...

November 2, 2022 at 10:05 am - Permalink

In reply to You probably have picked the… by chozo

A is complex.

What is wrong is that if I do MatrixOp A x A and then powC(A,2) I get different answers.

powC(A,2) gives me the same result as typing A^2

A x A gives me the same result as MatrixMultiply A x A

For some reason powC and ^ are not doing the expected matrix multiplication.

November 2, 2022 at 11:53 am - Permalink

I see your point now. powC / powR does indeed seem to do only scalar multiplication. I have no idea if there is a command for proper repeated matrix multiplication. Maybe someone else knows. By the way, I find a factor of 50000 also a bit extreme. What is your use case here?

November 2, 2022 at 07:04 pm - Permalink

Here is a snippet of a working example using MatrixEigenV. A caveat is that I have tested it only with a small Hermitian matrix. The example uses a power of 8 for demonstration. Note the use of powC( ) on the eigenvalues. (IP 8.04, Win10). In principle, after diagonalization, the entire power operation could be written in one MatrixOP line.

Make/D/N=(2,2)/O/C H={{cmplx(3,0),cmplx(2,-1)} , {cmplx(2,1),cmplx(4,0)}} // Hermitian

MatrixEigenV/R/C H // here is the core of the algorithm

Wave W_eigenvalues, M_R_eigenVectors // created by MatrixEigenV

MatrixOP/O/C MdiagH = diagonal(W_eigenvalues) // square matrix

MatrixOP/O/C prodH = M_R_eigenVectors x (MdiagH x inv(M_R_eigenVectors))) // checks with H

MatrixOP/O/C HprodH = H x H x H x H x H x H x H x H

MatrixOP/O/C HprodHdiag = powC(W_eigenvalues,8) // EIGENVALUES raised to power

MatrixOP/O/C MdiagHprod = diagonal(HprodHdiag)

MatrixOP/O/C prodHpwr = M_R_eigenVectors x (MdiagHprod x inv(M_R_eigenVectors))) //checks with HprodH

end

November 3, 2022 at 08:59 am - Permalink

In reply to A is complex. What is wrong… by sleepingawake86

s.r.chinn posted the proper way to calculate the power.

Powc() PowR() work on (matrix) element-by-element basis. This is NOT matrix multiplication. Try, for example:

Make/O/N=(2,2) yyy={{1,2},{3,4}}

Matrixop/o/p=1 aa=powr(ddd,yyy)

aa[0]={2,4,8,16}

November 3, 2022 at 09:32 am - Permalink

It seems to me that one plausible scenario for expanding to a power of 50000 is that the matrix is unitary. Then the magnitude of the determinant must equal 1, and each eigenvalue is of the form exp( i lambda), also magnitude 1. This type of matrix can be thought of as a generalized rotation. If the matrix is not unitary, I don't see how numeric overflow or underflow can be avoided for arbitrarily large powers.

November 5, 2022 at 06:30 am - Permalink

This is actually a hard task, numerically. There is a classic article on this, "Nineteen Dubious Ways to Compute the Exponential of a Matrix, Twenty-Five Years Later" by Cleve Moler and Charles Van Loan. (Moler was a creator of Matlab.)

See https://epubs.siam.org/doi/abs/10.1137/S00361445024180

November 9, 2022 at 05:41 pm - Permalink