NumericalMethods[PowerIterate] - approximates the largest eigenvalue and its eigenvector of a matrix

Calling Sequence
    PowerIterate(A, xinit=v, tolerance=t, options)

Parameters
    A - Matrix; a square numerical matrix
    v - Vector; initial approximate eigenvector of appropriate dimension
    t - positive; tolerance for the stopping condition
    options - (optional) parameters of the form keyword = value, where keyword is one of method, matrix, stoppingcriterion, output, vector, maxiterations, maxnumber, digits.

Description
The PowerIterate command approximates the largest eigenvalue (real or complex) and its associated eigenvector, assuming that the largest eigenvalue is unique. It also allows for using a shifted matrix and for inverse power iteration. Its main purpose is to exhibit behaviour of basic methods and be accessible to students not versed in Maple. Acordingly, the main focus is on ease of use and informative output, not performace. Controls are similar to other iterative procedures in the package.

Two sequences are constructed iteratively, vector approximations  x_k and eigenvalue approximations l_k.

There is no default tolerance, it must be user specified.

It is assumed that the package LinearAlgebra has ben loaded before this procedure is read.

This command is part of the NumericalMethods package, so it can be used in the form PowerIterate(..) only after executing the command with(NumericalMethods). However, it can always be accessed through the long form of the command by using NumericalMethods[PowerIterate](..).

Options

The options argument can be one or more of the following equations.

method = Rayleigh, normalized, or RayleighQuotient
- Rayleigh: constructs new iterations in three steps: The first step sets up an auxiliary vector y = Ax_k / l_k. Then it is normalized with respect to infinity (max) norm: x_k+1 = y / |y|. Finally, a new approximation of eigenvalue is determined using the Rayleigh quotient: l_k+1 = (x_k+1^T *Ax_k+1) / (x_k+1^T * x_k+1), for complex matrices and/or vectors Hermitian transpose is used.
- normalized: the traditional plain iteration: x_k+1 = Ax_k / |Ax_k |. Then a new approximation of eigenvalue is determined using the Rayleigh quotient: l_k+1 = (x_k+1^T *Ax_k+1) / (x_k+1^T * x_k+1), for complex matrices and/or vectors Hermitian transpose is used. Calculations are fast, but this iteration does not converge when the largest eigenvalue is negative or complex.
- RayleighQuotient: another three step iteration, it uses inverse matrix for faster convergence. The first step sets up an auxiliary vector y = (A - l_k E_n )^-1 x_k / l_k. Then it is normalized with respect to infinity norm: x_k+1 = y / |y|. Finally, a new approximation of eigenvalue is determined using the Rayleigh quotient: l_k+1 = (x_k+1^T *Ax_k+1) / (x_k+1^T * x_k+1), for complex matrices and/or vectors Hermitian transpose is used.
The default value is method = Rayleigh.

matrix = [shifted,numeric] or [inverse,numeric]

The iteration is applied to a transformed matrix, the resulting eigenvalue is then transformed back to yield approximation for an eigenvalue of A. In both cases eigenvectors are common to the original and transformed matrix.
- [shifted,mu], where mu is numeric: power iteration is applied to (A - mu*E_n ), it yields the largest eigenvalue of this shifted matrix. Can be used for instance to break a tie when l and -l are largest eigenvalues. This happens in particular when a real matrix has complex eigenvalues, then one has to shift by a complex number.
- [inverse,mu], where mu is numeric: power iteration is applied to (A - mu*E_n ) ^-1, resulting in the inverse power iteration. It yields the eigenvalue of A that is closest to mu.
By default, no transformation is applied.

stoppingcriterion = vectors or lambdas
- vectors: the absolute difference between last two eigenvector approximations in maximum norm is compared to tolerance.
- lambdas: the relative difference between the last two eigenvalue approximations is compared to tolerance.
The default value is stoppingcriterion = vectors.

output = information, lambdas, iterates, value, or sequence
Specifies what information is given by the procedure. In all cases, the iteration index k is printed at every iteration step, so that the user can see the progress and knows how many iterations it took to produce the final approximation.
- information returns the final numerical approximation of the eigenvalue and also the final approximation of its eigenvector if specified by vector. During iterations, generated eigenvector and eigenvalue approximations are shown.
- lambdas returns the final numerical approximation of the eigenvalue and also the final approximation of its eigenvector if specified by vector. During iterations, generated eigenvalue approximations are shown.
- iterates returns the final numerical approximation of the eigenvalue and also the final approximation of its eigenvector if specified by vector. During iterations, generated eigenvector approximations are shown.
- value returns the final numerical approximation of the eigenvalue and also the final approximation of its eigenvector if specified by vector.
- sequence: returns all eigenvalue approximations produced by iterations as a sequence. Does not return the final approximation of eigenvector even if vector = show.
The default value is output = information if matrix dimension is less than 7 for real valued approximations and if matrix dimension is less than 4 for complex valued approximations. Otherwise the default is output = lambdas.

vector = show or noshow
Determines whether the final approximation of eigenvector should be included in the procedure's return.
The default value is vector = >show unless the output is set to output = sequence.

maxiterations = posint
The maximum number of iterations to perform. If it is reached, the algorithm is interrupted and a message is displayed. The return as defined by output is still produced.
The default value is maxiterations = 50.

maxnumber = positive
The maximal number that is not to be exceeded by approximations or function values during iterations. If it does happen, the algorithm is interrupted and a message is displayed.  The return as defined by output is still produced.
The default value is maxnumber = 1000000.

digits = posint
Specifies how many decimal places should be shown when output = information, output = lambdas or output = iterates is specified.
If digits is set to a value greater than Digits, the latter is set to digits for the duration of this procedure.
The default value is digits = Digits.

Note: When the chosen method is RayleighQuotient, it may happen that l_k is an eigenvalue at some step, so the appropriate inverse matrix cannot be calculated. If the vector option is noshow, then the algorithm terminates and shows lambda. However, if the vector option is show, then the algorithm cannot stop yet, because x_k need not be even close to the appropriate eigenvector. It therefore keeps iterating and the next step is done using normalized iteration.

Examples

> with(LinearAlgebra):
  A:=<<1,2,2>|<1,2,4>|<-1,1,-1>>;

We start with default settings.
> PowerIterate(A, xinit=<1.,13.,14.>, tolerance=0.01);
k=00  x=[  1.0000000000, 13.0000000000, 14.0000000000],    lambda=  3.0218579230 k=01  x=[  0.0000000000,  1.0000000000,  0.9523809524],    lambda=  3.0701545770
k=02  x=[  0.0156250000,  0.9687500003,  1.0000000000],    lambda=  2.9823699800
k=03  x=[ -0.0052631578,  0.9999999996,  0.9789473682],    lambda=  3.0206232230
k=04  x=[  0.0052447552,  0.9860139861,  1.0000000000],    lambda=  2.9891242520
k=05  x=[ -0.0029308323,  0.9999999996,  0.9906213361],    lambda=  3.0081112460
k=06  x=[  0.0021467603,  0.9937548792,  1.0000000000],    lambda=  2.9948495660

If the output turns out to be too wide, one can use the digits option, digits = 5 is often enough for educational purposes. In particular it does not make much sense to have digits significantly higher than tolerance.
Since A has a positive largest eigenvalue, the plain iteration with normalization should work as well. On the other hand, applying it to -A with negative largest eigenvalue means that the vector iterates will not converge. However, the Rayleigh estimates for eigenvalue still should. We therefore change the stopping condition. Note how vector iterates switch between two good candidates for eigenvectors.
> PowerIterate(-A, xinit=<1.,13.,14.>, tolerance=0.01, method=normalized, stoppingcriterion=lambdas, vector=noshow, digits=5);
k=00  x=[  1.00000, 13.00000, 14.00000],    lambda= -3.02186
k=01  x=[  0.00000, -1.00000, -0.95238],    lambda= -3.07015
k=02  x=[  0.01562,  0.96875,  1.00000],    lambda= -2.98237
k=03  x=[  0.00526, -1.00000, -0.97895],    lambda= -3.02062
k=04  x=[  0.00524,  0.98601,  1.00000],    lambda= -2.98912
k=05  x=[  0.00293, -1.00000, -0.99062],    lambda= -3.00811

The return is in the precision given by Digits.
We return to A and try another form of output and Rayleigh Quotient Iteration.
> PowerIterate(A, xinit=<1.,13.,14.>, tolerance=0.01, method=RayleighQuotient, output=value);
k=00 k=01 k=02 k=03 k=04

We still get to see how many iterations it took. The output = sequence option does not show the eigenvector even if we ask for it. We also try increased precision.
> PowerIterate(A,xinit=<1.,13.,14.>, tolerance=0.01, method=RayleighQuotient, vector=show, output=sequence, digits=15);
k=00 k=01 k=02 k=03 k=04

We can now work with this output.
> %[4];

Shifting the matrix "to the left" by three turns 3 into 0, so it will definitely not be the largest eigenvalue any more and thus we get to see another one.
> PowerIterate(A, xinit=<1.,13.,14.>, tolerance=0.01, matrix=[shift,3], digits=7);
Iterating with the matrix (A-3*E)
k=00  x=[  1.0000000, 13.0000000, 14.0000000],    lambda=  0.0218579
k=01  x=[ -1.0000000,  1.0000000, -0.6666667],    lambda= -4.2727273
k=02  x=[ -0.7857143,  0.7857143, -1.0000000],    lambda= -4.8538813
k=03  x=[ -0.6025641,  0.6025641, -1.0000000],    lambda= -4.9756237
k=04  x=[ -0.5394089,  0.5394089, -1.0000000],    lambda= -4.9960730
k=05  x=[ -0.5155189,  0.5155189, -1.0000000],    lambda= -4.9993710
k=06  x=[ -0.5061693,  0.5061693, -1.0000000],    lambda= -4.9998993

Since this is a school example, we expect "nice numbers", which allows us to make a good guess of what the precise eigenvalue and eigenvector is. Let's check:
> v:=<1,-1,2>: l:=-2:
  A.v-l*v;

Yes, we were right.
Now we know that the third eigenvalue is somewhere between -2 and 3. We will try the inverse matrix iteration and ask for the eigenvalue closest to mu=0.5.
> PowerIterate(A, xinit=<1.,13.,14.>, tolerance=0.01, matrix=[inverse,1/2], digits=7);
Iterating with the matrix (A-1/2*E)^(-1)
k=00  x=[  1.0000000, 13.0000000, 14.0000000],    lambda=  0.3125683
k=01  x=[  0.2727273,  1.0000000,  0.9090909],    lambda= -0.0034783
k=02  x=[ -1.0000000,  0.1764706,  0.1176471],    lambda=  2.3284768
k=03  x=[ -1.0000000,  0.8313253,  0.8433735],    lambda=  2.1887130
k=04  x=[ -1.0000000,  0.9664269,  0.9640288],    lambda=  2.0361903
k=05  x=[ -1.0000000,  0.9932789,  0.9937590],    lambda=  2.0071853
k=06  x=[ -1.0000000,  0.9986560,  0.9985600],    lambda=  2.0014342

As before we can confirm that the third eigenvalue is 1 with eigenvector [-1,1,1].
Now we try another matrix.
> B:=<<4,13>|<-1,0>>;

It has complex eigenvalues and being a real matrix, there are two of equal magnitude, meaning that iteration does not work.
Shifting by a real constant does not help, so we have to shift by a complex number to get one of the eigenvalues.
For complex eigenvalues it is recommended to start with a complex initial vector and we do it here just to show that we can, although in this case it is not needed as the complex shift is enough.
> PowerIterate(B, xinit=<1.3+I,1.4-I>, tolerance=0.001, matrix=[shift,2*I], digits=4);
Iterating with the matrix (A-2*I*E)
k=00  x=[  1.3000 +1.0000I,  1.4000 -1.0000I],    lambda=  3.6460 +4.6903I
k=01  x=[  0.3021 -0.1720I,  0.9524 -0.3048I],    lambda=  4.0728 -2.8964I
k=02  x=[  0.0943 -0.1613I,  0.9618 -0.2737I],    lambda=  1.6986 -3.7494I
k=03  x=[  0.0333 -0.2593I,  0.9669 -0.2553I],    lambda=  1.3614 -5.1745I
k=04  x=[  0.0898 -0.2872I,  0.9633 -0.2685I],    lambda=  2.1331 -5.2416I
k=05  x=[  0.0944 -0.2618I,  0.9623 -0.2720I],    lambda=  2.0922 -4.9397I
k=06  x=[  0.0848 -0.2600I,  0.9629 -0.2697I],    lambda=  1.9727 -4.9631I
k=07  x=[  0.0845 -0.2643I,  0.9631 -0.2693I],    lambda=  1.9856 -5.0127I
k=08  x=[  0.0862 -0.2643I,  0.9629 -0.2697I],    lambda=  2.0058 -5.0055I
k=09  x=[  0.0861 -0.2636I,  0.9629 -0.2698I],    lambda=  2.0021 -4.9974I