iterative estimation of eigenvalues and eigenvectors by the inverse power method. Your MATLAB/Octave code shouldn’t be hardcoded. ● Calculate the eigenvectors and eigenvalues

function [lambda, v] = inverse_power_method(A, tol, max_iter)

% This function uses the inverse power method to estimate the smallest

% eigenvalue and corresponding eigenvector of matrix A.

% Inputs:

% A - The input square matrix

% tol - The tolerance for convergence (e.g., 1e-6)

% max_iter - The maximum number of iterations

% Outputs:

% lambda - The estimated smallest eigenvalue

% v - The estimated corresponding eigenvector

n = size(A, 1); % Get the size of the matrix

v = rand(n, 1); % Initialize a random vector

v = v / norm(v); % Normalize the vector

% Perform the inverse power iteration

for k = 1:max_iter

w = (A \ v); % Solve A * w = v

v_new = w / norm(w); % Normalize the vector

% Check for convergence

if norm(v_new - v) < tol

break;

end

v = v_new; % Update the vector

end

% Calculate the corresponding eigenvalue

lambda = v' * A * v;

% Normalize the eigenvector to have unit length

v = v / norm(v);

end

% Define the matrix A

A = [1, -3, 3;

3, -5, 3;

6, -6, 4];

% Set the tolerance and maximum number of iterations

tol = 1e-6;

max_iter = 1000;

% Call the inverse power method function

[lambda, v] = inverse_power_method(A, tol, max_iter);

% Display the results

fprintf('Estimated smallest eigenvalue: %.6f\n', lambda);

fprintf('Estimated corresponding eigenvector:\n');

disp(v);

[John D'Errico]

Exactly what is your question?

Answer 1

[John D'Errico]

Let me guess, you have a problem with this matrix.

Does the matrix have any special properties? What, for example, are the eigenvalues of A?

A = [1, -3, 3;
     3, -5, 3;
     6, -6, 4];
[V,D] = eig(A)

V = 3x3

-0.4082 -0.8103 0.1933 -0.4082 -0.3185 -0.5904 -0.8165 0.4918 -0.7836

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

D = 3x3

4.0000 0 0 0 -2.0000 0 0 0 -2.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Do you see the smallest eigenvalue is actually a double one, and it is negative. Negative is not an issue. But the higher multiplicity throws a curve at you. And, while your code does work, returning -2 as the eigenvalue, it did not generate the same eigenvector as one of those from eig?

[lambda, vec] = inverse_power_method(A, 1e-6,10000)

lambda = -2

vec = 3x1

0.5289 0.8031 0.2742

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Now you need to read your class notes again. I assume it was in there. Well, maybe not. Or not yet. This may be in the notes from tomorrow. You should have learned that if there is an eigenvalue of multiplicity higher than 1, then the eigenvectors form a subspace. So the eigenvector you will find will generally be some linear combination of the two eigenvectors eig would return. HOWEVER, this matrix is what we would call a "defective" matrix. That may have been a mean thing to do to you as homework, to give you such a matrix, and then have you try to figure out what is wrong. If that is the case, then the eigenvectors of A will not form a complete basis of linearly independent vectors. Do some reading here:

https://en.wikipedia.org/wiki/Defective_matrix#:~:text=In%20linear%20algebra%2C%20a%20defective,linearly%20independent%20eigenvectors.

We can do this, and if the product V'*V is not effectively an identity matrix, then there is an issue with A.

V'*V

ans = 3x3

1.0000 0.0593 0.8020 0.0593 1.0000 -0.3540 0.8020 -0.3540 1.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Ah yes. So A is in fact a defective matrix. But is lambda an eigenvalue? Yes.

[A*vec,lambda*vec]

ans = 3x2

-1.0578 -1.0578 -1.6063 -1.6063 -0.5484 -0.5484

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

And we see that vec as found does have the property of being an eigenvector, lambda an eigenvalue. But if you re-run your code, you will get a totally different eigenvector, at random. And neither of them will match the vectors that will be returned from eig. You did start with a random vector, so that is not surprising, if the eigenvector is not unique.

Does this resolve your question? As I said, it may have been a slightly mean trick to throw such a matrix at you, without giving you some background for that class of matrix, and what to expect as a result. Or maybe it is in your notes.

The classic example of a defective matrix is this one:

Adef = [1 1; 0 1];

Clearly the eigenvalues are 1 and 1.

[V,D] = eig(Adef)

V = 2x2

1.0000 -1.0000 0 0.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

D = 2x2

1 0 0 1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

You can see the eigenvectors returned from eig are not at all orthogonal.

We get spoiled, and so often expect all matrices to be simple things. For example...

Asimpl = rand(3);Asimpl = Asimpl + Asimpl'

Asimpl = 3x3

0.6317 0.7976 0.7983 0.7976 0.9042 0.2476 0.7983 0.2476 0.0954

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

So Asimpl is a nice, well-behaved symmetric matrix with distinct eigenvalues, and your code has no problems...

[V,D] = eig(Asimpl)

V = 3x3

-0.6504 0.3736 0.6613 0.2349 -0.7290 0.6429 0.7223 0.5735 0.3864

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

D = 3x3

-0.5429 0 0 0 0.3008 0 0 0 1.8735

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[lambda, v] = inverse_power_method(Asimpl,1e-8,1000)

lambda = 0.3008

v = 3x1

0.3736 -0.7290 0.5735

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

function [lambda, v] = inverse_power_method(A, tol, max_iter)
% This function uses the inverse power method to estimate the smallest
% eigenvalue and corresponding eigenvector of matrix A.
% Inputs:
%   A        - The input square matrix
%   tol      - The tolerance for convergence (e.g., 1e-6)
%   max_iter - The maximum number of iterations
% Outputs:
%   lambda   - The estimated smallest eigenvalue
%   v        - The estimated corresponding eigenvector
n = size(A, 1);  % Get the size of the matrix
v = rand(n, 1);  % Initialize a random vector
v = v / norm(v); % Normalize the vector
% Perform the inverse power iteration
for k = 1:max_iter
    w = (A \ v);     % Solve A * w = v
    v_new = w / norm(w); % Normalize the vector
    % Check for convergence
    if norm(v_new - v) < tol
        break;
    end
    v = v_new; % Update the vector
end
% Calculate the corresponding eigenvalue
lambda = v' * A * v;
% Normalize the eigenvector to have unit length
v = v / norm(v);
end