- Email: [email protected]
A method for the solution of the complete eigenvalue problem
A METHOD FOR THE SOLUTION OF THE COMPLETE EIGENVALUE PROBLEM+ V. V. VOYEVODIN Moscow (Received 22 September 1961)
IN THIS paper we shall give some iterative methods for the solution of the complete eigenvalue problem of a matrix. These can be easily carried out on electronic computers and the main portion of the work can be done with fixed-point arithmetic without materially reducing the accuracy. The method for solving the complete eigenvalue problem for normal matrices is given in Q 1. It resembles Jacobi’s process in its structure (see [I], 5 Sl). The method of solution in the case of an arbitrary matrix with a simple structure is given in 5 2. It resembles the triangular power method (see [l], 9 78). For simplicity we shall consider real matrices only, although all the methods given can be applied to complex matrices also without particular difficulty. fr.I. METHOD OF SOLVING THE COMPLETE EIGENVALUE A NORMAL
PROBLEM FOR
MATRIX
1. Determination of’ the principal values oj’ a matrix We shall first give a brief account of the method of determining the principal vaIues of a matrix (see [2]), because this will be used in solving the complete problem of eigenvalues of a normal matrix. (Certain changes have been introduced in the description of this method.) The principal values of the matrix A are the positive values of the square root of the eigenvalues of the matrix A*A (or AA*). The method is based on the following result. Let U and V be two orthogonal matrices such that II = UAV, where II is a diagonal matrix. Then the moduli of the diagonal matrix II will give the principal values of the matrix A, the rows of the eigenvectors of the matrix AA*, and the columns of the matrix vectors of the matrix A*A. All this follows from the two obvious ing from (1): II II* = UAA*U*,
I-I*II = V*A*AV
The computing scheme for the method is as follows. 7 Zh. vych. math. 2: No. I, 15-24, 1962.
13
(1) elements the matrix V are the equations
of the U are eigenresult(2)
14
V. V. VOYEVODIN
In the matrix A we choose a pair of non-diagonal elements a, and a,, such that aTj+uji
=
ytlx
and construct the two simple rotation
(&+a&)#
matrices:
cos p
. ..-sin v . ..
u&P =
(3) 1
/1. ..cos y . .
VW = 11
I;
sinv...
i’l
. ..i
i *.
--sinv . .
\
0
cos I .. .
j
. ..I
‘** 1 /
\O
.
The angles 47and y will be defined by the condition that the elements a,j and uji of the matrix Al = U$)AVij) become zero. This condition gives a system of two equations for q and y: COS?/J(Uij
COS lJ?+Ujj
sin y)-siny
(UiiCOS
q+Uji
sin $0)= 0, (4)
COS?#(--U~ising,+UjiCOSpl)+sin~(U~jsinpl+UjjCOS~)=0, or COSQl(UjiCOS
y+Ujj
sin y) -sin
Q1(UiiCOS
sin y) = 0,
y+Uij
(5)
COSpl(-U~iSin?fJ+UijCOSy)+sing,(-Ujisin~+UjjCOS~)=O. These
systems
of. equations In order
algebraic
systems
with respect to cos y, sin y and cos p, sin p with a normalized
can
be considered
solution.
that we may have non-zero
the determinants
of the systems
is easy to obtain formulae
as linear solutions
should
become
for the determination UiiajifUijajj
tan 29 = 2 ___--___.
&-ajj+a:j-aji'
tan
homogeneous it is necessary zero.
From
and sufficient that these
conditions
it
of the angles p and y:
aii aij+ajiajj 2y = 2 2
2.
(6)
aii-ajj+U,Ti-Uij
Those values of cp and y are taken which are determined from (6) and satisfy the system (4) or (5). The matrix A, is constructed after U$) and V$) a similar process being used. The required matrices U and V are determined as infinite products of all the matrices Uiy) and VjT) (m = 1, 2, . . .) respectively, obtained in each step. All theorems on the convergence of this method are proved in exactly the same way as in Jacobi’s method (see [I]).
A method for the solution of the complete eigenvalue problem
1.5
We now pass directly to the solution of the complete eigenvalue problem of normal matrices. The matrix A is called normal if it can be commuted with its conjugate, i.e. AA* = A*A.
(7)
The following proposition was proved for normal matrices (see [3]): in order that a matrix A of’the n-th order may be normal it is necessary and suficient that it should have orthogonal eigenveetors. It follows from this proposition that a normal matrix A can always be represented
in the form A = QAQ*,
(8)
where A is a diagonal and Q is a unitary matrix. Hence A*A =QA*Q*QAQ’
= QA*AQ*;
(9)
therefore the principal values of a normal matrix are equal to the moduli of the eigenvalues of the matrix itself. This and the fact that the principal values can be determined by the method described above make it possible to solve effectively the complete eigenvalue problem for a normal matrix. We consider a few cases of the distribution of eigenvalues. 2. All the eigenvalues of a normal matrix A are simple and di@er in their moduli, except the complex conjugate values
In this case we can obtain an expansion of the matrix A into factors by the method of determination of the principal values (see [I]): A = U’DV’,
(10)
The matrices U and V can always be chosen in such a way that diagonal elements of the matrix D with equal moduli are situated side by side and have identical signs. It follows from (10) that the matrix A will be similar to the matrix D = V’AV = V’U’IIV’V = V’U’II = RD.
(11)
Let us investigate the form of the matrix D. Since the matrix A is normal, the matrix D will also be normal, i.e. the equation DD’ = D’D. will be satisfied. Since the matrix R is orthogonal
(12) we obtain from (12)
RII (RII)’ = (RTI)’ RI-I,
(13)
or RIIII’R
= D’R’RD.
Hence RIP = IPR.
(14)
16
V. V. VOYEVODIN
Let rij be the elements of the matrix R andpFi the diagonal elements of the matrix I12. Then from (13) we obtain that for any pair of indices (k, r) the relation rk; 6
=
dkrklr
WJ)
will be satisfied. Hence if p& # ~7~then rkr = 0. From the assumption regarding the distribution of the eigenvalues of the matrix A it follows that the multiplicity of the diagonal elements of the matrix II2 is not greater than 2. Therefore it is evident that the matrix R will be quasi-diagonal with sub-matrices whose diagonals are not higher than the second order. The matrix D will also have the same form: 0 DI D=
D2.
(\O
(16)
I'hn
Thus
the solution of the complete problem for the matrix A is reduced to the solution of the complete problem of the matrix D which has a very simple form. Let Dj be a sub-matrix of the second order. Then its matrix is equal to the product of an orthogonal and a diagonal matrix:
Since sub-matrices of the second order correspond to complex-conjugate values, therefore ($
$)
eigen-
(17)
= (I;;;-::;;).
Hence we have the final expression for the eigenvalues ;Ifi and eigenvectors zei of the matrix with the sub-matrix Dj: n$jg= pjj(cos~j&isin~j),
z(j) 1.2--
+=; 1.
+J.
Thus with the assumption made regarding the distribution of the eigenvalues of a normal matrix the complete eigenvalue problem is fully solved by using the method of principal values once. If the matrix has more than two eigenvalues with the same modulus, sub-matrices of an order higher than the second may appear in the matrix D. For orthogonal matrices, in particular, this method will not make the solution of the problem any simpler. 3. All the eigenvalues of’ a normal matrix are simple; any number oj them may have the same modulus
Let us represent the matrix A in the form A=B+C,
(19,
A method for the solution of the complete eigenvalue problem where the matrix B is symmetrical and C is skew-symmetrical. tively
the elements
Jacobi’s
aij, bij, cij are respec-
of the matrices A, B, C, then b,
Using
17
=
bji =
method
a-
;
Cij
-cji
=
C
aij-ajj 2 -_
(see [l]) we reduce the symmetrical
matrix B to the
diagonal matrix% and the matrix C to a similar matrix ii,. Then for the matrix A we shall have A = Q(ii+i?)Q’, where Q is an orthogonal
transformation
Since the matrix A is normal the matrix B+6 have the equation (B-j-C)(B+C)’
= (B tZ$‘(BS_C),
W)
matrix.
or
will also be normal. Hence we
(i&j-C)(El-cc),=
(B--C)(B+C),
Hence SC = CB.
(22)
As in the case of eqn. (14) it can be inferred that the matrix B+ecan be reduced to the quasi-diagonal form (16) by transposing the rows and columns, only the sub-matrices Dj(jj = 1, 2, . . . , m) in this case will have an arbitrary order lj, but these must be of the form: Dj 1~fDj-tdjE,
(23)
where aj is a skew-symmetrica matrix, and dj is a real number. The eigenvaiues A$$ (k z- 1,2 , . . . , lj) of the matrix Dj are retated to the eigenvalues @ of the matrix 6j by the equation @‘=
$‘+&
J’
(24)
and their eigenvectors coincide so that further solution of the complete problem for normat matrices is reduced to a solution of the complete problem for skewsymmetrical matrices By hypothesis the matrix A has no multiple eigenvalues; the skew-symmetrical 5j.
matrix fij also will not have any. Hence the method described in section 2 can be successfully used for the solution of the complete eigenvalue problem of the latter niatrix. A special case is of particular interest. Let the matrix A be orthogonal and have no multiple eigenvalues. In that case the sub-matrices Dj(j = 1, 2, . . . . m) will also be orthogonal matrices. It is known (see [3]) that all the eigenvalues of an orthogonal matrix have a modulus equa1 to unity and all eigenvalues of a skew-symmetrical matrix are pure imaginary numbers. Hence from equation (24) it follows J$‘!‘+$=
I,
k=
1,2, . . ..ij.
G5)
will be satisfied. Since there can be no more than two different pure imaginary numbers @ satisfying the equation
(29, it follows that none of the sub-matrices
v. v.
18
VoYEvoDIN
ofD,(j = 1,2, . . . m) will be of an order higher than the second. Thus if an orthogonal matrix has no multiple eigenvalues the solution of its complete eigenvalue problem can be obtained by using Jacobi’s method once. Finally it should be noted that when this method is used in practice the expansion (19) is unnecessary, and the orthogonal transformations should be applied with it to the matrix A after calculating the angle of rotation from the formula afj+aji tan29 = ~ #ii-ajj
*
These transformations are equivalent to reducing the matrix B by Jacobi’s method to the diagonal form, but the matrix A can be reduced directly to the matrix B+C. To ensure maximum speed of convergence at each step the pair of indices (i,j) should be so chosen that the sum Gfj+aji has the m~imum modulus. 4. The normal matrix A has multiple eigenvalues
It can be easily seen that if the method described in section 3 is used in this case, we shall ultimately have to solve the complete problem for matrices N which are both skew-symmetrical and orthogonal, i.e. for which N’=
N’N=NN’=E.
-N,
(27)
from (27) we easily get NZ= -E.
(28)
The eigenvalues of the matrix N are equal to ki, and the eigenvectors can be easily found by remembering the following. Let p be an arbitrary non-zero vector, then the vector p-iNp will be an eigenvector of the matrix N corresponding to an eigenvalue equal to +i. Bearing (28) in mind we have N(p-iNp)
= Np--iNap
0 2. METHOD OF SOtVING
=Np+ip
= i(p-iip).
THE COMPLETE EIGENVALUE ARBITRARY MATRIX
(29)
PROBLEM OF AN
Let a real non-specialized matrix A be given, corresponding to an operator which is real in n-dimensional space R and has n linear independent eigenvectors Xl. . * . . x,. Let us divide the eigenvalues Ai (i = 1,2, . . . , n) of the matrix A into m groups, each group containing only eigenvalues with the same modulus, and suppose that the condition (30) is satisfied, where k, = n. We construct m sub-spaces PX, Pz, . . . . P,,, of the space R invariant with respect to A by defining Pj G = 1,2, . . . , m) as a sub-space formed by the eigenvectors x1, . ..*Q. Evidently the following relations hold (31)
A method for the solution of the complete eigenvalue problem
19
with the vectors ul, uZ, . . . , u,, form a basis of the space R having the property that the VeCtOrS II,, . . . , Ukj (j = 1, 2, . . . , m), form a basis of the sub-space Pj (j = 1, 2, .** , m). Let us determine the form of the matrix of an operator given on this base, (see [l], 5 7). Since the vectors AU,, . . . AUkI (j = 1, 2, . . . , m) belong to Pj i.e. are linear combinations of the vectors ul, . . . , ut, only, their coordinates on the base selected beginning from kj+ 1 are equal to zero. Consequently the matrix A on the new base has a quasi-triangular form AllA,, . . . A,, A=
0 A22...A2m . . . . . +. . . . . . . . . (0
(32)
0 e.. A,,,,,,1.
where the sub-matrix Ajj(j = 1,2, . . . , m) has the order kj-kj_1, k. = 0. If Q is a matrix the columns of which are respectively the coordinates of the vectors ul, . . . u, given on the old base, the matrices A and A are connected by the relation A===Q-‘AQ.
(33)
The problem of finding the eigenvalues of the matrix A is considerably simpler than for A, because it amounts to finding the eigenvalues of the matrices Ajj (j = 1,2, . . . ) m), which are usually of a low order. In particular, if the eigenvalues of the matrix A all have different moduli, the sub-matrices Ajj(j = 1,2, . . . , m) will be of the first order. Thus the solution of the complete eigenvalue problem of the matrix A or its simplification when there are eigenvalues with muftiple moduli is reduced to finding a matrix Q, which satisfies the eqn. (33) where A is the matrix in (32). The matrix Q is found as follows. An arbitrary orthogonal matrix Q. is taken and the sequence of orthogonal matrices Qk(k = 1, 2, . . . ,) is constructed by a method similar to the triangular power method;
AQo=Q,A,, AQ,=Qa&, . . . . . . . . *..a . . . . . . . . .
(34)
AQc=Qk+~h+~r Here A, (k = 1,2, .a.) are right triangular matrices. It is best to expand the matrix into a product of orthogonal and right triangular matrices by the reflection method (see PI, § 16). We shall prove that the sequence of matrices Q,(k = 1,2, . , .) has the same property as the sequence of matrices Ak = Q;AQk tend to a quasi-triangular matrix of the form (32). To do this we have to evaluate the minors of the matrix & = Q;AtkAkQO.
(35)
Let A = QOCAC-“Q;,
R=[l,,...,J,]
(36)
V. V.
20
VOYEVODIN
then B, = C-l’AkC’CAkC-l
(37)
are introduced for the minors of the matrix B,:
the following abbreviations
bfk! btk! se,1 ‘a,2
b:,k,!, ......... ... ......
b!k! bfk!
llJa.*.
1111
+.-
(38)
btk> ‘IJI
Using Binet Cauchy’s formula (see [3]) for the determinant of the product of two rectangular matrices, we obtain Bk(;;;;:
::: j
= I_~~~~~,cn,_~C-l~AkC’ xCnkC-l
=
c
I
x
C-1’
c
(;;,;:,“.‘.:f’,,)
(3g)
x
Pl,Pa, ...>Pr (h,j2, .. ..A ) (;I:
:/Jfq;:;
:::; ;;)“‘(“b,;
:::;,“:)
x
l~BIC...~B~~”
,_.c,, . .
’
in q-I:: 1’
::::Epk(:::
:p’(:::::::;)*
::::
since
we shall have, taking into account (30), (40) where O$j are quantities bounded above with respect to the modulus for all i, j, k, I by a constant independent of k. The expression for principal minors of the matrix B, is Bk(;:;::::::) =
x l
[
c
c-l(;:~.:::;“‘)c~::
:::J;)A:‘“:“...XL
]
l~al~...G+” =
pppp
. . . pp@p.
(41)
The matrix Q,, can always be so chosen that all the quantities Oj”) will be bounded below for all k by a strictly positive number independent of k. From (34) we obtain AkQ, = Qk’bk’&-I...&
= Qk.Vk,
(42)
whence B=V;Vk.
(43)
A method
for the solution
of the complete
eigenvalue
21
problem
Let 8$) be the elements of the matrix Vk. Then we know (see [3]) that
#),.j(&) iJ
j=i+l,i+2
II
,...,
n.
It can be easily deduced from (40) (41) (44) that the matrix bk can be represented in the form Vk= Dk.V/‘,
(45)
where D, = [I1llk,...,lil,[k] and V, is a right triangular matrix with elements bounded above with respect to the modulus for all values of k by a constant independent of k. The diagonal elements of the matrix Vk are bounded below for all values of k by a strictly positive number independent of k. Hence vi’ is a right triangular matrix with elements bounded above with respect to the modulus for all values of k by a constant independent of k. From (42) we obtain Qk = AkQOVkl.
(46)
If we form the matrix A, and take into account (45) and (46), we shall have Ak = Q;AQk = QklAQk = VkQolA-kAAkQ,,Vil = VkQ;‘AQ,,V;l
= DkVkQ,lAQ,V;lD+=
D,N,D;l.
(47)
The elements of the matrix Nk are numbers bounded above with respect to the modulus for all values of k by a constant independent of k. Let us divide up the matrix A, into submatrices Ai;) as in (32) and let Zj(j = 1, 2, . . ., m) be a number equal to the modulus of the eigenvalues in the jth group. Then it can be deduced from (47) that the elements of the submatrix A$) for i > j are quantities
Hence the matrices A, tend to a quasi-triangular form. The chief merit of this method is that it is free from errors of rounding off, because orthogonal transformations of matrices are used and this makes it selfcorrecting. But when eigenvalues with nearly equal moduli are present this method gives a very weak convergence. This difficulty can be overcome by using the following analogue of the M-algorithm (see [I], 0. 80). Let AzpQ P V P’
(48)
22
V.V.
where Q, is an orthogonal
VOYEVODlh'
and e,, a right triangular
matrix. With the notations
introduced above we shall evidently have QP = QLp, and fi, = WpP when Q. = E. Squaring equation (48) we obtain A @+I = ijpijptjpi?,.
(4%
If the matrix e7,Q, is expanded as a product of an orthogonal a right triangular matrix Lp, then A ap+1 = i&R,L,i$
matrix R, and
= ~,,,vp,,,
(50)
where ii p+l = QPRP,
and eP+r
= Pep.
(51)
The matrices Q1, Qz, Qa, QB, . . . . Q2p can be determined successively by this method. The method will converge under the conditions of the principal method, the convergence being of the order of
i.e. quadratic. It can be inferred from (45) that, as a rule, the elements of the matrix GP will either increase or decrease rapidly with p. To reduce the effect of errors it is advisable to combine the method with the one described above as follows. We construct the matrix QaP and specify it by means of a number of iterations, obtaining the matrix QI, where 2P < 1. Let us construct AI = Q;AQt and apply to it the process of acceleration etc. A high-speed steady process can be obtained with this combination. It can be proved similarly that the complete eigenvalue problem of any arbitrary matrix can be solved with the help of the triangular power method. Trunslufe~ by
PRASENJIT
BASLJ
REFERENCES 1. FADDEYEV, K. D. and FADDEYEVA, V. N., Vychislitel’nye metody lineinoi algebry. (Computational methods of linear algebra.) Moscow, Fizmatgiz, 1960. 2. FORRYTIiK, G. E. and HENRICI, P., The cyclic Jacobi method for computing the principal values of a complex matrix. Trans. Math. Sot. 94: No. 1, l-23, 1960. 3. GANTMAKHER, F. R., Teorya matrits. (Theory of matrices.) Gostekhizdat, MOSCOW,1954.
IN THIS paper we shall give some iterative methods for the solution of the complete eigenvalue problem of a matrix. These can be easily carried out on electronic computers and the main portion of the work can be done with fixed-point arithmetic without materially reducing the accuracy. The method for solving the complete eigenvalue problem for normal matrices is given in Q 1. It resembles Jacobi’s process in its structure (see [I], 5 Sl). The method of solution in the case of an arbitrary matrix with a simple structure is given in 5 2. It resembles the triangular power method (see [l], 9 78). For simplicity we shall consider real matrices only, although all the methods given can be applied to complex matrices also without particular difficulty. fr.I. METHOD OF SOLVING THE COMPLETE EIGENVALUE A NORMAL
PROBLEM FOR
MATRIX
1. Determination of’ the principal values oj’ a matrix We shall first give a brief account of the method of determining the principal vaIues of a matrix (see [2]), because this will be used in solving the complete problem of eigenvalues of a normal matrix. (Certain changes have been introduced in the description of this method.) The principal values of the matrix A are the positive values of the square root of the eigenvalues of the matrix A*A (or AA*). The method is based on the following result. Let U and V be two orthogonal matrices such that II = UAV, where II is a diagonal matrix. Then the moduli of the diagonal matrix II will give the principal values of the matrix A, the rows of the eigenvectors of the matrix AA*, and the columns of the matrix vectors of the matrix A*A. All this follows from the two obvious ing from (1): II II* = UAA*U*,
I-I*II = V*A*AV
The computing scheme for the method is as follows. 7 Zh. vych. math. 2: No. I, 15-24, 1962.
13
(1) elements the matrix V are the equations
of the U are eigenresult(2)
14
V. V. VOYEVODIN
In the matrix A we choose a pair of non-diagonal elements a, and a,, such that aTj+uji
=
ytlx
and construct the two simple rotation
(&+a&)#
matrices:
cos p
. ..-sin v . ..
u&P =
(3) 1
/1. ..cos y . .
VW = 11
I;
sinv...
i’l
. ..i
i *.
--sinv . .
\
0
cos I .. .
j
. ..I
‘** 1 /
\O
.
The angles 47and y will be defined by the condition that the elements a,j and uji of the matrix Al = U$)AVij) become zero. This condition gives a system of two equations for q and y: COS?/J(Uij
COS lJ?+Ujj
sin y)-siny
(UiiCOS
q+Uji
sin $0)= 0, (4)
COS?#(--U~ising,+UjiCOSpl)+sin~(U~jsinpl+UjjCOS~)=0, or COSQl(UjiCOS
y+Ujj
sin y) -sin
Q1(UiiCOS
sin y) = 0,
y+Uij
(5)
COSpl(-U~iSin?fJ+UijCOSy)+sing,(-Ujisin~+UjjCOS~)=O. These
systems
of. equations In order
algebraic
systems
with respect to cos y, sin y and cos p, sin p with a normalized
can
be considered
solution.
that we may have non-zero
the determinants
of the systems
is easy to obtain formulae
as linear solutions
should
become
for the determination UiiajifUijajj
tan 29 = 2 ___--___.
&-ajj+a:j-aji'
tan
homogeneous it is necessary zero.
From
and sufficient that these
conditions
it
of the angles p and y:
aii aij+ajiajj 2y = 2 2
2.
(6)
aii-ajj+U,Ti-Uij
Those values of cp and y are taken which are determined from (6) and satisfy the system (4) or (5). The matrix A, is constructed after U$) and V$) a similar process being used. The required matrices U and V are determined as infinite products of all the matrices Uiy) and VjT) (m = 1, 2, . . .) respectively, obtained in each step. All theorems on the convergence of this method are proved in exactly the same way as in Jacobi’s method (see [I]).
A method for the solution of the complete eigenvalue problem
1.5
We now pass directly to the solution of the complete eigenvalue problem of normal matrices. The matrix A is called normal if it can be commuted with its conjugate, i.e. AA* = A*A.
(7)
The following proposition was proved for normal matrices (see [3]): in order that a matrix A of’the n-th order may be normal it is necessary and suficient that it should have orthogonal eigenveetors. It follows from this proposition that a normal matrix A can always be represented
in the form A = QAQ*,
(8)
where A is a diagonal and Q is a unitary matrix. Hence A*A =QA*Q*QAQ’
= QA*AQ*;
(9)
therefore the principal values of a normal matrix are equal to the moduli of the eigenvalues of the matrix itself. This and the fact that the principal values can be determined by the method described above make it possible to solve effectively the complete eigenvalue problem for a normal matrix. We consider a few cases of the distribution of eigenvalues. 2. All the eigenvalues of a normal matrix A are simple and di@er in their moduli, except the complex conjugate values
In this case we can obtain an expansion of the matrix A into factors by the method of determination of the principal values (see [I]): A = U’DV’,
(10)
The matrices U and V can always be chosen in such a way that diagonal elements of the matrix D with equal moduli are situated side by side and have identical signs. It follows from (10) that the matrix A will be similar to the matrix D = V’AV = V’U’IIV’V = V’U’II = RD.
(11)
Let us investigate the form of the matrix D. Since the matrix A is normal, the matrix D will also be normal, i.e. the equation DD’ = D’D. will be satisfied. Since the matrix R is orthogonal
(12) we obtain from (12)
RII (RII)’ = (RTI)’ RI-I,
(13)
or RIIII’R
= D’R’RD.
Hence RIP = IPR.
(14)
16
V. V. VOYEVODIN
Let rij be the elements of the matrix R andpFi the diagonal elements of the matrix I12. Then from (13) we obtain that for any pair of indices (k, r) the relation rk; 6
=
dkrklr
WJ)
will be satisfied. Hence if p& # ~7~then rkr = 0. From the assumption regarding the distribution of the eigenvalues of the matrix A it follows that the multiplicity of the diagonal elements of the matrix II2 is not greater than 2. Therefore it is evident that the matrix R will be quasi-diagonal with sub-matrices whose diagonals are not higher than the second order. The matrix D will also have the same form: 0 DI D=
D2.
(\O
(16)
I'hn
Thus
the solution of the complete problem for the matrix A is reduced to the solution of the complete problem of the matrix D which has a very simple form. Let Dj be a sub-matrix of the second order. Then its matrix is equal to the product of an orthogonal and a diagonal matrix:
Since sub-matrices of the second order correspond to complex-conjugate values, therefore ($
$)
eigen-
(17)
= (I;;;-::;;).
Hence we have the final expression for the eigenvalues ;Ifi and eigenvectors zei of the matrix with the sub-matrix Dj: n$jg= pjj(cos~j&isin~j),
z(j) 1.2--
+=; 1.
+J.
Thus with the assumption made regarding the distribution of the eigenvalues of a normal matrix the complete eigenvalue problem is fully solved by using the method of principal values once. If the matrix has more than two eigenvalues with the same modulus, sub-matrices of an order higher than the second may appear in the matrix D. For orthogonal matrices, in particular, this method will not make the solution of the problem any simpler. 3. All the eigenvalues of’ a normal matrix are simple; any number oj them may have the same modulus
Let us represent the matrix A in the form A=B+C,
(19,
A method for the solution of the complete eigenvalue problem where the matrix B is symmetrical and C is skew-symmetrical. tively
the elements
Jacobi’s
aij, bij, cij are respec-
of the matrices A, B, C, then b,
Using
17
=
bji =
method
a-
;
Cij
-cji
=
C
aij-ajj 2 -_
(see [l]) we reduce the symmetrical
matrix B to the
diagonal matrix% and the matrix C to a similar matrix ii,. Then for the matrix A we shall have A = Q(ii+i?)Q’, where Q is an orthogonal
transformation
Since the matrix A is normal the matrix B+6 have the equation (B-j-C)(B+C)’
= (B tZ$‘(BS_C),
W)
matrix.
or
will also be normal. Hence we
(i&j-C)(El-cc),=
(B--C)(B+C),
Hence SC = CB.
(22)
As in the case of eqn. (14) it can be inferred that the matrix B+ecan be reduced to the quasi-diagonal form (16) by transposing the rows and columns, only the sub-matrices Dj(jj = 1, 2, . . . , m) in this case will have an arbitrary order lj, but these must be of the form: Dj 1~fDj-tdjE,
(23)
where aj is a skew-symmetrica matrix, and dj is a real number. The eigenvaiues A$$ (k z- 1,2 , . . . , lj) of the matrix Dj are retated to the eigenvalues @ of the matrix 6j by the equation @‘=
$‘+&
J’
(24)
and their eigenvectors coincide so that further solution of the complete problem for normat matrices is reduced to a solution of the complete problem for skewsymmetrical matrices By hypothesis the matrix A has no multiple eigenvalues; the skew-symmetrical 5j.
matrix fij also will not have any. Hence the method described in section 2 can be successfully used for the solution of the complete eigenvalue problem of the latter niatrix. A special case is of particular interest. Let the matrix A be orthogonal and have no multiple eigenvalues. In that case the sub-matrices Dj(j = 1, 2, . . . . m) will also be orthogonal matrices. It is known (see [3]) that all the eigenvalues of an orthogonal matrix have a modulus equa1 to unity and all eigenvalues of a skew-symmetrical matrix are pure imaginary numbers. Hence from equation (24) it follows J$‘!‘+$=
I,
k=
1,2, . . ..ij.
G5)
will be satisfied. Since there can be no more than two different pure imaginary numbers @ satisfying the equation
(29, it follows that none of the sub-matrices
v. v.
18
VoYEvoDIN
ofD,(j = 1,2, . . . m) will be of an order higher than the second. Thus if an orthogonal matrix has no multiple eigenvalues the solution of its complete eigenvalue problem can be obtained by using Jacobi’s method once. Finally it should be noted that when this method is used in practice the expansion (19) is unnecessary, and the orthogonal transformations should be applied with it to the matrix A after calculating the angle of rotation from the formula afj+aji tan29 = ~ #ii-ajj
*
These transformations are equivalent to reducing the matrix B by Jacobi’s method to the diagonal form, but the matrix A can be reduced directly to the matrix B+C. To ensure maximum speed of convergence at each step the pair of indices (i,j) should be so chosen that the sum Gfj+aji has the m~imum modulus. 4. The normal matrix A has multiple eigenvalues
It can be easily seen that if the method described in section 3 is used in this case, we shall ultimately have to solve the complete problem for matrices N which are both skew-symmetrical and orthogonal, i.e. for which N’=
N’N=NN’=E.
-N,
(27)
from (27) we easily get NZ= -E.
(28)
The eigenvalues of the matrix N are equal to ki, and the eigenvectors can be easily found by remembering the following. Let p be an arbitrary non-zero vector, then the vector p-iNp will be an eigenvector of the matrix N corresponding to an eigenvalue equal to +i. Bearing (28) in mind we have N(p-iNp)
= Np--iNap
0 2. METHOD OF SOtVING
=Np+ip
= i(p-iip).
THE COMPLETE EIGENVALUE ARBITRARY MATRIX
(29)
PROBLEM OF AN
Let a real non-specialized matrix A be given, corresponding to an operator which is real in n-dimensional space R and has n linear independent eigenvectors Xl. . * . . x,. Let us divide the eigenvalues Ai (i = 1,2, . . . , n) of the matrix A into m groups, each group containing only eigenvalues with the same modulus, and suppose that the condition (30) is satisfied, where k, = n. We construct m sub-spaces PX, Pz, . . . . P,,, of the space R invariant with respect to A by defining Pj G = 1,2, . . . , m) as a sub-space formed by the eigenvectors x1, . ..*Q. Evidently the following relations hold (31)
A method for the solution of the complete eigenvalue problem
19
with the vectors ul, uZ, . . . , u,, form a basis of the space R having the property that the VeCtOrS II,, . . . , Ukj (j = 1, 2, . . . , m), form a basis of the sub-space Pj (j = 1, 2, .** , m). Let us determine the form of the matrix of an operator given on this base, (see [l], 5 7). Since the vectors AU,, . . . AUkI (j = 1, 2, . . . , m) belong to Pj i.e. are linear combinations of the vectors ul, . . . , ut, only, their coordinates on the base selected beginning from kj+ 1 are equal to zero. Consequently the matrix A on the new base has a quasi-triangular form AllA,, . . . A,, A=
0 A22...A2m . . . . . +. . . . . . . . . (0
(32)
0 e.. A,,,,,,1.
where the sub-matrix Ajj(j = 1,2, . . . , m) has the order kj-kj_1, k. = 0. If Q is a matrix the columns of which are respectively the coordinates of the vectors ul, . . . u, given on the old base, the matrices A and A are connected by the relation A===Q-‘AQ.
(33)
The problem of finding the eigenvalues of the matrix A is considerably simpler than for A, because it amounts to finding the eigenvalues of the matrices Ajj (j = 1,2, . . . ) m), which are usually of a low order. In particular, if the eigenvalues of the matrix A all have different moduli, the sub-matrices Ajj(j = 1,2, . . . , m) will be of the first order. Thus the solution of the complete eigenvalue problem of the matrix A or its simplification when there are eigenvalues with muftiple moduli is reduced to finding a matrix Q, which satisfies the eqn. (33) where A is the matrix in (32). The matrix Q is found as follows. An arbitrary orthogonal matrix Q. is taken and the sequence of orthogonal matrices Qk(k = 1, 2, . . . ,) is constructed by a method similar to the triangular power method;
AQo=Q,A,, AQ,=Qa&, . . . . . . . . *..a . . . . . . . . .
(34)
AQc=Qk+~h+~r Here A, (k = 1,2, .a.) are right triangular matrices. It is best to expand the matrix into a product of orthogonal and right triangular matrices by the reflection method (see PI, § 16). We shall prove that the sequence of matrices Q,(k = 1,2, . , .) has the same property as the sequence of matrices Ak = Q;AQk tend to a quasi-triangular matrix of the form (32). To do this we have to evaluate the minors of the matrix & = Q;AtkAkQO.
(35)
Let A = QOCAC-“Q;,
R=[l,,...,J,]
(36)
V. V.
20
VOYEVODIN
then B, = C-l’AkC’CAkC-l
(37)
are introduced for the minors of the matrix B,:
the following abbreviations
bfk! btk! se,1 ‘a,2
b:,k,!, ......... ... ......
b!k! bfk!
llJa.*.
1111
+.-
(38)
btk> ‘IJI
Using Binet Cauchy’s formula (see [3]) for the determinant of the product of two rectangular matrices, we obtain Bk(;;;;:
::: j
= I_~~~~~,cn,_~C-l~AkC’ xCnkC-l
=
c
I
x
C-1’
c
(;;,;:,“.‘.:f’,,)
(3g)
x
Pl,Pa, ...>Pr (h,j2, .. ..A ) (;I:
:/Jfq;:;
:::; ;;)“‘(“b,;
:::;,“:)
x
l~BIC...~B~~”
,_.c,, . .
’
in q-I:: 1’
::::Epk(:::
:p’(:::::::;)*
::::
since
we shall have, taking into account (30), (40) where O$j are quantities bounded above with respect to the modulus for all i, j, k, I by a constant independent of k. The expression for principal minors of the matrix B, is Bk(;:;::::::) =
x l
[
c
c-l(;:~.:::;“‘)c~::
:::J;)A:‘“:“...XL
]
l~al~...G+” =
pppp
. . . pp@p.
(41)
The matrix Q,, can always be so chosen that all the quantities Oj”) will be bounded below for all k by a strictly positive number independent of k. From (34) we obtain AkQ, = Qk’bk’&-I...&
= Qk.Vk,
(42)
whence B=V;Vk.
(43)
A method
for the solution
of the complete
eigenvalue
21
problem
Let 8$) be the elements of the matrix Vk. Then we know (see [3]) that
#),.j(&) iJ
j=i+l,i+2
II
,...,
n.
It can be easily deduced from (40) (41) (44) that the matrix bk can be represented in the form Vk= Dk.V/‘,
(45)
where D, = [I1llk,...,lil,[k] and V, is a right triangular matrix with elements bounded above with respect to the modulus for all values of k by a constant independent of k. The diagonal elements of the matrix Vk are bounded below for all values of k by a strictly positive number independent of k. Hence vi’ is a right triangular matrix with elements bounded above with respect to the modulus for all values of k by a constant independent of k. From (42) we obtain Qk = AkQOVkl.
(46)
If we form the matrix A, and take into account (45) and (46), we shall have Ak = Q;AQk = QklAQk = VkQolA-kAAkQ,,Vil = VkQ;‘AQ,,V;l
= DkVkQ,lAQ,V;lD+=
D,N,D;l.
(47)
The elements of the matrix Nk are numbers bounded above with respect to the modulus for all values of k by a constant independent of k. Let us divide up the matrix A, into submatrices Ai;) as in (32) and let Zj(j = 1, 2, . . ., m) be a number equal to the modulus of the eigenvalues in the jth group. Then it can be deduced from (47) that the elements of the submatrix A$) for i > j are quantities
Hence the matrices A, tend to a quasi-triangular form. The chief merit of this method is that it is free from errors of rounding off, because orthogonal transformations of matrices are used and this makes it selfcorrecting. But when eigenvalues with nearly equal moduli are present this method gives a very weak convergence. This difficulty can be overcome by using the following analogue of the M-algorithm (see [I], 0. 80). Let AzpQ P V P’
(48)
22
V.V.
where Q, is an orthogonal
VOYEVODlh'
and e,, a right triangular
matrix. With the notations
introduced above we shall evidently have QP = QLp, and fi, = WpP when Q. = E. Squaring equation (48) we obtain A @+I = ijpijptjpi?,.
(4%
If the matrix e7,Q, is expanded as a product of an orthogonal a right triangular matrix Lp, then A ap+1 = i&R,L,i$
matrix R, and
= ~,,,vp,,,
(50)
where ii p+l = QPRP,
and eP+r
= Pep.
(51)
The matrices Q1, Qz, Qa, QB, . . . . Q2p can be determined successively by this method. The method will converge under the conditions of the principal method, the convergence being of the order of
i.e. quadratic. It can be inferred from (45) that, as a rule, the elements of the matrix GP will either increase or decrease rapidly with p. To reduce the effect of errors it is advisable to combine the method with the one described above as follows. We construct the matrix QaP and specify it by means of a number of iterations, obtaining the matrix QI, where 2P < 1. Let us construct AI = Q;AQt and apply to it the process of acceleration etc. A high-speed steady process can be obtained with this combination. It can be proved similarly that the complete eigenvalue problem of any arbitrary matrix can be solved with the help of the triangular power method. Trunslufe~ by
PRASENJIT
BASLJ
REFERENCES 1. FADDEYEV, K. D. and FADDEYEVA, V. N., Vychislitel’nye metody lineinoi algebry. (Computational methods of linear algebra.) Moscow, Fizmatgiz, 1960. 2. FORRYTIiK, G. E. and HENRICI, P., The cyclic Jacobi method for computing the principal values of a complex matrix. Trans. Math. Sot. 94: No. 1, l-23, 1960. 3. GANTMAKHER, F. R., Teorya matrits. (Theory of matrices.) Gostekhizdat, MOSCOW,1954.