1.	$d (1) = 0$ . This setting should be used if $λ$ is an ill-conditioned eigenvalue (provided the matrix elements do not vary widely in order of magnitude). In this case it is essential to accept only a vector found after one half iteration, and $μ$ must be a very good approximation to $λ$ . If acceptable growth is achieved in the solution of $U x = e$ , then the normalized $x$ is accepted as the eigenvector. If not, columns of an orthogonal matrix are tried in turn in place of $e$ . If none of these give acceptable growth, the function fails, indicating that $μ$ was not a sufficiently good approximation to $λ$ .
2.	$d (1) > 0$ . This setting should be used if $μ$ is moderately close to an eigenvalue which is not ill-conditioned (provided the matrix elements do not differ widely in order of magnitude). If acceptable growth is achieved in the solution of $U x = e$ , the normalized $x$ is accepted as the eigenvector. If not, inverse iteration is performed. Up to $30$ iterations are allowed to achieve a vector and a correction to $μ$ which together give acceptably small residuals.
3.	$d (1) < 0$ . This setting should be used if the elements of $A$ and $B$ vary widely in order of magnitude. Inverse iteration is performed, but a different convergence criterion is used.

See Algorithmic Details for further details.

Note that the bandwidth of the matrix

A

must not be less than the bandwidth of

B

. If this is not so, either

A

must be filled out with zeros, or matrices

A

and

B

may be reversed and

1 / μ

supplied as an approximation to the eigenvalue

1 / λ

. Also it is assumed that

A

and

B

each have the same number of subdiagonals as superdiagonals. If this is not so, they must be filled out with zeros. If

A

and

B

are both symmetric, only the upper triangles need be supplied.

References

Peters G and Wilkinson J H (1979) Inverse iteration, ill-conditioned equations and Newton's method SIAM Rev. 21 339–360

Wilkinson J H (1965) The Algebraic Eigenvalue Problem Oxford University Press, Oxford

Wilkinson J H (1972) Inverse iteration in theory and practice Symposia Mathematica Volume X 361–379 Istituto Nazionale di Alta Matematica, Monograf, Bologna

Wilkinson J H (1974) Notes on inverse iteration and ill-conditioned eigensystems Acta Univ. Carolin. Math. Phys. 1–2 173–177

Wilkinson J H (1979) Kronecker's canonical form and the

Q Z

algorithm Linear Algebra Appl. 28 285–303

Parameters

Compulsory Input Parameters

1: $ma1$ – int64int32nag_int scalar

The value

m_{A} + 1

, where

m_{A}

is the number of nonzero lines on each side of the diagonal of

A

. Thus the total bandwidth of

A

2 m_{A} + 1

Constraint:

1 \leq ma1 \leq n

2: $mb1$ – int64int32nag_int scalar

mb1 \leq 0

B

is assumed to be the unit matrix. Otherwise mb1 must specify the value

m_{B} + 1

, where

m_{B}

is the number of nonzero lines on each side of the diagonal of

B

. Thus the total bandwidth of

B

2 m_{B} + 1

Constraint:

mb1 \leq ma1

3: $a (lda, n)$ – double array

lda, the first dimension of the array, must satisfy the constraint

lda \geq 2 \times ma1 - 1

The

n

n

band matrix

A

. The

m_{A}

subdiagonals must be stored in the first

m_{A}

rows of the array; the diagonal in the (

m_{A} + 1

)th row; and the

m_{A}

superdiagonals in rows

m_{A} + 2

2 m_{A} + 1

. Each row of the matrix must be stored in the corresponding column of the array. For example, if

n = 6

and

m_{A} = 2

the storage scheme is:

\begin{array}{l} * & * & a_{31} & a_{42} & a_{53} & a_{64} \\ * & a_{21} & a_{32} & a_{43} & a_{54} & a_{65} \\ a_{11} & a_{22} & a_{33} & a_{44} & a_{55} & a_{66} \\ a_{12} & a_{23} & a_{34} & a_{45} & a_{56} & * \\ a_{13} & a_{24} & a_{35} & a_{46} & * & * \end{array} .

Elements of the array marked

*

need not be set. The following code assigns the matrix elements within the band to the correct elements of the array:

 for j=1:n
 for i=max(1,j-ma1+1):min(n,j+ma1-1)
 a(i-j+ma1,j) = matrix(i,j);
 end
 end

sym = true

(i.e., both

A

and

B

are symmetric), only the lower triangle of

A

need be stored in the first ma1 rows of the array.

4: $b (ldb, n)$ – double array

ldb, the first dimension of the array, must satisfy the constraint

if $sym = false$ , $ldb \geq 2 \times mb1 - 1$ ;
if $sym = true$ , $ldb \geq mb1$ .

mb1 > 0

, b must contain the

n

n

band matrix

B

, stored in the same way as

A

. If

sym = true

, only the lower triangle of

B

need be stored in the first mb1 rows of the array.

mb1 \leq 0

, the array is not used.

5: $sym$ – logical scalar

sym = true

, both

A

and

B

are assumed to be symmetric and only their upper triangles need be stored. Otherwise sym must be set to false.

6: $relep$ – double scalar

The relative error of the coefficients of the given matrices

A

and

B

. If the value of relep is less than the machine precision, the machine precision is used instead.

7: $rmu$ – double scalar

μ

, an approximation to the eigenvalue for which the corresponding eigenvector is required.

8: $d (30)$ – double array

d (1)

must be set to indicate the type of problem (see Description):

$d (1) > 0.0$: Indicates a well-conditioned eigenvalue.
$d (1) = 0.0$: Indicates an ill-conditioned eigenvalue.
$d (1) < 0.0$: Indicates that the matrices have elements varying widely in order of magnitude.

Optional Input Parameters

1: $n$ – int64int32nag_int scalar: Default: the second dimension of the arrays a, b. (An error is raised if these dimensions are not equal.)
$n$ , the order of the matrices $A$ and $B$ .

Constraint: $n \geq 1$ .

Output Parameters

1: $a (lda, n)$ – double array: Details of the factorization of $A - \bar{λ} B$ , where $\bar{λ}$ is an estimate of the eigenvalue.
2: $b (ldb, n)$ – double array: Elements in the top-left corner, and in the bottom right corner if $sym = false$ , are set to zero; otherwise the array is unchanged.
3: $vec (n)$ – double array: The eigenvector, normalized so that the largest element is unity, corresponding to the improved eigenvalue $rmu + d (30)$ .
4: $d (30)$ – double array: If $d (1) \neq 0.0$ on entry, the successive corrections to $μ$ are given in $d (i)$ , for $i = 1, 2, \dots, k$ , where $k + 1$ is the total number of iterations performed. The final correction is also given in the last position, $d (30)$ , of the array. The remaining elements of d are set to zero.
If $d (1) = 0.0$ on entry, no corrections to $μ$ are computed and $d (i)$ is set to $0.0$ , for $i = 1, 2, \dots, 30$ . Thus in all three cases the best available approximation to the eigenvalue is $rmu + d (30)$ .
5: $ifail$ – int64int32nag_int scalar: $ifail = 0$ unless the function detects an error (see Error Indicators and Warnings).

Error Indicators and Warnings

Errors or warnings detected by the function:

$ifail = 1$

On entry,	$n < 1$ ,
or	$ma1 < 1$ ,
or	$ma1 > n$ ,
or	$lda < 2 \times ma1 - 1$ ,
or	$ldb < mb1$ when $sym = true$ ,
or	$ldb < 2 \times mb1 - 1$ when $sym = false$ (ldb is not checked if $mb1 \leq 0$ ).

$ifail = 2$

On entry,

ma1 < mb1

. Either fill out a with zeros, or reverse the roles of a and b, and replace rmu by its reciprocal, i.e., solve

B x = λ^{- 1} A x .

$ifail = 3$

On entry,	$lwork < 2 \times n$ when $d (1) = 0.0$ ,
or	$lwork < n \times (ma1 + 1)$ when $d (1) \neq 0.0$ .

$ifail = 4$: $A$ is null. If $B$ is nonsingular, all the eigenvalues are zero and any set of n orthogonal vectors forms the eigensolution.

$ifail = 5$: $B$ is null. If $A$ is nonsingular, all the eigenvalues are infinite, and the columns of the unit matrix are eigenvectors.

$ifail = 6$

On entry,

A

and

B

are both null. The eigensolution is arbitrary.

$ifail = 7$: $d (1) \neq 0.0$ on entry and convergence is not achieved in $30$ iterations. Either the eigenvalue is ill-conditioned or rmu is a poor approximation to the eigenvalue. See Algorithmic Details.

$ifail = 8$: $d (1) = 0.0$ on entry and no eigenvector has been found after $\min (n, 5)$ back-substitutions. rmu is not a sufficiently good approximation to the eigenvalue.

$ifail = 9$: $d (1) < 0.0$ on entry and rmu is too inaccurate for the solution to converge.

$ifail = - 99$: An unexpected error has been triggered by this routine. Please contact NAG.

$ifail = - 399$: Your licence key may have expired or may not have been installed correctly.

$ifail = - 999$: Dynamic memory allocation failed.

Accuracy

The eigensolution is exact for some problem

(A + E) x = μ (B + F) x,

where

‖E‖, ‖F‖

are of the order of

η (‖A‖ + μ ‖B‖)

, where

η

is the value used for relep.

Further Comments

Timing

The time taken by nag_eigen_real_band_geneig (f02sd) is approximately proportional to

n {(2 m_{A} + 1)}^{2}

for factorization, and to

n (2 m_{A} + 1)

for each iteration.

Storage

The storage of the matrices

A

and

B

is designed for efficiency on a paged machine.

nag_eigen_real_band_geneig (f02sd) will work with full matrices but it will do so inefficiently, particularly in respect of storage requirements.

Algorithmic Details

Inverse iteration is performed according to the rule

(A - μ B) y_{r + 1} = B x_{r}

x_{r + 1} = \frac{1}{α_{r + 1}} y_{r + 1}

where

α_{r + 1}

is the element of

y_{r + 1}

of largest magnitude.

Thus:

(A - μ B) x_{r + 1} = \frac{1}{α_{r + 1}} B x_{r} .

Hence the residual corresponding to

x_{r + 1}

is very small if

|α_{r + 1}|

is very large (see Peters and Wilkinson (1979)). The first half iteration,

U y_{1} = e

, corresponds to taking

L^{- 1} P B x_{0} = e

μ

is a very accurate eigenvalue, then there should always be an initial vector

x_{0}

such that one half iteration gives a small residual and thus a good eigenvector. If the eigenvalue is ill-conditioned, then second and subsequent iterated vectors may not be even remotely close to an eigenvector of a neighbouring problem (see pages 374–376 of Wilkinson (1972) and Wilkinson (1974)). In this case it is essential to accept only a vector obtained after one half iteration.

However, for well-conditioned eigenvalues, there is no loss in performing more than one iteration (see page 376 of Wilkinson (1972)), and indeed it will be necessary to iterate if

μ

is not such a good approximation to the eigenvalue. When the iteration has converged,

y_{r + 1}

will be some multiple of

x_{r}

y_{r + 1} = β_{r + 1} x_{r}

, say.

Therefore

(A - μ B) β_{r + 1} x_{r} = B x_{r},

giving

(A - (μ + \frac{1}{β_{r + 1}}) B) x_{r} = 0 .

Thus

μ + \frac{1}{β_{r + 1}}

is a better approximation to the eigenvalue.

β_{r + 1}

is obtained as the element of

y_{r + 1}

which corresponds to the element of largest magnitude,

+ 1

, in

x_{r}

. The function terminates when

‖(A - (μ + \frac{1}{β_{r}}) B) x_{r}‖

is of the order of the machine precision relative to

‖A‖ + |μ| ‖B‖

If the elements of

A

and

B

vary widely in order of magnitude, then

‖A‖

and

‖B‖

are excessively large and a different convergence test is required. The function terminates when the difference between successive corrections to

μ

is small relative to

μ

In practice one does not necessarily know if the given problem is well-conditioned or ill-conditioned. In order to provide some information on the condition of the eigenvalue or the accuracy of

μ

in the event of failure, successive values of

\frac{1}{β_{r}}

are stored in the vector d when

d (1)

is nonzero on input. If these values appear to be converging steadily, then it is likely that

μ

was a poor approximation to the eigenvalue and it is worth trying again with

rmu + d (30)

as the initial approximation. If the values in d vary considerably in magnitude, then the eigenvalue is ill-conditioned.

A discussion of the significance of the singularity of

A

and/or

B

is given in relation to the

Q Z

algorithm in Wilkinson (1979).

Example

Given the generalized eigenproblem

A x = λ B x

where

A = (\begin{array}{r} 1 & 1 & 2 \\ - 1 & 2 & 1 & 2 \\ - 1 & 3 & 1 & 2 \\ - 1 & 4 & 1 \\ - 1 & 5 \end{array}) and B = (\begin{array}{l} 5 & 1 \\ 1 & 4 & 2 \\ 2 & 3 & 2 \\ 2 & 2 & 1 \\ 1 & 1 \end{array})

find the eigenvector corresponding to the approximate eigenvalue

- 12.33

Although

B

is symmetric,

A

is not, so sym must be set to false and all the elements of

B

in the band must be supplied to the function.

A

(as written above) has

1

subdiagonal and

2

superdiagonals, so ma1 must be set to

3

and

A

filled out with an additional subdiagonal of zeros. Each row of the matrices is read in as data in turn.

Open in the MATLAB editor: f02sd_example

function f02sd_example


fprintf('f02sd example results\n\n');

% A and B in non-symmmetric banded format:
% A has 2 super-diagonals, B has 1 sub/super-diagonal;
% a row of zeros is added to A to give same number of subdiagonals.
a = [0, 0, 0, 0, 0;
     0,-1,-1,-1,-1;
     1, 2, 3, 4, 5;
     1, 1, 1, 1, 0;
     2, 2, 2, 0, 0];
b = [0, 1, 2, 2, 1;
     5, 4, 3, 2, 1;
     1, 2, 2, 1, 0];

ma1 = int64(2+1);
mb1 = int64(1+1);
sym = false;

% Find eigenvalue closest to rmu and associated eigenvector.
rmu   = -12.33;
relep = 0;
d     = zeros(30, 1);
% Assume the eigenvalue is well-conditioned
d(1)  = 1;

[a, b, vec, d, ifail] = f02sd( ...
                               ma1, mb1, a, b, sym, relep, rmu, d);

eval = rmu + d(30);
disp('Corrected eigenvalue');
disp(eval);
disp('Eigenvector')
disp(vec/norm(vec));

f02sd example results

Corrected eigenvalue
  -12.3394

Eigenvector
   -0.0377
    0.2606
   -0.5560
    0.6598
   -0.4315

PDF version (NAG web site, 64-bit version, 64-bit version)

Chapter Contents

Chapter Introduction

NAG Toolbox