How to properly do matrix exponents?

Share icon

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

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.

function bar()

 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

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

sleepingawake86 wrote:

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

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

Make/O/N=(2,2) ddd=2
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