I want to calculate the nth power of a matrix A. However, when I try A^n, it gives me the element-by-element power and not the matrix power. Is there any easy way to calculate A^n rather than A*A*...*A?

Thanks a lot for any help!

Szabolcs

## 2 Answers

0

**Update:**

GAUSS now has a function to compute matrix powers more efficiently. You can find the documentation here: powerM.

**Original post:**

The easiest way to calculate the power of a matrix is probably using the eigenvalues (which we will call `lambda`) and the eigenvectors (which we will call `e_vec`). The steps are:

- Calculate the eigenvalues and eigenvectors:
`{ lambda, e_vec } = eigv(A);`

- Diagonalize the eigenvalues:
`D = lambda .* eye(rows(A));`

- Multiply:
`power = 3; A_to_3 = e_vec * D.^power * inv(e_vec);`

Here is a GAUSS procedure that with a simple example:

```
A = rndn(4, 4);
power = 3;
A_3 = A*A*A;
A_3_powmat = powMat(A, power);
print "A_3 = " A_3;
print "A_3_powmat = " A_3_powmat;
proc (1) = powMat(A, power);
local e_val, e_vec, A_pow;
//Calculate eigenvalues
{ e_val, e_vec } = eigv(A);
//Calculate power of matrix
A_pow = e_vec * (e_val .* eye(rows(A))).^power * inv(e_vec);
//The eigenvalues for many matrices are complex
//and will return a complex portion that is round-off error
//remove the round-off error for real matrices
if iscplx(A);
retp(A_pow);
else;
retp(real(A_pow));
endif;
endp;
```

0

Thanks a lot for the answer, that is one possible way to calculate it. After digging into the language reference, I have found a command that seems even easier.

For example, if I want to calculate the *n*th power of a matrix A, I can do the following:

`A_pow = polyeval(A, 1|zeros(n,1))`

I do not know how this function is evaluated, so I can't say anything about how round-off errors are accumulated with the two methods. If the matrix contains very small numbers, the eigenvalue decomposition and `polyeval`

can lead to different results, but no idea which is more reliable.

Thanks again for your help.

Szabolcs

## Your Answer

## 2 Answers

**Update:**

GAUSS now has a function to compute matrix powers more efficiently. You can find the documentation here: powerM.

**Original post:**

The easiest way to calculate the power of a matrix is probably using the eigenvalues (which we will call `lambda`) and the eigenvectors (which we will call `e_vec`). The steps are:

- Calculate the eigenvalues and eigenvectors:
`{ lambda, e_vec } = eigv(A);`

- Diagonalize the eigenvalues:
`D = lambda .* eye(rows(A));`

- Multiply:
`power = 3; A_to_3 = e_vec * D.^power * inv(e_vec);`

Here is a GAUSS procedure that with a simple example:

```
A = rndn(4, 4);
power = 3;
A_3 = A*A*A;
A_3_powmat = powMat(A, power);
print "A_3 = " A_3;
print "A_3_powmat = " A_3_powmat;
proc (1) = powMat(A, power);
local e_val, e_vec, A_pow;
//Calculate eigenvalues
{ e_val, e_vec } = eigv(A);
//Calculate power of matrix
A_pow = e_vec * (e_val .* eye(rows(A))).^power * inv(e_vec);
//The eigenvalues for many matrices are complex
//and will return a complex portion that is round-off error
//remove the round-off error for real matrices
if iscplx(A);
retp(A_pow);
else;
retp(real(A_pow));
endif;
endp;
```

Thanks a lot for the answer, that is one possible way to calculate it. After digging into the language reference, I have found a command that seems even easier.

For example, if I want to calculate the *n*th power of a matrix A, I can do the following:

`A_pow = polyeval(A, 1|zeros(n,1))`

I do not know how this function is evaluated, so I can't say anything about how round-off errors are accumulated with the two methods. If the matrix contains very small numbers, the eigenvalue decomposition and `polyeval`

can lead to different results, but no idea which is more reliable.

Thanks again for your help.

Szabolcs