Rigorous Asymptotic Perturbation Bounds for Hermitian Matrix Eigendecompositions

1. Introduction

The eigenvalue perturbation analysis of real symmetric and complex Hermitian matrices is a well established part of matrix analysis which is presented in depth in the fundamental works of Kato [12], Wilkinson [31], Parlett [20], Stewart and Sun [27], Bhatia [4], Stewart [26] and Chatelin [7], the surveys of Sun [29] and Li [14], as well as in numerous papers, see for instance [2,11,16,25,28,30] and the references therein. The perturbation analysis of Hermitian decompositions is usually simpler than the analysis of non-symmetric matrix decompositions. (For a survey on the perturbation theory of non-symmetric matrix decompositions, see [10]). The aim of such an analysis is to find perturbation bounds of the eigenvalues and eigenvectors of Hermitian matrices in the case when these matrices are subject to Hermitian or non-Hermitian perturbations. Depending on the problem solved, such an analysis can be a priori, when we want to predict in advance the changes in the eigendecomposition before the perturbation is applied, or a posteriori, when we attempt to find the errors in eigenvalues and eigenvectors based on the computed perturbed quantities [7]. Here we will be concerned with the first type of analysis, the second one being related to the determination of residual bounds as it is considered by several authors, see Davis and Kahan [8], Parlett [Ch. 11][20], Stewart and Sun [Ch. 5][27], Chatelin [7], and Nakatsukasa [18]. On the other hand, the perturbation analysis can be asymptotic, when we are interested in the effect of vanishing perturbations or global, when we want to find guaranteed bounds in the case of sufficiently large perturbations. We note that in most cases the authors are usually concerned with the eigenvalue perturbations and the eigenvector perturbation analysis is reduced to the sensitivity analysis of the corresponding invariant subspace. In several cases this sensitivity analysis is restricted to bounding the distance between the unperturbed and perturbed subspace [28] or finding the angles between these subspaces [14,26,27,29]. However, in some applications it is necessary to have bounds on the individual elements of specific eigenvector, not only on the sensitivity of the corresponding invariant subspace. This leads to the necessity to perform a componentwise perturbation analysis of the matrix eigensystem.

In this paper we are interested to carry out a componentwise perturbation analysis of indefinite Hermitian matrices that allows to bound each element of every corresponding eigenvector. This analysis is done only for regular eigenvalue problems, i.e., perturbation problems for matrices with distinct eigenvalues. The singular problems, corresponding to multiple eigenvalues, are treated by different techniques, see [6,15]. The regular problems can be solved simply and efficiently by using the method of splitting operators [13] which is already implemented to solve several other matrix perturbation problems. It is important that the bounds on eigenvector elements, obtained by this method, can also be used to obtain a bound on the sensitivity of the invariant subspace spanned by the corresponding eigenvector.

Theoretically, the perturbation analysis of Hermitian decompositions can be done by using some of the methods intended for non-Hermitian problems. This, however, may require much larger memory and volume of computations. For instance, the similar method for perturbation analysis of the Schur form proposed in [21], requires to construct large matrices by implementing the Kronecker product. The method presented in this paper avoids the using of such product which makes possible to analyze much larger problems.

Consider a Hermitian matrix $A\in {\mathbb{C}}^{n\times n}$. Using a unitary transformation matrix $U\in {\mathbb{C}}^{n\times n}$, the matrix A is reduced to the diagonal form

$$\Lambda =U^{H}AU=\mathrm{diag}({\lambda }_1,{\lambda }_2,\dots ,{\lambda }_n),$$

(1)

where ${\lambda }_1,{\lambda }_2,\dots ,{\lambda }_n$ are the eigenvalues of A and $U=[U_1,U_2,\dots ,U_n]$ is the matrix of the corresponding orthonormal eigenvectors. Note that the eigenvalues of a Hermitian matrix are always real. If A is real, then U is orthogonal [9,20]. The diagonal decomposition (1) is also referred to as symmetric Schur decomposition. Further on, without loss of generality, we will assume that the eigenvalues of A are ordered such that ${\lambda }_1\ge {\lambda }_2\ge \dots \ge {\lambda }_n$.

If the matrix A is subject to a Hermitian perturbation $\delta A\in {\mathbb{C}}^{n\times n},\delta A^H=\delta A$, then instead of the decomposition (1) we have the new diagonal decomposition

$$\tilde{\Lambda }={\tilde{U}}^H\tilde{A}\tilde{U}=\mathrm{diag}({\tilde{\lambda }}_1,{\tilde{\lambda }}_2,\dots ,{\tilde{\lambda }}_n),$$

(2)

where ${\tilde{\lambda }}_1\ge {\tilde{\lambda }}_2\ge \dots \ge {\tilde{\lambda }}_n$ are the perturbed eigenvalues and $\tilde{U}=U+\delta U=[{\tilde{U}}_1,{\tilde{U}}_2,\dots ,{\tilde{U}}_n]$ is the matrix of the perturbed eigenvectors.

The aim of the perturbation analysis is to determine bounds on the eigenvalue perturbations $\delta \Lambda =\mathrm{diag}\left(\delta {\lambda }_j\right),\delta {\lambda }_j={\tilde{\lambda }}_j-{\lambda }_j,i=1,2,\dots ,n$ and the corresponding eigenvector perturbations $\delta U_j={\tilde{U}}_j-U_j$. The perturbation problem for Hermitian matrices eigendecompositions is regular, if and only if the eigenvalues of A are distinct. An additional important problem is to determine the sensitivity of the one-dimensional invariant subspaces spanned by the eigenvectors $U_j,j=1,2,\dots ,n$. According to the Wielandt-Hofmann theorem [9], the eigenvalues of a Hermitian matrix are always well conditioned since their perturbations obey

$$|{\tilde{\lambda }}_j-{\lambda }_j{|\le \parallel \delta A\parallel }_2.$$

However, the eigenvectors can be very sensitive to perturbations of A, depending on the separation of the eigenvalues, see for instance [25].

According to the splitting operator method [13], it is possible to derive separately perturbation bounds on $\left|\delta U\right|$ and $\left|\delta \Lambda \right|$ which allows to obtain tighter perturbation bounds on the eigenvector. For this aim, it is appropriate to introduce the matrix

$$\delta W=U^H\delta U=\left[\begin{array}{cccc}{U}_1^H\delta U_1& U_1^H\delta U_2& \dots & U_1^H\delta U_n\\ U_2^H\delta U_1& U_2^H\delta U_2& \dots & U_2^H\delta U_n\\ ⋮& ⋮& \ddots & ⋮\\ U_n^H\delta U_1& U_n^H\delta U_2& \dots & U_n^H\delta U_n\end{array}\right]\in {\mathbb{C}}^{n\times n},$$

which is unitary equivalent to the unknown matrix $\delta U$. Finding estimates of the entries of this matrix allows to find sharp bounds on the entries of $\delta U$ thanks to the orthogonality of the matrix U.

Further on, we shall use the perturbation parameter vector

$$\begin{array}{ccc}\hfill x& =& \mathrm{vec}\left(\mathrm{Low}\right(\delta W\left)\right)\hfill \\ & =& {\left[U_2^H\delta U_1\dots U_n^H\delta U_1|U_3^H\delta U_2\dots U_n^H\delta U_2|\dots |U_n^H\delta U_{n-1}\right]}^T\in {\mathbb{C}}^{\nu },\nu =n(n-1)/2,\hfill \end{array}$$

where the components of x are the entries of the strictly lower triangular part of $\delta W$. Using this vector, it is convenient to find bounds on various quantities related to the perturbation analysis of A.

The paper is organized as follows. In sect. Section 2, we derive asymptotic bounds on the perturbation parameters which are then used in the perturbation analysis of eigenvalues and eigenvectors of a symmetric matrix. The eigenvector perturbation bounds are used to find asymptotic bounds on the perturbations of the invariant subspaces spanned by the eigenvectors. In sect. Section 3, we briefly present some results concerning the determination of realistic probabilistic bounds on the entries of a random matrix using its Frobenius norm. We describe the implementation of such bounds in the derivation of probabilistic perturbation bounds on the eigenvector entries, eigenvalues and invariant subspaces. We illustrate the theoretical results in Sect. Section 4 by presenting two numerical experiments with matrices of order 5000 demonstrating the performance of the derived bounds. Some conclusions are made in sect. Section 5.

All computations in the paper are done with MATLAB® Version 9.9 (R2020b) [17] using IEEE double precision arithmetic on a machine equipped with an 12th Gen Intel(R) Core(TM) i5-1240P CPU running at 1.70 GHz and with 32 GB of RAM. M-files implementing the perturbation bounds described in the paper can be obtained from the authors.

2. Asymptotic Perturbation Bounds

2.1. Asymptotic Bounds for the Perturbation Parameters

Theorem 1.

Let $A\in {\mathbb{C}}^{n\times n}$ be a Hermitian matrix ($A=A^H$) with distinct eigenvalues ${\lambda }_1,{\lambda }_2,\dots ,{\lambda }_n$, which is decomposed as $A=U\Lambda U^H$, where U is unitary and Λ is diagonal. Assume that the matrix A is perturbed by a Hermitian perturbation $\delta A$, so that

$$\tilde{A}=A+\delta A=\tilde{U}\tilde{\Lambda }{\tilde{U}}^H.$$

Denote by $F=-U^H\delta AU$ the transformed perturbation matrix and construct the vector

$$f=\mathrm{vec}\left(\mathrm{Low}\left(F\right)\right)\in {\mathbb{C}}^{\nu },f_{\ell }=-U_i^H\delta A{U}_j,\ell =i+(j-1)n-j(j+1)/2,1\le j<i\le n.$$

Then the perturbation parameters

$$x_{\ell }=U_i^H\delta U_j,\ell =i+(j-1)n-j(j+1)/2,1\le j<i\le n,$$

satisfy the equation

$$Mx=f+{\Delta }^x,$$

(3)

where $M\in {\mathbb{C}}^{\nu \times \nu }$ is a diagonal matrix of the form

with ${\mu }_{ij}={\lambda }_i-{\lambda }_j$ and ${\Delta }^x\in {\mathbb{C}}^{\nu }$ is a vector containing second order terms in $\delta U_i,\delta U_j$.

Proof. From (1) and (2) we have that

$${\tilde{U}}_i^H(A+\delta A){\tilde{U}}_j=U_i^{H}A{U}_j=0,j=1,2,\dots ,n-1,i=j+1,j+2,\dots ,n.$$

Hence

$$U_i^{H}A\delta U_j+\delta U_i^{H}A{U}_j+\delta U_i^{H}A\delta U_j=-{\tilde{U}}_i^H\delta A{\tilde{U}}_j$$

(4)

From the diagonal decomposition (1) we have that

$$\begin{array}{ccc}\hfill U_i^{H}A& =& {\lambda }_i{U}_i^H,i=1,2,\dots ,n,\hfill \\ \hfill A{U}_j& =& {\lambda }_j{U}_j,j=1,2,\dots ,n.\hfill \end{array}$$

(5)

The substitution of these expressions in (4) gives

$${\lambda }_i{U}_i^H\delta U_j+{\lambda }_j\delta U_i^H{U}_j+\delta U_i^{H}A\delta U_j=-{\tilde{U}}_i^H\delta A{\tilde{U}}_j$$

(6)

Since the matrices U and $\tilde{U}$ are unitary, it follows that

$$U^H\delta U=-\delta U^{H}U-\delta U^H\delta U$$

and hence

$$\delta U_i^H{U}_j=-U_i^H\delta U_j-\delta U_i^H\delta U_j,i=1,2,\dots ,n-1,j=i,i+1,\dots ,n.$$

(7)

Replacing this expression for $\delta U_i^H{U}_j$ into (6), we obtain that

$$({\lambda }_i-{\lambda }_j)U_i^H\delta U_j={\lambda }_j\delta U_i^H\delta U_j-\delta U_i^{H}A\delta U_j-{\tilde{U}}_i^H\delta A{\tilde{U}}_j$$

(8)

or

$$\begin{array}{ccc}\hfill ({\lambda }_i-{\lambda }_j)U_i^H\delta U_j& =& {\lambda }_j\delta U_i^H\delta U_j-\delta U_i^{H}A\delta U_j-U_i^H\delta A{U}_j-U_i^H\delta A\delta U_j\hfill \\ & & -\delta U_i^H\delta A{U}_j-\delta U_i^H\delta A\delta U_j.\hfill \end{array}$$

(9)

The expression (9) can be rewritten as a system of $\nu =n(n-1)/2$ nonlinear algebraic equations for the unknown quantities

$$x_{\ell }=U_i^H\delta U_j,\ell =i+(j-1)n-j(j+1)/2,1\le j<i\le n$$

in the form (3), where the component ${\Delta }_{\ell }^x,\ell =i+(j-1)n-j(j+1)/2$ of the nonlinear term ${\Delta }^x$ is equal to

$$\begin{array}{ccc}\hfill {\Delta }_{\ell }^x& =& {\lambda }_j\delta U_i^H\delta U_j-\delta U_i^{H}A\delta U_j-U_i^H\delta A\delta U_j-\delta U_i^H\delta A{U}_j\hfill \\ & & -\delta U_i^H\delta A\delta U_j.\hfill \end{array}$$

(10)

❒

We note that the solution of the system of equations (3) does not require the explicit formation of the matrix M. In particular, the asymptotic bound $x_{\ell }^{lin}$ in the inequality

$$|x_{\ell }|⪯x_{\ell }^{lin},\ell =1,2,\dots ,\nu ,$$

can be obtained easily by using the expression

$$x_{\ell }^{lin}={\parallel \delta A\parallel }_F/|{\lambda }_i-{\lambda }_j|,\ell =i+(j-1)n-j(j+1)/2,1\le j<i\le n.$$

(11)

Example 1.

Consider the $5\times 5$ matrix

$$A=\left[\begin{array}{ccccc}\hfill 2.401616& \hfill 1.198384& \hfill -0.398384& \hfill 0.402384& \hfill -0.398424\\ \hfill 1.198384& \hfill 2.201616& \hfill -0.001616& \hfill -0.002384& \hfill -0.001576\\ \hfill -0.398384& \hfill -0.001616& \hfill 1.801616& \hfill -0.797616& \hfill 0.801576\\ \hfill 0.402384& \hfill -0.002384& \hfill -0.797616& \hfill 1.803616& \hfill -0.797576\\ \hfill -0.398424& \hfill -0.001576& \hfill 0.801576& \hfill -0.797576& \hfill 1.801636\end{array}\right].$$

The eigenvalues of this matrix are

$${\lambda }_1=4.0,{\lambda }_2=3.0,{\lambda }_3=1.01,{\lambda }_4=1.0001,{\lambda }_5=1.0.$$

Note the closeness of the last three eigenvalues which prompts that the corresponding eigenvectors can be ill conditioned.

The perturbation matrix is taken as

$$\delta A={10}^{-9}\times (\delta A_0+\delta A_0^T),$$

where $\delta A_0$ is a matrix whose elements are random normal numbers.

For this perturbation problem, the matrix $M^{-1}$, which determines the perturbation parameter vector x in (3) has a 2-norm equal to $1\times {10}^4$, which confirms that the problem is relatively ill conditioned.

The exact perturbation parameters $x_{\ell },\ell =1,2,\dots ,\nu$ and their asymptotic approximations $x_{\ell }^{lin}$ computed by using (11), are shown to eight decimal digits in Table 1. Note the good coherence between the magnitude of the corresponding elements of both vectors. The differences between the values of $x_{\ell }^{lin}$ and $x_{\ell }$ are due to the bounding of the elements of the vector f by the value of ${\parallel \delta A\parallel }_F$ and the neglecting of the nonlinear term ${\delta }^x$.

2.2. Asymptotic Componentwise Eigenvector Bounds

Theorem 2.

Under the conditions of Theorem 1, a strict asymptotic bound on the perturbation of each eigenvector $U_j,j=1,2,\dots ,n$ of A under a perturbation $\delta A$, is given by

$$|{\tilde{U}}_j-U_j|⪯\delta U_j^{lin}=\left|U\right|W_j^{lin},$$

where

$$\begin{array}{ccc}\hfill \delta W^{lin}& =& \left[\begin{array}{cccccc}0& x_1^{lin}& x_2^{lin}& \dots & x_{n-2}^{lin}& x_{n-1}^{lin}\\ x_1^{lin}& 0& x_n^{lin}& \dots & x_{2n-4}^{lin}& x_{2n-3}^{lin}\\ x_2^{lin}& x_n^{lin}& 0& \dots & x_{3n-7}^{lin}& x_{3n-6}^{lin}\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ x_{n-2}^{lin}& x_{2n-4}^{lin}& x_{3n-7}^{lin}& \dots & 0& x_{\nu }^{lin}\\ x_{n-1}^{lin}& x_{2n-3}^{lin}& x_{3n-6}^{lin}& \dots & x_{\nu }^{lin}& 0\end{array}\right]\in {\mathbb{C}}^{n\times n},\hfill \end{array}$$

and the elements of the vector $x^{lin}={[x_1^{lin},x_2^{lin},\dots ,x_{\nu }^{lin}]}^T$ are determined form

$$x_{\ell }^{lin}={\parallel \delta A\parallel }_F/|{\lambda }_i-{\lambda }_j|,\ell =i+(j-1)n-j(j+1)/2,1\le j<i\le n.$$

Proof. Consider the auxiliary matrix

$$\delta W=U^H\delta U:=[\delta W_1,\delta W_2,\dots ,\delta W_n],\delta W_j\in {\mathbb{C}}^n$$

introduced above. As already noted, the strictly lower part of this matrix contains elements of the form

$$U_i^H\delta U_j,j=1,2,\dots ,n-1,i=j+1,j+1,\dots ,n$$

which can be substituted by the corresponding elements $x_{\ell },\ell =i+(j-1)n-j(j+1)/2$ of the vector x. The elements of the strictly upper part of $\delta W$ are of the form

$$U_i^H\delta U_j,i=1,2,\dots ,n-1,j=i+1,i+2,\dots ,n$$

which, according to the unitary condition (7), can be represented as

$$U_i^H\delta U_j=\delta U_i^H{U}_j-\delta U_i^H\delta U_j$$

or

$$U_i^H\delta U_j=-\overline{U_j^H\delta U_i}-\delta U_i^H\delta U_j,$$

(12)

where the term $\overline{U_j^H\delta U_i},j>i$ is the conjugate value of the element $x_{\ell }$. Hence the matrix $\delta W$ can be represented as

$$\delta W=\delta V-\delta D-\delta Y,$$

(13)

where

$$\begin{array}{ccc}\hfill \delta V& =& \left[\begin{array}{ccccc}0& -{\overline{x}}_1& -{\overline{x}}_2& \dots & -{\overline{x}}_{n-1}\\ x_1& 0& -{\overline{x}}_n& \dots & -{\overline{x}}_{2n-3}\\ x_2& x_n& 0& \dots & -{\overline{x}}_{3n-6}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ x_{n-1}& x_{2n-3}& x_{3n-6}& \dots & 0\end{array}\right]\in {\mathbb{C}}^{n\times n},\hfill \end{array}$$

and the matrices

$$\delta D=\left[\begin{array}{cccc}\delta U_1^H\delta U_1/2& 0& \dots & 0\\ 0& \delta U_2^H\delta U_2/2& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & \delta U_n^H\delta U_n/2\end{array}\right]\in {\mathbb{C}}^{n\times n},$$

$$\delta Y=\left[\begin{array}{ccccc}0& \delta U_1^H\delta U_2& \delta U_1^H\delta U_3& \dots & \delta U_1^H\delta U_n\\ 0& 0& \delta U_2^H\delta U_3& \dots & \delta U_2^H\delta U_n\\ 0& 0& 0& \dots & \delta U_3^H\delta U_n\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \dots & \delta U_{n-1}^H\delta U_n\end{array}\right]\in {\mathbb{C}}^{n\times n}$$

contain second order terms in $\delta U_j,j=1,2,\dots ,n$. Note that for definiteness, we have assumed that $\tilde{U}$ is chosen so that the diagonal elements of $\delta D$ are real and

$$U_i^H\delta U_i+\overline{U_i^H\delta U_i}=2U_i^H\delta U_i=-\delta U_i^H\delta U_j,$$

since the condition (12) does not restrict the imaginary part of $U_i^H\delta U_i$.

According to (13), the matrix $\left|\delta W\right|$ can be estimated as

$$\left|\delta W\right|⪯\left|\delta V\right|+{\delta }^W,$$

(14)

where

$${\delta }^W=\left|\delta D\right|+\left|\delta Y\right|$$

contains second order terms in the elements of x. Thus, an asymptotic (linear) approximation of the matrix $\left|\delta U\right|$ can be determined as

$$|\delta U^{lin}|⪯|U\left|\right|U^H\delta U|=|U|\delta W^{lin}.$$

(15)

❒

Note that the bounds derived are also valid for real symmetric matrices.

Example 2.

For the same matrix A and perturbation $\delta A$ as in Example 1, the absolute values of the exact changes of the entries of the matrix U and their linear estimates, obtained by using (15) are, respectively

$$\left|\delta U\right|={10}^{-4}\times \left[\begin{array}{ccccc}\hfill 0.0000045& \hfill 0.0000028& \hfill 0.0025383& \hfill 0.0928278& \hfill 0.0953600\\ \hfill 0.0000037& \hfill 0.0000020& \hfill 0.0025204& \hfill 0.0928299& \hfill 0.0953488\\ \hfill 0.0000025& \hfill 0.0000070& \hfill 0.0010687& \hfill 0.0928181& \hfill 0.1381715\\ \hfill 0.0000016& \hfill 0.0000139& \hfill 0.0025323& \hfill 0.0942843& \hfill 0.0904961\\ \hfill 0.0000071& \hfill 0.0000066& \hfill 0.0023290& \hfill 0.1407080& \hfill 0.0953601\end{array}\right]$$

and

$$\delta U^{lin}={10}^{-4}\times \left[\begin{array}{ccccc}\hfill 0.0000879& \hfill 0.0001319& \hfill 0.0088754& \hfill 0.4437818& \hfill 0.4438262\\ \hfill 0.0001099& \hfill 0.0001099& \hfill 0.0088791& \hfill 0.4437855& \hfill 0.4438298\\ \hfill 0.0000952& \hfill 0.0001209& \hfill 0.0110648& \hfill 0.4437745& \hfill 0.6634911\\ \hfill 0.0000953& \hfill 0.0001210& \hfill 0.0088680& \hfill 0.4459712& \hfill 0.4460378\\ \hfill 0.0000952& \hfill 0.0001209& \hfill 0.0110869& \hfill 0.6634467& \hfill 0.4438189\end{array}\right]$$

For all $i,j=1,2,\dots ,5$ it is fulfilled that $|\delta U_{i,j}|<\delta U_{i,j}^{lin}$.

2.3. Eigenvalue Sensitivity

For the changes in the elements of the perturbed diagonal form $\tilde{\Lambda }$ one has the following expressions

$$\delta {\lambda }_i={\tilde{\lambda }}_i-{\lambda }_i={\tilde{U}}_i^H(A+\delta A){\tilde{U}}_i-U_i^{H}A{U}_i,i=1,2,\dots ,n.$$

Hence

$$\begin{array}{ccc}\hfill \delta {\lambda }_i& =& U_i^{H}A\delta U_i+\delta U_i^{H}A{U}_i+\delta U_i^{H}A\delta U_i+U_i^H\delta A{U}_j\hfill \end{array}$$

(16)

$$\begin{array}{ccc}& & +U_i^H\delta A\delta U_i+\delta U_i^H\delta A{U}_i+\delta U_i^H\delta A\delta U_i.\hfill \end{array}$$

(17)

Thus, we obtain

$$\delta {\lambda }_i=U_i^H\delta A{U}_i+{\delta }_i^d,i=1,2,\dots ,n,$$

(18)

where, taking into account equations (5), we have that

$${\delta }_i^d={\lambda }_i(U_i^H\delta U_i+\delta U_i^H{U}_i)+\delta U_i^{H}A\delta U_i+U_i^H\delta A\delta U_i+\delta U_i^H\delta A{U}_i+\delta U_i^H\delta A\delta U_i.$$

According to (7), it is fulfilled that

$$\delta U_i^H{U}_j+U_i^H\delta U_j=-\delta U_i^H\delta U_j,$$

which gives

$${\delta }_i^d=-{\lambda }_i\delta U_i^H\delta U_i+\delta U_i^{H}A\delta U_i+U_i^H\delta A\delta U_i+\delta U_i^H\delta A{U}_i+\delta U_i^H\delta A\delta U_i.$$

(19)

The quantity ${\delta }_i^d$ contains second order terms in $\delta U_i$.

If the perturbation $\delta A$ is known, expression (19) allows to obtain bound on ${\delta }_i^d$ using the eigenvector bound (15).

In the asymptotic eigenvalue analysis the higher order terms are neglected and one has

$$\delta {\lambda }_i^{lin}=U_i^H\delta A{U}_i,i=1,2,\dots ,n.$$

Hence,

$$|\delta {\lambda }_i|\le \delta {\lambda }^{lin}={\parallel \delta A\parallel }_2,i=1,2,\dots ,n$$

which reduces to the well known corollary of the Wielandt-Hoffman theorem [9].

2.4. Sensitivity of One Dimensional Invariant Subspaces

The estimate of the eigenvector perturbation $\delta U_j$ can be used to find an estimate of the one dimensional (simple) invariant subspace, associated with the eigenvector $U_j$.

Consider the one dimensional invariant subspace ${\mathcal{X}}_j=\mathcal{R}\left(U_j\right),j=1,2,\dots ,n$. The sensitivity of this subspace is measured by the angle between the perturbed and unperturbed invariant subspace. Since

$${\tilde{U}}_j=U_j+\delta U_j,$$

we have that [5,26]

$$sin\left(\mathrm{\Theta }({\tilde{\mathcal{X}}}_j,{\mathcal{X}}_j)\right)=\parallel {U_j^{\perp }}^H{\tilde{U}}_j{\parallel }_2={\parallel {U_j^{\perp }}^H\delta U_j\parallel }_2,$$

(20)

where

$$U_j^{\perp }=\left[U_1,U_2,\dots ,U_{j-1},U_{j+1},\dots ,U_n\right]\in {\mathbb{C}}^{n\times (n-1)}$$

is the orthogonal complement of $U_j,{U_j^{\perp }}^H{U}_j=0_{(n-1)\times 1}$.

Equation (20) shows that the sensitivity of the one dimensional invariant subspace $\mathcal{X}$ is connected to the values of the perturbation parameters $x_{\ell }=U_i^H\delta U_j$. Consequently, if bounds on the perturbation parameters are known, it is possible to find the sensitivity estimates of all invariant subspaces. More specifically, we have that

$${U_1^{\perp }}^H\delta U_1=\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_{n-1}\end{array}\right],{U_2^{\perp }}^H\delta U_2=\left[\begin{array}{c}-{\overline{x}}_1\\ x_n\\ ⋮\\ x_{2n-3}\end{array}\right],\dots ,{U_n^{\perp }}^H\delta U_n=\left[\begin{array}{c}-{\overline{x}}_{n-1}\\ -{\overline{x}}_{2n-3}\\ -{\overline{x}}_{3n-6}\\ ⋮\\ -{\overline{x}}_{\nu }\end{array}\right].$$

Then we obtain that the angle between the jth perturbed and unperturbed invariant subspace has the following asymptotic estimate

$$\mathrm{\Theta }({\tilde{\mathcal{X}}}_j,{\mathcal{X}}_j)\le {\mathrm{\Theta }}^{lin}({\tilde{\mathcal{X}}}_j,{\mathcal{X}}_j)=arcsin(\parallel \delta W^{lin}(1:n,j){\parallel }_2),j=1,2,\dots ,n.$$

(21)

We note that this estimate produces practically the same results as the linear bound derived in [3,28,29].

Example 3.

For the same matrix and perturbation, as in the previous examples, the exact angles between the perturbed and unperturbed one-dimensional invariant subspaces and their linear estimates, are shown in Table 2.

3. Probabilistic Asymptotic Bounds

The numerical experiments with symmetric matrices show that the estimates obtained by using Theorem 2 can become very pessimistic with the increasing of matrix order. For instance, if the matrix is of order 2000, then the ratio between the bound (15) and the actual values of the entries of $\left|\delta U\right|$ may become of order ${10}^5$ which makes the computed bound useless. A further reduction of the perturbation bounds can be achieved by implementing probabilistic perturbation bounds which for large n allow to obtain sufficiently close bounds with a high probability. For this aim we will make use of the probabilistic matrix bounds proposed in [23,24] that are based on the Markoff inequality [19].

Consider briefly the essence of the approach presented in [23]. The aim is to reduce the size of entries of an estimate $\Delta A$ of the matrix perturbation $\delta A$ at the price that some entries of $|\Delta A|$ are smaller then the corresponding entries of $\left|\delta A\right|$. We have the following result.

Theorem 3.

Let $\Delta A=[\Delta a_{ij}]$ be an estimate of the $m\times n$ random perturbation $\delta A$ and $\mathcal{P}\left\{\right|\delta a_{ij}|<\Delta a_{ij}\}$ be the probability that $|\delta a_{ij}|<\Delta a_{ij}$. If the entries of $\Delta A$ are chosen as

$$\Delta a_{ij}={\displaystyle \frac{{\parallel \delta A\parallel }_F}{\mathrm{\Xi }}},$$

where

$$\mathrm{\Xi }=(1-{\mathcal{P}}^{ref})\sqrt{mn},$$

(22)

and $0<{\mathcal{P}}^{ref}<1$ is a desired probability, then it is fulfilled that

$$\mathcal{P}\left\{\right|\delta a_{ij}|<\Delta a_{ij}\}\ge {\mathcal{P}}^{ref},i=1,2,\dots ,m,j=1,2,\dots ,n.$$

If the number $mn$ is sufficiently large so that $\mathrm{\Xi }>1$, Theorem 3 allows to decrease the mean value of the bound $\Delta A$ and hence the magnitude of its entries by the scaling factor $\mathrm{\Xi }$, choosing the desired probability ${\mathcal{P}}^{ref}$ less than 1. Note that the probability bound produced by the Markoff inequality is conservative, the actual results usually being much better than the results predicted by the quantity ${\mathcal{P}}^{ref}$. This is due to the fact that the Markoff inequality is valid for the worst possible distribution of the random variables, but in the given case we do not impose a restriction on the probability distribution of the entries of $\delta A$.

According to Theorem 3, the using of the scaling factor (22) guarantees that the inequality

$$|\delta a_{ij}|<\Delta a_{ij}$$

holds for each i and j with probability no less than ${\mathcal{P}}^{ref}$.

Using Theorem 3, a probability bound on the perturbation parameters $|x_{\ell }|,\ell =1,2,\dots ,\nu$ can be determined by the following theorem.

Theorem 4.

If the asymptotic estimate of the perturbation parameter vector x is chosen as

$$x_{\ell }^{est}={\displaystyle \frac{{\parallel \delta A\parallel }_F}{\mathrm{\Xi }|{\lambda }_i-{\lambda }_j|}},\ell =i+(j-1)n-j(j+1)/2,1\le j<i\le n,$$

(23)

where Ξ is determined according to

$$\mathrm{\Xi }=n(1-{\mathcal{P}}^{ref}),$$

(24)

then

$$\mathcal{P}\left\{\right|x_{\ell }|\le \parallel x^{est}{\parallel }_2\}\le {\mathcal{P}}^{ref}.$$

(25)

The inequality (25) shows that the probability estimate of the component $|x_{\ell }|$ can be determined, if in the linear estimate (11) we replace the perturbation norm ${\parallel \delta A\parallel }_F$ by the probability estimate ${\parallel \delta A\parallel }_F/\mathrm{\Xi }$, where the scaling factor $\mathrm{\Xi }$ is taken as shown in (22) for a specified probability ${\mathcal{P}}^{ref}$.

The probabilistic perturbation bounds on the elements of the vector x allow to find probabilistic bounds on the elements of $\delta U$.

Theorem 5.

An asymptotic bound with probability ${\mathcal{P}}^{ref}$ on the perturbation of each eigenvector $U_j,j=1,2,\dots ,n$ of A under a perturbation $\delta A$, is given by

$$|{\tilde{U}}_j-U_j|⪯|\delta U_j^{est}|=|U|W_j^{est},$$

where

$$\begin{array}{ccc}\hfill \delta W^{est}& =& \left[\begin{array}{cccccc}0& x_1^{est}& x_2^{est}& \dots & x_{n-2}^{est}& x_{n-1}^{est}\\ x_1^{est}& 0& x_n^{est}& \dots & x_{2n-4}^{est}& x_{2n-3}^{est}\\ x_2^{est}& x_n^{est}& 0& \dots & x_{3n-7}^{est}& x_{3n-6}^{est}\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ x_{n-2}^{est}& x_{2n-4}^{est}& x_{3n-7}^{est}& \dots & 0& x_{\nu }^{est}\\ x_{n-1}^{est}& x_{2n-3}^{est}& x_{3n-6}^{est}& \dots & x_{\nu }^{est}& 0\end{array}\right]\in {\mathbb{C}}^{n\times n},\hfill \end{array}$$

and the elements of the vector $x^{est}={[x_1^{est},x_2^{est},\dots ,x_{\nu }^{est}]}^T$ are determined form

$$x_{\ell }^{est}={\displaystyle \frac{{\parallel \delta A\parallel }_F}{\mathrm{\Xi }|{\lambda }_i-{\lambda }_j|}},\ell =i+(j-1)n-j(j+1)/2,1\le j<i\le n.$$

The probabilistic bound on $\delta A$ produces the following asymptotic perturbation bounds on the eigenvalues of A:

$$|\delta {\lambda }_i|\le \delta {\lambda }_i^{est}={\parallel \delta A\parallel }_F/\mathrm{\Xi },i=1,2,\dots ,n.$$

Finally, we obtain that the angle between the jth perturbed and unperturbed one-dimensional invariant subspace satisfies the probabilistic asymptotic estimate

$${\mathrm{\Theta }}_j({\tilde{\mathcal{X}}}_j,{\mathcal{X}}_j)\le {\mathrm{\Theta }}_j^{est}({\tilde{\mathcal{X}}}_j,{\mathcal{X}}_j)=arcsin(\parallel \delta W^{est}(1:n,j){\parallel }_2),j=1,2,\dots ,n.$$

(26)

4. Numerical Experiments

In this section, we present two numerical experiments for determining the perturbation bounds for the eigendecompositions of matrices of order 5000. The matrices are constructed in the form $A=VD{V}^T$ where V is orthogonal and D is a diagonal matrix containing the desired eigenvalues. The perturbation $\delta A$ is taken as ${10}^{-9}\times (\delta A_0+\delta A_0^T)$, where $\delta A_0$ is a matrix whose elements are normal random numbers.

Example 4.

In this examples the eigenvalues of A are taken uniformly between 1 and 5.999 with an increment equal to 0.001, which results in matrix $M^{-1}$ with a 2-norm equal to ${10}^3$.

In Figure 1 we show the entries of the exact perturbation $|\delta U=\tilde{U}-U|$ and the entries of the asymptotic $\delta U^{lin}$ and probabilistic $\delta U^{est}$ bounds. The desired lower bound probability is set to $P^{ref}=0.8$ and the scaling factor in finding the probabilistic estimates is equal to 1000. In a result, the actual probability is equal to $100\%$, since the bounds of all entries of $\delta U^{est}$ exceed the corresponding entries of $\left|\delta U\right|$. As already mentioned, this phenomenon is due to the conservatism of the Markoff inequality. The eigenvalue perturbations $\delta {\lambda }_i={\tilde{\lambda }}_i-{\lambda }_i$ along with their asymptotic $\delta {\lambda }^{lin}=2.8287199\times {10}^{-6}$ and probability $\delta {\lambda }^{est}=7.0717998\times {10}^{-9}$ estimates, are shown in Figure 2. There are no eigenvalues whose probability perturbation bound is smaller than the actual $\delta \lambda$, which gives $100\%$ accuracy instead of the reference value $80\%$. Due to the equal space between the eigenvalues, in the given case the sensitivity of all invariant subspaces is the same (Figure 3). Clearly, the probabilistic estimates are much closer to the actual perturbations than the linear deterministic bounds.

Example 5.

In this example the eigenvalues of A are taken as ${\lambda }_i=1+500/i$ and are spread between $1.1$ and 501, the norm of the matrix $M^{-1}$ being equal to $4.999\times {10}^4$. The eigenvalues of A are shown in Figure 4.

In Figure 5, Figure 6 and Figure 7 we show the matrix $|{\delta }_U|$, the eigenvalue perturbations and the angles between the perturbed and unperturbed invariant subspaces, respectively, as well as their asymptotic and probabilistic bounds for $P^{ref}=80\%$. Due to the decreasing distance between the eigenvalues, the sensitivity of the invariant subspace is increasing with the sequence number i. The actual probabilistic bounds for the eigenvector matrix and for the eigenvalues are again equal to $100\%$ as in the previous example.

Note that in both examples, the linear estimates reflect correctly the changes of the corresponding exact quantities.

5. Conclusions

In this paper, we present strict asymptotic perturbation bounds of the eigenvalues, eigenvectors and invariant subspaces of symmetric and Hermitian matrices. The implementation of the splitting operator method makes possible to unify the analysis with the perturbation analysis of the Schur form [21], the generalized Schur form [32], the singular value decomposition [1] and the QR decomposition of a matrix [22]. Since the deterministic bounds derived can be conservative especially for high-order matrices, we present probabilistic bounds based on the Markoff inequality. The probabilistic bounds are found easily from the asymptotic bounds and allow to decrease significantly the perturbation bounds with a guaranteed probability which is near to 1 in practice. An attractive future of the bounds for symmetric matrices is that these bounds can be determined with low requirements in respect to the necessary memory and volume of computations.

6. Notation

$\mathbb{C}$,	the set of complex numbers;
${\mathbb{C}}^{n\times m}$,	the space of $n\times m$ complex matrices;
$A=\left[a_{ij}\right]$,	a matrix with entries $a_{ij}$;
$A_j$,	the jth column of A;
$A_{i,1:n}$,	the ith row of an $m\times n$ matrix A;
$A_{1:m,j}$,	the jth column of an $m\times n$ matrix A;
$\mathrm{Low}\left(A\right)$,	the strictly lower triangular part of A;
$\left|A\right|$,	the matrix of absolute values of the elements of A;
$A^H$,	the Hermitian transposed of A;
$0_{m\times n}$,	the zero $m\times n$ matrix;
$I_n$,	the unit $n\times n$ matrix;
$\delta A$,	the perturbation of A;
${\parallel A\parallel }_2$,	the spectral norm of A;
${\parallel A\parallel }_F$,	the Frobenius norm of A;
$:=$,	equal by definition;
⪯,	relation of partial order. If $a,b\in {\mathbb{R}}^n$, then $a⪯b$ means
	$a_i\le b_i,i=1,2,\dots ,n$;
$\mathcal{X}=\mathcal{R}\left(X\right)$,	the subspace spanned by the columns of X;
$U^{\perp }$,	the orthogonal complement of U, $U^H{U}^{\perp }=0$;
❒,	the end of a proof.