How can I multiply an NxN matrix A in Fortran x times to get its power without amplifying rounding errors?
Asked
Active
Viewed 330 times
0
-
How precise does your computation needs to be? Have you tried using 64, 128 bit floating point numbers? It would be great if you could show a test case with results and expectations – solalito Nov 16 '15 at 11:02
2 Answers
3
If A can be diagonalized as
A P = P D,
where P is some NxN matrix (each column is called 'eigenvector'), and D is an NxN diagonal matrix (the diagonal elements are called 'eigenvalues'), then
A = P D P^{-1},
where P^{-1} is the inverse matrix of P. Therefore the second power of A results in
A A= P D P^{-1} P D P^{-1} = P D D P^{-1}.
Repeating multiplication of A for x times yields
A^x = P D^x P^{-1}.
Note here that D^x is still a diagonal matrix. Let the i-th diagonal element of D be D_{ii}. Then, the i-th diagonal element of D^x is
[D^x]_{ii} = (D_{ii})^x.
That is, the elements of D^x is simply x-th power of the elements of D and can be computed without much rounding error, I guess. Now, you multiply P and P^{-1} from left and right, respectively, to this D^x to obtain A^x. The error in A^x depends on the error of P and P^{-1}, which can be calculated by some subroutines in numerical packages such as LAPACK.
norio
- 3,652
- 3
- 25
- 33
2
as alluded to in the answer by norio, one can employ in general the Jordan (or alternatively Schur) decomposition and proceed in a similar fashion - for details (including brief error analysis) see, e.g., Chapter 11 of Matrix computations by Golub and Loan.
ewcz
- 12,819
- 1
- 25
- 47