Perturbation of invariant subspaces^{∗}

Ninoslav Truhar^{†}

Abstract. We consider two different theoretical approaches for the problem of the perturbation of invariant subspaces. The first approach belongs to the standard theory. In that approach the bounds for the norm of the perturbation of the projector are proportional to the norm of perturbation matrix, and inversely proportional to the distance between the corresponding eigenvalues and the rest of the spectrum. The second approach belongs to the relative theory which deals only with Hermitian matrices. The bounds which result from this approach are proportional to the size of relative perturbation of matrix elements and the condition number of a scaled matrix, and inversely proportional to the relative gap between the corresponding eigenvalue and the rest of the spectrum. Because of relative gap these bounds are in some cases less pessimistic than the standard norm estimates.

Key words: perturbation bound, invariant subspace, orthogonal projection

Sažetak. Prikazali smo dva teorijska pristupa koji se bave problemom perturbacijama invarijantnih potprostora. Veličina perturbacije mjeri se normom pripadnih projektora, jer projektori ne ovise o izboru baza promatranih potprostora. Prvi pristup pripada standardnoj perturbacijskoj teoriji. U tom pristupu su ocjene za normu perturbacije projektora proporcionalne normi matrice perturbacije, a obrnuto proporcionalne udaljenosti odgovarajuće vlastite vrijednosti od ostatka spektra. Drugi pristup pripada tzv. relativnoj perturbacijskoj teoriji, koja promatra samo perturbacije invarijantnih potprostora hermitskih matrica. Ocjene koje slijede iz tog pristupa su proporcionalne veličini relativne promjene matričnih elemenata i uvjetovanosti skalirane matrice, a obrnuto proporcionalne s relativnom udaljenošću odgovarajuće vlastite vrijednosti od ostatka spektra. Zbog relativne udaljenosti u nekim situacijama su relativne perturbacijske ocjene manje pesmistične od ocjena standardne perturbacijske teorije.

Ključne riječi: perturbacijska granica, invarijantni potprostor, ortogonalni projektor

∗Predavanje održano u okviru ZNANSTVENOG KOLOKVIJA, HMD – Podružnica Osijek, 8.

prosinca 1995.

†Faculty of Civil Engineering, Department of Mathematics, Drinska 16 a, HR-31 000 Osijek, e-mail: truhar@oliver.unios.hr

One of the important topics in matrix theory is an invariant subspace of the matrix M. We have the following definition:

the subspace X is an invariant subspace of the matrix M if

M X ⊂ X .

Let the columns of the matrix X form the basis for the invariant subspace X of the matrix M. Then there is a unique matrix L such that (see [11, p. 22]).

M X = XL.

The matrix L is the representation of the matrix M on the subspace X with respect to the basis X.

Let us denote the spectrum of M by α(M) = { λ_{1}, λ_{2}, ..., λ_{l} }. We can easily show that α(L) ⊂ α(M), that is, if L is the representation of M on X then the eigenvalues of L are also the eigenvalues of M. Therefore, we shall often say that the invariant subspace X = R(X) corresponds to the eigenvalues from α(L), where R(X) denotes the column space of X (i.e. R(X) = {X x : x ∈ C^{n}}).

For a Hermitian matrix H we can find a spectral decomposition H = U Λ U^{∗}, where U is a unitary matrix, and Λ is a real diagonal matrix, whose diagonal elements are the eigenvalues of H. In other words, the Hermitian matrix of order n has the system of orthonormal eigenvectors that span C^{n}. Thus, every invariant subspace of the Hermitian matrix H has an orthonormal basis of eigenvectors of H. This means that the invariant subspace X of H which corresponds to the eigenvalues from the set α(L) will be spanned by the eigenvectors which correspond to the eigenvalues from α(L).

In many applications where we need to determine eigenvalues and the corresponding eigenvectors of a given matrix, and one of the most important criteria which must be fulfilled is the stability with respect refers to the perturbation of matrix elements. By stability we colloquially mean that small changes in matrix elements cause changes in the eigenvalues and invariant subspaces of similar order.

In order to be able to measure perturbations of invariant subspaces, we need to define the distance between subspaces. Let the columns of orthonormal matrices X and Y form the orthonormal basis for invariant subspaces X = R(X) i Y = R(Y), respectively. The matrices

P_{X} = X X^{∗} P_{Y} = Y Y^{∗}

are called the orthogonal projections onto X and Y, respectively. The distance between two subspaces, is defined as norm of difference between the corresponding orthogonal projections ([4]):

dist(X , Y) =^{def} ∥P_{X} − P_{Y}∥.

Throughout this paper ∥ · ∥ denotes the spectral norm. Moreover, we have

∥P_{X} − P_{Y}∥ = sin θ_{1},

where θ_{1} is the greatest canonical angle between X and Y, which is defined by [11]:

cos θ_{1} = min_{{x∈X \atop x≠0}} max_{{y∈Y \atop y≠0}} \frac{y^{∗}x}{∥x∥ ∥y∥}.

We can now define our problem: let X be the invariant subspace of the matrix M which corresponds to eigenvalues from α(L) = { λ_{1}, ..., λ_{i} } ⊂ α(M), and let P be the orthogonal projection onto X . Let X̃ be the invariant subspace of matrix M + δM, which corresponds to eigenvalues from α(L̃) = { λ̃_{1}, ..., λ̃_{i} } ⊂ α(M̃), and let P̃ = P + δP be the orthogonal projection on to X̃. Here δM is a perturbation, and λ_{i} and λ̃_{i} are given in the same order. We want to determine the distance between subspaces X and X̃, dist(X , X̃) = ∥P̃ − P ∥ = ∥δP ∥.

There are two different approaches to the problem of the perturbation of invariant subspaces. The first approach gives the perturbation bound in terms of ∥ δM ∥. This approach is based upon the standard perturbation theory, and it can be applied to all quadratic matrices. In general, the perturbation bounds which follow from this approach have the form [11, 4]

∥ δP ∥ ≤ \frac{∥ δM ∥}{min ; \begin{matrix}λ ∈ α(L) \\ μ ∈ α(M)\α(L) \end{matrix} ; |λ − μ|},

The second approach is based upon the relative perturbation theory and it considers the relative perturbation of matrix elements of Hermitian matrices. For example, with this approach we can consider perturbations δH such that

|δH_{ij}| ≤ ε |H_{ij}|.  (1)

Such perturbations occur, for example, whenever the matrix H is stored into the computer memory, that is, whenever H is given with a certain number of correct digits. The bounds which follow from this approach have the following form [3, 14, 12]:

∥ δP ∥ ≤ \frac{ε ∥|A|∥ ∥Â^{−1}∥}{relgap(λ)}.  (2)

Here relgap(λ) is the relative distance between λ from the rest of the spectrum of H, and the matrices A and Â are defined as follows: let H = QΛQ^{∗} be the eigenvalue decomposition of H. The spectral absolute value H is defined by

H = Q |Λ| Q^{∗} = \sqrt{H^{2}}.  (3)

Further,

Â = D H D,  A = D H D,  (4)

where D is a diagonal positive definite matrix such that Â_{ii} = 1, that is, D =

( diag( H ) )^{1/2} . Also, the matrix |A| is defined by (|A|)_{ij} = |a_{ij}|. Finally, note that H needs to be non-singular.

Let us explain the basic characteristics of each approach: the first approach obviously produces better bounds if |λ| is close to ∥M∥, and λ is well separated from the rest of the spectrum of M. The bound (2) will e.g. also be good for tiny λ if all other eigenvalues are of O(∥M∥) in absolute value. The second approach gives an almo st uniform bound for all kinds of subspaces, independently of the magnitude of eigenvalues. But for the matrix H which has a condition number much greater than a condition number of the scaled matrix Â (i.e. κ(H) >> κ(Â)), and the relative distance is not very small, the bound can be much better than bound (1).

Within the framework of the standard perturbation theory we shall give two different types of results, by Stewart [10] and Davis and Kahan [2].

Bounds which follow from the results of Stewart [10] are called the approximation bounds since they follow from the approximation theorem. These bounds hold for all kinds of quadratic matrices. We can describe the basic approximation result in the following way: let [ X Y ] be a unitary matrix, and let

[ X Y ]^{∗} M [ X Y ] = [ \begin{matrix}L_{1} & F \\ G & L_{2} \end{matrix} ] ,  (5)

where X ∈ C^{n×l}, and M is a quadratic matrix. The space R(X) is an invariant subspace for matrix M if and only if G ≡ Y^{∗}MX = 0. If G is sufficiently small there will be an invariant subspace R(X̃) of matrix M, which approaches R(X) as G approaches zero. Stewart [10] shows that there exists an invariant subspace R(X̃) of the matrix M for which

∥ δP ∥ ≤ 2 \frac{∥G∥}{sep(L_{1}, L_{2})},  (6)

where P and P̃ = P + P̃ are orthogonal projections onto R(X) and R(X̃), respectively, and sep(L_{1}, L_{2}) is a function defined by

sep(L_{1}, L_{2}) =^{def} inf_{∥X∥=1} ∥ L_{1}X − XL_{2} ∥.

Obviously, the bound (6) is meaningful whenever sep(L_{1}, L_{2}) > 0. This condition is fulfilled whenever T = L_{1}X − XL_{2} is a non-singular operator. In [11] it was shown that T is non-singular whenever the spectrums of L_{1} and L_{2} are disjoint, α(L_{1}) ∩ α(L_{2}) = ∅. Because of this, we shall consider only perturbation of invariant subspaces whose representations have a disjoint spectra. The subspaces of this kind are called simple invariant subspaces.

Now we consider a unitary matrix [ X Y ] , and assume that the columns of X span a simple invariant subspace X of M. Stewart and Sun [11] showed that there exists a matrix X̃ whose columns form the orthonormal basis for an invariant subspace X̃ of M̃ = M + δM, for which

dist(X , X̃) = ∥ δP ∥ ≤ 2 \frac{∥ M̃ X − X (X^{∗}M̃ X) ∥}{sep(L_{1}, L_{2}) − ∥X^{∗}δM X∥ − ∥Y^{∗}δM Y ∥},  (7)

Here L_{1} and L_{2} are defined by (5). The numerator in the bound (7) is the norm of the residual R = M̃ X −X (X^{∗}M̃ X). This norm ∥R∥ is the smallest norm referring to all residuals of the form M X − X L (see [11, p. 176]). The perturbation bound in the Hermitian case follows directly from the bound (7).

From a different point of view, Davis and Kahan [2] develop the so called direct bounds for the perturbations of Hermitian matrices. Their approach uses the fact that for Hermitian matrices the existence of the perturbed invariant subspace, can be assumed under mild conditions. We state one of the ”sin θ” theorems which gives the residual bound for sines of the greatest canonical angle between an invariant subspace of the Hermitian matrix H and its perturbation H + δH. Briefly let [ X Y ] and [ X̃ Ỹ ] be unitary and

[ X Y ]^{∗} H [ X Y ] = [ \begin{matrix}L_{1} & 0 \\ 0 & L_{2} \end{matrix} ] ,  (8)

and

[ X̃ Ỹ ]^{∗} (H + δH) [ X̃ Ỹ ] = [ \begin{matrix}L̃_{1} & 0 \\ 0 & L̃_{2} \end{matrix} ] .

Let

R = (H + δH) X − X L_{1}

be the residual, and let the spectrums α(L_{1}) and α(L̃_{2}) be separate as in figure:

α(L̃_{2})  α(L_{1})  α(L̃_{2})  .

If we set δ = min |α(L_{1}) − α(L̃_{2})| ≡ min{|λ_{i} − λ_{k}| : λ_{i} ∈ α(L_{1}), λ_{k} ∈ L̃_{2}}, than we have

∥ δP ∥ ≤ \frac{∥R∥}{δ}.  (9)

For Hermitian matrices with

sep(L_{1}, L_{2}) − ∥X^{∗}δH X∥ − ∥Y^{∗}δH Y ∥ ≈ δ,

the perturbation bound (7) from the approximation theorem can be much better than the direct bound (9) because the norm of the residual in (7) can be much smaller than the norm of the residual in (9). But for

δ > sep(L_{1}, L_{2}) − ∥X^{∗}δH X∥ − ∥Y^{∗}δH Y ∥,

and when the norms of both residuals are close, than the direct bound (9) is at least two times better than approximation bound (7).

Notice, that for both bounds, (9) and (7), ∥δP ∥ is bounded by a function of the norm of the perturbation. But these bounds sometimes can be useless: to illustrate this, consider the symmetric positive define matrix H = D A D from [3], where D = diag(10^{20}, 10^{10}, 1),

H = [ \begin{matrix}10^{40} & 10^{29} & 10^{19} \\ 10^{29} & 10^{20} & 10^{9} \\ 10^{19} & 10^{9} & 1 \end{matrix} ] , A = [ \begin{matrix}1 & 0.1 & 0.1 \\ 0.1 & 1 & 0.1 \\ 0.1 & 0.1 & 1 \end{matrix} ] .

To six correct figures, H′s eigenvalue matrix Λ and eigenvector matrix V (normalized to have the largest entry each eigenvector equal to 1) are

Λ = diag(1.00000 · 10^{40}, 9.90000 · 10^{19}, 9.81818 · 10^{−1})

and

V = [ \begin{matrix}1.00000 & −1.00000 · 10^{−11} & −9.09091 · 10^{−22} \\ 1.00000 · 10^{−11} & 1.00000 & −9.09091 · 10^{−12} \\ 1.00000 · 10^{−21} & 9.09091 · 10^{−12} & 1.00000 \end{matrix} ] .

Here κ(H) ≈ 10^{40} and κ(A) ≈ 1.33. Let δH be Hermitian perturbation

δH = [ \begin{matrix}2.2 · 10^{33} & 3.6 · 10^{22} & 5.9 · 10^{12} \\ 3.6 · 10^{22} & 9.3 · 10^{13} & 6.1 · 10^{2} \\ 5.9 · 10^{12} & 6.1 · 10^{2} & 3.5 · 10^{−8} \end{matrix} ] .

where |δH_{ij}| ≤ ε|H_{ij}|, ε = 10^{−6}. For perturbation of invariant subspace spanned by second eigenvector (i.e. the second column of the matrix V), we can use neither. The bound (9) since

sep(L_{1}, L_{2}) − ∥X^{∗}δH X∥ − ∥Y^{∗}δH Y ∥ = −1.2 · 10^{−33},

nor the bound (7) which gives

∥ δP ∥ ≤ \frac{1.2998 · 10^{23}}{9.9 · 10^{19}} ≤ 1.32 · 10^{3}.

Because of this, it was necessary to change the approach to this problem. Let us give some examples of the relative approach. First, Barlow and Demmel [1] considered the scaled diagonally dominant matrices, that is matrices of the form

H = D A D, A = E + M,

where D is diagonal and non-singular, E is diagonal with elements ±1, diag(M) = 0 and ∥M∥ = ζ < 1. They showed that for such matrices relative changes of matrix elements (1) imply

∥v′_{i} − v_{i}∥ ≤ \frac{(n − 1) ε}{(1 − ζ) relgap(λ_{i})} + O(ε^{2}),  (10)

where v_{i} is the eigenvector corresponding to the simple eigenvalue λ_{i}, v′_{i} is the eigenvector corresponding to the simple eigenvalue λ′_{i}. The function relgap(λ_{i}) in (10) is the relative gap between the eigenvalue λ_{i} and the rest of the spectrum α(H),

relgap(λ_{i}) = min_{j≠i} \frac{|λ_{i} − λ_{j}|}{\sqrt{|λ_{i} λ_{j}|}}.  (11)

Further, Demmel and Veselić [3] considered perturbation as in (1) for positive definite Hermitian matrices. They proved the following result: write H = D A D, where D = ( diag(H) )^{1/2} is a scaling so that A_{ii} = 1. Than:

∥v′_{i} − v_{i}∥ ≤ \frac{(n − 1)^{1/2} κ(A)ε}{relgap(λ_{i})} + O(ε^{2}),  (12)

where κ(A) = ∥A∥ · ∥A^{−1}∥ is the condition number of A, and relgap(λ_{i}) is defined by (11). Here it is also assume that λ_{i} and λ′_{i} are simple.

By the theorem of Van der Sluis [9]

κ(A) ≤ n · min_{D} κ(D H D) ≤ n · κ(H),

that is κ(A) nearly minimizes the condition number of the positive definite H over all possible diagonal scalings. Clearly, it is possible that κ(A) << κ(H), so the bound (12) is always at least about as good and can be much better than the bounds from the standard perturbation theory, whenever the denominator in (7) or in (9) is less than or equal to the relative gap (11).

Veselić and Slapničar [14] generalized the above results to non-singular Hermitian matrices. They proved the following result: let H be defined by (3) and let the perturbation δH satisfy

|x^{∗}δH x| ≤ η x^{∗} H x,  for all x ∈ C^{n}.  (13)

Than the perturbation bound for the eigenprojection corresponding to possibly multiple eigenvalue λ, is given by

∥δP ∥ ≤ \{ \begin{matrix}\frac{η}{rg(λ)} \frac{1}{1 − (1 + \frac{1}{rg(λ)}) η} & for λ_{L} > 0, 2\sqrt{λ} − \sqrt{λ_{L}} < \sqrt{λ_{R}}, \\ \frac{η}{rg(λ)} \frac{1}{1\frac{η}{rg(λ)}} & u otherwise, \end{matrix} (14)

provided that the right side is positive. The function of the relative gap rg(λ) for the positive eigenvalue λ is defined by

rg(λ) = \{ \begin{matrix}min \{\frac{\sqrt{λ} − \sqrt{λ_{L}}}{\sqrt{λ_{L}}} , \frac{\sqrt{λ_{R}} − \sqrt{λ}}{\sqrt{λ_{R}}}\} , & λ_{L} > 0, \\ min \{2(\sqrt{2} − 1), \frac{λ_{R} − λ}{λ_{R} + λ}\} , & otherwise. \end{matrix}  (15)

In (14) and (15) λ_{R} and λ_{L} denote the left and right neighbour of λ in the spectrum α(H) of H, respectively. Negative eigenvalues are considered as the positive eigenvalues of −H.

The perturbation of the type (13) occurs, for example whenever H is given with floating-point error as in (1). Namely, from (1) it follows

|x^{∗}δH x| ≤ |x|^{T} ε|H| |x| ≤ ε|x|^{T} D |A| D x ≤ ∥|A|∥·ε x^{∗}D^{2} x ≤ ε∥|A|∥·∥Â^{−1}∥x^{∗} H x,

where A and Â are defined by (3) and (4). We see that the perturbation of the form (1) implies (13), with η = ε ∥|A|∥ · ∥Â^{−1}∥.

This result generalizes the above results of Barlow and Demmel and Demmel and Veselić. First, in [14] it is shown that for Hermitian matrices of the form

H = D (E + N) D,

where D is diagonal and positive definite, E = E^{∗} = E^{−1} and ∥N∥ = 1, the perturbation of the type (1) implies (13) with

η = ε n \frac{1 + ∥ |N| ∥}{1 − ∥N∥}.

This, together with the fact that for a simple eigenvalue λ and its eigenvector v we have

∥δv∥ ≤ \sqrt{2} ∥δP ∥,  (16)

generalizes (10). Further, if H is positive definite, then H = H, and (14) together with (16) is similar to (12) with slightly different relative gap.

Slapničar and Veselić [14, 8] also proved the relative perturbation bound for the Hermitian matrix H given in the factorized form

H = G J G^{∗},

where G is a n×r matrix of the full column rank, and J is a non-singular Hermitian matrix, under the perturbation of the factor G. More precisely, the perturbed matrix H̃ is defined by

H̃ = (G + δG) J (G + δG)^{∗},

where

∥δG x∥ ≤ η∥G x∥,

for all x ∈ C^{r} and some η < 1. The most common J is of the form

J = [ \begin{matrix}I_{m} & 0 \\ 0 & −I_{r−m} \end{matrix} ]

in which case m, r − m, and n − r is the number of the positive, negative and zero eigenvalues of H, respectively. The perturbation bound for the eigenprojection reads

∥δP ∥ ≤ \frac{4 η̃}{rg_{G}(λ)} · \frac{1}{1 − \frac{3 η̃}{rg_{G}(λ)}} ,  (17)

where η̃ = η (2 + η), provided that the right hand side in (17) is positive. The relative gap rg_{G}(λ) is for the positive eigenvalue λ defined by

rg_{G}(λ) = min \{\frac{λ − λ_{L}}{λ + λ_{L}} , \frac{λ_{R} − λ}{λ_{R} + λ}, 1\} .  (18)

Here λ_{L} and λ_{R} denote the left and right neighbour of λ in the spectrum α(H) of H, respectively.

In practice the best bounds for invariant subspaces are obtained by combining both approaches, that is by combining (7), (9) (14), (17). Some interesting examples are given in [12] and [13].

Singer [6], improved slightly the bound (14). Ren-Cang Li [5] generalized the result of Demmel and Veselić by establishing the perturbation bound for the simple invariant subspace which corresponds to the set of the leftmost or the rightmost eigenvalues of the positive definite Hermitian matrix. Lastly, Truhar [12] and Truhar and Slapničar [13] generalized the results from [1, 3, 14, 8, 5] by giving the perturbation bounds for the simple invariant subspaces which correspond to a set of neighboring eigenvalues of an nonsingular (possibly indefinite) Hermitian matrix.

References

[1] J. Barlow, J. Demmel, Computing Accurate Eigensystems of Scaled Diagonally Dominant Matrices, SIAM Journal on Numerical Analysis Vol. 27, No. 3, (762-791) 1990.

[2] C. Davis, W. M. Kahan, The Rotation of Eigenvectors by a Perturbation. III, SIAM Journal on Numerical Analysis Vol. 7, No. 1. (1-46) 1970.

[3] J. Demmel, K. Veselić, Jacobi’s Method is More Accurate than QR, SIAM Journal on Matrix Analysis and Applications, Vol. 13. No. 4, (1204-1244) 1992.

[4] G. H. Golub, C. F. Van Loan, Matrix Computation, Johns Hopkins, Baltimore, 1989.

[5] Ren-Cang Li Relative Perturbation Theory: (ii) eigenspace variations. Technical report, Department of Mathematics, University of California at Berkeley, 1994.

[6] S. Singer, Računanje spektralne dekompozicije pozitivno definitne matrice, M. S. thesis, University of Zagreb, Zagreb, 1993.

[7] I. Slapničar, Accurate Symmetric Eigenreduction by Jacobi Method, Ph. D. thesis, Fernuniversität, Hagen, 1992.

[8] I. Slapničar, K.Veselić, Perturbations of the Eigenprojections of a Factorised Hermitian Matrix, Linear Algebra Appl., 218 (273-280) 1995.

[9] A. van der Sluis, Condition Numbers and Equilibration of Matrices, Numerische Matematik, 14. (14-23) 1969.

[10] G. W. Stewart, Error and Perturbation Bounds for Subspaces Associated with Certain Eigenvalue Problem, SIAM Review Vol 15. No. 4. (727-764) 1973.

[11] G. W. Stewart, Ji-Guang Sun, Matrix Perturbation Theory, Academic Press, 1990.

[12] N. Truhar, Perturbacije invarijantnih potprostora, M. S. thesis, University of Zagreb, Zagreb, 1995.

[13] N. Truhar, I. Slapničar, Relative Perturbation of Invariant Subspaces, Technical report, University of Osijek and University of Split, 1995, in preparation.

[14] K .Veselić, I. Slapničar, Floating - point Perturbation of Hermitian Matrices, Linear Algebra Appl., 195 (81-116) 1993.