Power Method MATLAB Program

Power Method, used in mathematics and numerical methods, is an iteration method to compute the dominant eigenvalue and eigenvector of a matrix. It is a simple algorithm which does not compute matrix decomposition, and hence it can be used in cases of large sparse matrices. Power method gives the largest eigenvalue and it converges slowly.

In earlier tutorials, we discussed algorithm/flowchart and C program for power method. Here, we are going to write a program source code for Power method in MATLAB and go through its theoretical background along with a numerical example.

Theoretical Background of Power Method:

In linear algebra, eigenvalue of a square matrix is defined as a vector which points invariant direction under the associated linear transformation. If λ be the eigenvalue of a square matrix A, it can be mathematically depicted as:

Av = λv

where, ‘v’ is a vector (called eigenvector) which is non-zero. This vector ‘v’ is said to be normalized if its coordinate of largest magnitude is unity.

Theorem Derivation:

Consider a matrix A of order mxn which has n distinct eigenvalues λ₁, λ₂, λ₃,… λ_n. The descending order of these eigenvalues can be given as:

|λ₁| > |λ₂| > |λ₃| …… > |λ_n|

Choose an appropriate value of x₀ . Now, the sequence can be given as:

{ xk = (x₁^(k), x₂^(k), x₃^(k), ………. x_n^(k))^T}

The sequence of  {C_k}  is generated by:

Y_n = AX_k and X_k+1 = 1/(C_k+1) * Y_k

where, C_k+1 = X_j^(k) and X_j ^(k) = max{|xi(k)|

These values will converge to the dominant eigenvalue λ₁ and eigenvector v_1.

Alternative Approach:

Power Method uses continuous guesses of X to the eigenvector and λ to the eigenvalue. In both manual calculation and programming approach (MATLAB code) of Power method, the matrix must be symmetric i.e. aij = aji.

Now, again consider a matrix A, whose eigenvector X and eigenvalue λ are to be found out, such that:

[A]{X} = λ{X}

Say, we’re given these three set of equations:

a11x1 + a12x2 + a13x3 = λx1
a12x1 + a22x2 + a23x3 = λx2
a31x1 + a32x2 + a33x3 = λx3

These equations can be re-written as:

(a11-λ)x1 +a12x2 +a13x3 = 0
a21x1 +(a22-λ)x2 +a23x3 = 0
a31x1 +a32x2 + (a33-λ)x3 = 0

Hence, the determinant of the 3*3 matrix thus formed out of the coefficients of ‘x’ terms gives three roots: λ1, λ2 and λ3. These roots are the required eigenvalues or characteristic values.

And, for each of these values, we obtain a set of column vector with elements x1, x2 and x3. This vector is the required eigenvector. This technique is utilized in the below program for Power method in MATLAB.

Did you know??

Google uses power iteration method to calculate the PageRank of webpages in their search engine! Also, the recommendations like “who to follow” in Twitter are based on this method!

Power Method in MATLAB:

function [ v d] = power_method( A )
% for finding the largest eigen value by power method

disp ( ' Enter the matrix whose eigen value is to be found')

% Calling matrix A
A = input ( ' Enter matrix A :   \n')
% check for matrix A
% it should be a square matrix
[na , ma ] = size (A);
if na ~= ma
    disp('ERROR:Matrix A should be a square matrix')
    return
end

% initial guess for X..?
% default guess is [ 1 1 .... 1]'
disp('Suppose X is an eigen vector corresponding to largest eigen value of matrix A')
r = input ( 'Any guess for initial value of X? (y/n):   ','s');
switch r 
    case 'y'
        % asking for initial guess
    X0 = input('Please enter initial guess for X :\n')
        % check for initial guess
    [nx, mx] = size(X0);
    if nx ~= na || mx ~= 1
        disp( 'ERROR: please check your input')
        return
    end
    otherwise
    X0 = ones(na,1);
end

%allowed error in final answer
t = input ( 'Enter the error allowed in final answer:  ');
tol = t*ones(na,1);
% initialing k and X
k= 1;
X( : , 1 ) = X0;
 %initial error assumption
err= 1000000000*rand(na,1);
% loop starts
while sum(abs(err) >= tol) ~= 0
    X( : ,k+ 1 ) = A*X( : ,k); %POWER METHOD formula
    % normalizing the obtained vector 
    [ v i ] = max(abs(A*X( : ,k+ 1 )));
    E = X( : ,k+ 1 );
    e = E( i,1);
    X(:,k+1) = X(:,k+1)/e;
    err = X( :,k+1) - X( :, k);% finding error
    k = k + 1;  
end

%display of final result
fprintf (' The largest eigen value obtained after %d itarations is  %7.7f  \n', k, e)
disp('and the corresponding eigen vector is ')
X( : ,k)

The above code for power method in MATLAB is used to calculate the eigenvalue and eigenvector of a square matrix of any order by using iteration principle of power method. In this program, the matrix whose eigenvalue is to be determined is the input and its corresponding eigenvalue and eigenvector are the output.

As the program is executed, it asks for the value of matrix whose eigenvalue and eigenvector are to be found out. After that, the input matrix is checked whether it is square or not. If the matrix is not a square matrix, an error message is displayed. Further, the matrix should be symmetric or positive definitive.

Then, an initial guess value is to be given to the code, and based on that value and following the theory explained above, the program calculates eigenvalue and eigenvector of the entered matrix. Here’s a sample output of the program:

Power Method Numerical Example:

Let’s now analyze numerically the above program for power method in MATLAB with this example. We have to find out the dominant eigenvalue of the following matrix using Power method:

First, assume an initial iteration vector as:

Now, the following approximation can be obtained:

Thus, the eigenvalue of the given matrix is 190 and the eigenvector is [3 1].

If you have any questions regarding Power method or its MATLAB code, bring them up the comments section. You can find more Numerical Methods tutorial using MATLAB here.