Componentwise Perturbation Analysis of the Singular Value Decomposition of a Matrix

1. Introduction

As it is known [5, Ch. 2], [13, Ch. 1], [4, Ch, 2], each real $m\times n,m\ge n$ matrix A can be represented by the singular value decomposition (SVD) in the factorized form

$$A=U\Sigma V^T,$$

(1)

where the $m\times m$ matrix U and the $n\times n$ matrix V are orthogonal and the $m\times n$ matrix $\Sigma$ is diagonal:

$$\Sigma =\left[\begin{array}{c}{\Sigma }_n\\ 0\end{array}\right],{\Sigma }_n=\mathrm{diag}\left({\sigma }_1,{\sigma }_2,\dots ,{\sigma }_n\right).$$

The numbers ${\sigma }_i\ge 0$ are called singular values of the matrix A. The columns of

$$U:=[u_1,u_2,\dots ,u_m],u_j\in {\mathbb{R}}^m,$$

are called left singular vectors and the columns of

$$V:=[v_1,v_2,\dots ,v_n],v_j\in {\mathbb{R}}^n$$

are the right singular vectors. The subspaces spanned by sets of left and right singular vectors are called left and right singular subspaces, respectively.

The usual assumption is that by an appropriate ordering of the columns $u_j$ and $v_j$, the singular values appear in the order

$${\sigma }_{max}:={\sigma }_1\ge {\sigma }_2\ge \dots \ge {\sigma }_n:={\sigma }_{min},$$

but in this paper we will not impose such a requirement.

The SVD has a lot of properties which make it an invaluable tool in matrix analysis and matrix computations, see the references cited above. Among them is the fact that the rank of A is equal to the number of its nonzero singular values as well as the equality ${\parallel A\parallel }_2={\sigma }_{max}$. The singular value decomposition has a long and interesting history which is described in [12].

In this paper we are interested in the case when the matrix A is subject to an additive perturbation $\delta A$. In such a case there exists another pair of orthogonal matrices $\tilde{U}$ and $\tilde{V}$ and a diagonal matrix $\tilde{\Sigma }$, such that

$$\tilde{A}:=A+\delta A=\tilde{U}\tilde{\Sigma }{\tilde{V}}^T.$$

(2)

The perturbation analysis of the singular value decomposition consists in determining the changes of the quantities related to the elements of the decomposition due to the perturbation $\delta A$. This includes determining bounds on the changes of the entries of the orthogonal matrices which reduce the original matrix to diagonal form and bounds on the perturbations of the singular values. Hence, the aim of the analysis is to find bounds on the sizes of $\delta U=\tilde{U}-U$, $\delta V=\tilde{V}-V$ and $\delta \Sigma =\tilde{\Sigma }-\Sigma$ as functions of the size of $\delta A$. According to the Weyl’s theorem [13, Ch. 1], we have that

$$|\delta {\sigma }_i|=|{\tilde{\sigma }}_i-{\sigma }_i{|\le \parallel \delta A\parallel }_2,i=1,2,\dots ,n$$

(3)

which shows that the singular values are perturbed by no more than the 2-norm of the perturbation of A, i.e., the singular values are always well conditioned. The SVD perturbation analysis is well defined if the matrix A is of full column rank n, i.e. ${\sigma }_{min}\ne 0$ since otherwwise the corresponding left singular vector is undetermined.

The size of the perturbations $\delta A$, $\delta U$, $\delta V$ and $\delta \Sigma$ is usually measured by using some of the matrix norms which leads to the so called normwise perturbation analysis. In several cases we are interested in the size of the perturbations of the individual entries of $\delta U$, $\delta V$ and $\delta \Sigma$, so that it is necessary to implement a componentwise perturbation analysis. This analysis has an advantage in the cases when the individual components of $\delta U$ and $\delta V$ differ very much in magnitude and the normwise estimates do not produce tight bounds on the perturbations.

The first results in the perturbation analysis of the singular value decomposition are obtained by Wedin [16] and Stewart [10], who developed estimates of the sensitivity of pairs of singular subspaces, see also [14, Ch. V]. Other results concerning the sensitivity of the SVD can be found in [3] and [15]. Several results concerning the sensitivity of the SVD are summarized in [6] and a survey on the perturbation theory of the singular value decomposition can be found in [11]. It should be pointed out that a componentwise perturbation analysis of the SVD apart from several results about the sensitivity of the singular values, is not available up to the moment.

In this paper we present a rigorous perturbation analysis of the orthogonal matrices, singular subspaces and singular values of a real matrix of full column rank. It is shown that the SVD perturbation problem is well posed only in case of distinct (simple) singular values. The analysis produces asymptotic (local) componentwise perturbation bounds of the entries of the orthogonal matrices U and V and of the singular values of the given matrix. Local bounds are derived for the sensitivity of a pair of singular subspaces measured by the angles between the unperturbed and perturbed subspaces. An iterative scheme is described to find global bounds on the respective perturbations and results of numerical experiments are presented. The analysis performed in the paper implements the same methodology as the one used previously in [8,9] to determine componentwise perturbation bounds of the Schur form and QR decomposition of a matrix. However, the SVD perturbation analysis has some distinctive features which makes it a self-dependent problem.

The paper is organized as follows. In Section 2 we derive the basic nonlinear algebraic equations used to perform the perturbation analysis of the SVD. After introducing in Section 3 the perturbation parameters that determine the perturbations of the matrices U and V, we derive a system of coupled equations for these parameters in Section 4. The solution of the equations for the first-order terms of the perturbation parameters allows to find asymptotic bounds on the parameters in Section 5, on the singular values in Section 6 and on the perturbations in the matrices U and V in Section 7. Using the bounds on the perturbation parameters, in Section 8 we derive bounds on the sensitivity of singular subspaces. In Section 9, we develop an iterative scheme for finding global bounds on the perturbations and in Section 10 we present the results of two higher order examples illustrating the proposed analysis. Some conclusions are drawn in Section 11.

2. Basic Equations

The perturbed singular value decomposition of A (2) can be written as

$${\tilde{U}}^T\tilde{A}\tilde{V}=\tilde{\Sigma },$$

(4)

where

$$\begin{array}{ccc}\hfill \tilde{U}& =& U+\delta U:=[{\tilde{u}}_1,{\tilde{u}}_2,\dots ,{\tilde{u}}_m],{\tilde{u}}_j\in {\mathbb{R}}^m,\hfill \\ \hfill \delta U& :=& [\delta u_1,\delta u_2,\dots ,\delta u_m],\delta u_j\in {\mathbb{R}}^m\hfill \\ \hfill {\tilde{u}}_j& =& u_j+\delta u_j,j=1,2,\dots ,m,\hfill \\ \hfill \tilde{V}& =& V+\delta V:=[{\tilde{v}}_1,{\tilde{v}}_2,\dots ,{\tilde{v}}_n],{\tilde{v}}_j\in {\mathbb{R}}^n\hfill \\ \hfill \delta V& :=& [\delta v_1,\delta v_2,\dots ,\delta v_n],\delta v_j\in {\mathbb{R}}^n\hfill \\ \hfill {\tilde{v}}_j& =& v_j+\delta v_j,j=1,2,\dots ,n,\hfill \end{array}$$

and

$$\begin{array}{ccc}\hfill \tilde{\Sigma }& =& \Sigma +\delta \Sigma =\left[\begin{array}{c}{\tilde{\Sigma }}_n\\ 0\end{array}\right],\hfill \\ \hfill {\tilde{\Sigma }}_n& =& {\Sigma }_n+\delta {\Sigma }_n=\mathrm{diag}({\tilde{\sigma }}_1,{\tilde{\sigma }}_2,\dots ,{\tilde{\sigma }}_n),\hfill \\ \hfill \delta {\Sigma }_n& =& \mathrm{diag}(\delta {\sigma }_1,\delta {\sigma }_2,\dots ,\delta {\sigma }_n),\hfill \\ \hfill {\tilde{\sigma }}_i& =& {\sigma }_i+\delta {\sigma }_i.\hfill \end{array}$$

Equation (4) is rewritten as

$$\delta U^{T}AV+U^{T}A\delta V+\delta F=\left[\begin{array}{c}\delta {\Sigma }_n\\ 0_{(m-n)\times n}\end{array}\right]+{\Delta }_0,$$

(5)

where $\delta F=U^T\delta AV$ and the matrix ${\Delta }_0$

$${\Delta }_0=-\delta U^{T}A\delta V-U^T\delta A\delta V-\delta U^T\delta AV-\delta U^T\delta A\delta V\in {\mathbb{R}}^{m\times n}$$

contains only higher-order terms in the elements of $\delta A$, $\delta U$ and $\delta V$.

Let the matrices U and $\delta U$ be divided as $U=[U_1,U_2],U_1\in {\mathbb{R}}^{m\times n}$ and $\delta U=[\delta U_1,\delta U_2],\delta U_1\in {\mathbb{R}}^{m\times n}$, respectively. Since the matrix A can be represented as $A=U_1{\Sigma }_n{V}^T$, the matrix $U_2$ is not well determined but should satisfy the orthogonality condition ${[U_1,U_2]}^T[U_1,U_2]=I_m$. The perturbation $\delta U_2$ is also undefined, so that we can bound only the perturbations of the entries in the first n columns of U, i.e., the entries of $\delta U_1$. Further on, we shall use (5) to determine componentwise bounds on $\delta U_1$, $\delta V$ and $\delta {\Sigma }_n$.

3. Perturbation Parameters and Perturbed Orthogonal Matrices

In the perturbation analysis of the SVD, it is convenient to find first componentwise bounds on the entries of the matrices $\delta W_U:=U^T\delta U_1$ and $\delta W_V:=V^T\delta V$, which are related to the corresponding perturbations $\delta U_1$ and $\delta V$ by orthogonal transformations. The implementation of the matrices $\delta W_U$ and $\delta W_V$ allows to find bounds on

$$\delta U_1=U\delta W_U$$

(6)

and

$$\delta V=V\delta W_V$$

(7)

using orthogonal transformations without increasing the norms of $\delta W_U$ and $\delta W_V$. This helps to determine bounds on $\delta U_1$ and $\delta V$ which are as tight as possible.

Consider first the matrix

$$\delta W_U=U^T\delta U_1=\left[\begin{array}{cccc}{u}_1^T\delta u_1& u_1^T\delta u_2& \dots & u_1^T\delta u_n\\ u_2^T\delta u_1& u_2^T\delta u_2& \dots & u_2^T\delta u_n\\ ⋮& ⋮& ⋮& ⋮\\ u_n^T\delta u_1& u_n^T\delta u_2& \dots & u_n^T\delta u_n\\ u_{n+1}^T\delta u_1& u_{n+1}^T\delta u_2& \dots & u_{n+1}^T\delta u_n\\ ⋮& ⋮& ⋮& ⋮\\ u_m^T\delta u_1& u_m^T\delta u_2& \dots & u_m^T\delta u_n\end{array}\right]\in {\mathbb{R}}^{m\times n}.$$

Further on, we shall use the vector of the subdiagonal entries of the matrix $\delta W_U$,

$$\begin{array}{ccc}\hfill x& :=& {[x_1,x_2,\dots ,x_p]}^T\in {\mathbb{R}}^p,\hfill \\ & =& {[\underset{m-1}{\underbrace{u_2^T\delta u_1,u_3^T\delta u_1,\dots ,u_m^T\delta u_1}},\underset{m-2}{\underbrace{u_3^T\delta u_2,\dots ,u_m^T\delta u_2}},\dots ,\underset{m-n}{\underbrace{u_{n+1}^T\delta u_n,\dots ,u_m^T\delta u_n}}]}^T,\hfill \end{array}$$

where

$$p=\sum _{i=1}^n(m-i)=n(n-1)/2+(m-n)n=n(2m-n-1)/2.$$

As it will become clear latter on, together with the orthogonality condition

$${[U_1+\delta U_1]}^T[U_1+\delta U_1]=I_n,$$

the vector x contains the whole information which is necessary to find the perturbation $\delta U_1$. This vector may be expressed as

$$x=\mathrm{vec}\left(\mathrm{Low}\left(\delta W_U\right)\right),$$

or, equivalently,

$$x={\Omega }_x\mathrm{vec}\left(\delta W_U\right),$$

where

$$\begin{array}{ccc}\hfill {\Omega }_x& :=& \left[\mathrm{diag}({\omega }_1,{\omega }_2,\dots ,{\omega }_n)\right]\in {\mathbb{R}}^{p\times mn},\hfill \\ \hfill {\omega }_i& :=& \left[0_{(m-i)\times i},I_{m-i}\right]\in {\mathbb{R}}^{(m-i)\times m},i=1,2,\dots ,n,\hfill \\ & & {\Omega }_x^T{\Omega }_x=I_p,{\parallel {\Omega }_x\parallel }_2=1\hfill \end{array}$$

is a matrix that “pulls out”the p elements of x from the $m.n$ elements of $\delta W_U$ (we consider $0_{0\times i}$ as a non-existing matrix).

We have that

$$\begin{array}{ccc}\hfill x_k& =& u_i^T\delta u_j,k=i+(j-1)m-\frac{j(j+1)}{2},\hfill \\ & & 1\le j\le n,j<i\le m,1\le k\le p.\hfill \end{array}$$

In a similar way, we introduce the vector of the subdiagonal entries of the matrix $\delta W_V=V^T\delta V$ (note that V is a square matrix),

$$\begin{array}{ccc}\hfill y& =& {[y_1,y_2,\dots ,y_q]}^T\in {\mathbb{R}}^q,\hfill \\ \hfill y_i& =& {[\underset{n-1}{\underbrace{v_2^T\delta v_1,v_3^T\delta v_1,\dots ,v_n^T\delta v_1}},\underset{n-2}{\underbrace{v_3^T\delta v_2,\dots ,v_n^T\delta v_2}},\dots ,\underset{1}{\underbrace{v_n^T\delta v_{n-1}}}]}^T,\hfill \end{array}$$

where $q=n(n-1)/2$. It is fulfilled that

$$y=\mathrm{vec}\left(\mathrm{Low}\left(\delta W_V\right)\right)$$

or, equivalently,

$$y={\Omega }_y\mathrm{vec}\left(\delta W_V\right),$$

where

$$\begin{array}{ccc}\hfill {\Omega }_y& :=& \left[\mathrm{diag}({\omega }_1,{\omega }_2,\dots ,{\omega }_n)\right]\in {\mathbb{R}}^{q\times n^2},\hfill \\ \hfill {\omega }_i& :=& \left[0_{(n-i)\times i},I_{n-i}\right]\in {\mathbb{R}}^{(n-i)\times n},i=1,2,\dots ,n,\hfill \\ & & {\Omega }_y^T{\Omega }_y=I_q,{\parallel {\Omega }_y\parallel }_2=1.\hfill \end{array}$$

In this case

$$\begin{array}{ccc}\hfill y_{\ell }& =& v_i^T\delta v_j,\ell =i+(j-1)n-\frac{j(j+1)}{2},\hfill \\ & & 1\le j\le n,j<i\le n-1,1\le \ell \le q.\hfill \end{array}$$

Further on the quantities $x_k,k=1,2,\dots ,p$ and $y_{\ell },\ell =1,2,\dots ,q$ will be referred to as perturbation parameters since they determine the perturbations $\delta U_1$ and $\delta V$, as well as the sensitivity of the singular values and singular subspaces.

Consider the matrix

$$\delta W_U=U^T\delta U_1:=[\delta {w_u}_1,\delta {w_u}_2,\dots ,\delta {w_u}_n],\delta {w_u}_j\in {\mathbb{R}}^m.$$

Using the vector x of the perturbation parameters, this matrix is written in the form

where the diagonal and super-diagonal entries that are not determined up to the moment, are highlighted in red boxes.

Consider how to determine the diagonal elements of the matrix $\delta W_U$,

$${d_u}_{jj}:=u_j^T\delta u_j,j=1,2,\dots ,n$$

by the elements of x. Since $\delta u_j^T{u}_j=u_j^T\delta u_j$, according to (9), we have that

$$2u_j^T\delta u_j=-\delta u_j^T\delta u_j,j=1,2,\dots ,n$$

or

$${d_u}_{jj}=-{\parallel \delta u_j\parallel }_2^2/2.$$

The above expression shows that ${d_u}_{jj}$ is always negative and is of second order of magnitude in $\parallel \delta u_j{\parallel }_2$.

Let us determine now the entries of the super-diagonal part of $W_U$. Since ${\tilde{U}}^T\tilde{U}=I_m$, it follows that

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

(8)

and

$$\delta u_i^T{u}_j=-u_i^T\delta u_j-\delta u_i^T\delta u_j,1\le i\le m,i<j\le m.$$

(9)

According to the orthogonality condition (9), the entries of the strictly upper triangular part of $\delta W_U$ can be represented as

$$u_i^T\delta u_j=-u_j^T\delta u_i-\delta u_i^T\delta u_j,1\le i\le n,i<j\le n.$$

Thus, the matrix $\delta W_u$ can be represented as the sum

$$\delta W_U=\delta Q_U+\delta D_U-\delta N_U,$$

(10)

where the matrix

$$\begin{array}{ccc}\hfill \delta Q_U& =& \left[\begin{array}{ccccc}0& -x_1& -x_2& \dots & -x_{n-1}\\ x_1& 0& -x_m& \dots & -x_{m+n-3}\\ x_2& x_m& 0& \dots & -x_{2m+n-6}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ x_{n-1}& x_{m+n-3}& x_{2m+n-6}& \dots & 0\\ x_n& x_{m+n-2}& x_{2m+n-5}& \dots & x_{(n-1)(2m-n)/2+1}\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ x_{m-1}& x_{2m-3}& x_{3m-6}& \dots & x_p\end{array}\right]\hfill \\ & :=& [\delta {q_u}_1,\delta {q_u}_2,\dots ,\delta {q_u}_n],\delta {q_u}_j\in {\mathbb{R}}^m\hfill \end{array}$$

has entries depending only on the perturbation parameters x,

$$\delta D_U=\left[\begin{array}{cccc}{d_u}_{11}& 0& \dots & 0\\ 0& {d_u}_{22}& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & {d_u}_{nn}\\ 0& 0& \dots & 0\\ ⋮& ⋮& ⋮& ⋮\\ 0& 0& \dots & 0\end{array}\right]\in {\mathbb{R}}^{m\times n},$$

and the matrix

$$\delta N_U=\left[\begin{array}{ccccc}0& \delta u_1^T\delta u_2& \delta u_1^T\delta u_3& \dots & \delta u_1^T\delta u_n\\ 0& 0& \delta u_2^T\delta u_3& \dots & \delta u_2^T\delta u_n\\ 0& 0& 0& \dots & \delta u_3^T\delta u_n\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \dots & \delta u_{n-1}^T\delta u_n\\ 0& 0& 0& \dots & 0\\ 0& 0& 0& \dots & 0\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ 0& 0& 0& \dots & 0\end{array}\right]\in {\mathbb{R}}^{m\times n}$$

contains second order terms in $\delta u_j,j=1,2,\dots ,n$.

Similarly, for the matrix

$$\delta W_V=V^T\delta V:=[\delta {w_v}_1,\delta {w_v}_2,\dots ,\delta {w_v}_n],\delta {w_v}_j\in {\mathbb{R}}^n,$$

like to the case of $\delta W_U$, it is possible to show that

$$\delta W_V=\delta Q_V+\delta D_V-\delta N_V,$$

(11)

where

$$\begin{array}{ccc}\hfill \delta Q_V& =& \left[\begin{array}{ccccc}0& -y_1& -y_2& \dots & -y_{n-1}\\ y_1& 0& -y_n& \dots & -y_{2n-3}\\ y_2& y_n& 0& \dots & -y_{3n-6}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ y_{n-1}& y_{2n-3}& y_{3n-6}& \dots & 0\end{array}\right]\hfill \\ & :=& [\delta {q_v}_1,\delta {q_v}_2,\dots ,\delta {q_v}_n],\delta {q_v}_j\in {\mathbb{R}}^n\hfill \end{array}$$

has elements depending only on the perturbation parameters y,

$$\delta D_V=\left[\begin{array}{cccc}{d_v}_{11}& 0& \dots & 0\\ 0& {d_v}_{22}& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & {d_v}_{nn}\end{array}\right]\in {\mathbb{R}}^{n\times n},$$

${d_v}_{jj}=v_j^T\delta v_j,j=1,2,\dots ,n$ and the matrix

$$\delta N_V=\left[\begin{array}{ccccc}0& \delta v_1^T\delta v_2& \delta v_1^T\delta v_3& \dots & \delta v_1^T\delta v_n\\ 0& 0& \delta v_2^T\delta v_3& \dots & \delta v_2^T\delta v_n\\ 0& 0& 0& \dots & \delta v_3^T\delta v_n\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \dots & \delta v_{n-1}^T\delta v_n\\ 0& 0& 0& \dots & 0\end{array}\right]\in {\mathbb{R}}^{n\times n},$$

contains second order terms in $\delta v_j,j=1,2,\dots ,n$. The diagonal entries of $\delta D_V$ are determined as in the case of the matrix $\delta D_U$.

4. Equations for the Perturbation Parameters

The elements of the perturbation parameter vectors x and y can be determined from equation (5). For this aim it is appropriate to transform this equation as follows. Taking into account that $AV=U\Sigma =U_1{\Sigma }_n$ and $U^{T}A=\Sigma V^T$, the equation is represented as

$$\delta U^T{U}_1{\Sigma }_n+\left[\begin{array}{c}{\Sigma }_n\\ 0_{(m-n)\times n}\end{array}\right]V^T\delta V+\delta F=\left[\begin{array}{c}{\tilde{\Sigma }}_n-{\Sigma }_n\\ 0_{(m-n)\times n}\end{array}\right]+{\Delta }_0.$$

(12)

where

$$\delta F:=\left[\begin{array}{c}\delta F_1\\ \delta F_2\end{array}\right]=\left[\begin{array}{cccc}\delta f_{11}& \delta f_{12}& \dots & \delta f_{1n}\\ \delta f_{21}& \delta f_{22}& \dots & \delta f_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ \delta f_{n1}& \delta f_{n2}& \dots & \delta f_{nn}\\ \delta f_{n+1,1}& \delta f_{n+1,2}& \dots & \delta f_{n+1,n}\\ ⋮& ⋮& \ddots & ⋮\\ \delta f_{m1}& \delta f_{m2}& \dots & \delta f_{mn}\end{array}\right].$$

According to (8), we have that

$$\delta U^T{U}_1=-U^T\delta U_1-\delta U^T\delta U_1.$$

(13)

Substituting in (12) the term $\delta U^T{U}_1$ by the sum in the right hand side of (13), we obtain

$$-U^T\delta U_1{\Sigma }_n+\left[\begin{array}{c}{\Sigma }_n\\ 0_{(m-n)\times n}\end{array}\right]V^T\delta V+\delta F=\left[\begin{array}{c}{\tilde{\Sigma }}_n-{\Sigma }_n\\ 0_{(m-n)\times n}\end{array}\right]+\Delta ,$$

(14)

where

$$\begin{array}{ccc}\hfill \Delta & =& {\Delta }_0+\delta U^T\delta U_1{\Sigma }_n\hfill \\ & =& \delta U^T\delta U_1{\Sigma }_n-\delta U^{T}A\delta V-U^T\delta A\delta V-\delta U^T\delta AV-\delta U^T\delta A\delta V\hfill \end{array}$$

(15)

contains higher order terms in the entries of $\delta A$, $\delta U$ and $\delta V$.

Replacing the matrices $U^T\delta U_1$ and $V^T\delta V$ by $\delta W_U$ and $\delta W_V$, respectively, (14) is rewritten as

$$-\delta W_U{\Sigma }_n+\left[\begin{array}{c}{\Sigma }_n\\ 0_{(m-n)\times n}\end{array}\right]\delta W_V+\delta F=\left[\begin{array}{c}{\tilde{\Sigma }}_n-{\Sigma }_n\\ 0_{(m-n)\times n}\end{array}\right]+\Delta$$

or

$$\begin{array}{ccc}\hfill -\delta Q_U{\Sigma }_n& +& \left[\begin{array}{c}{\Sigma }_n\\ 0_{(m-n)\times n}\end{array}\right]\delta Q_V+\delta F=\left[\begin{array}{c}{\tilde{\Sigma }}_n-{\Sigma }_n\\ 0_{(m-n)\times n}\end{array}\right]\hfill \\ & & +(\delta D_U-\delta N_U){\Sigma }_n-\Sigma (\delta D_V-\delta N_V)+\Delta .\hfill \end{array}$$

(16)

Note that the matrices $\delta D_U,\delta N_U,\delta D_V,\delta N_V$ and $\Delta$ contain only higher order terms in the entries of $\delta A,\delta U$ and $\delta V$.

The entries of the matrices $\delta Q_U$ and $\delta Q_V$ can be substituted by the corresponding elements $x_k,k=i+(j-1)m-\frac{j(j+1)}{2}$ and $y_{\ell },\ell =i+(j-1)n-\frac{j(j+1)}{2}$ of the vectors x and y as shown in the previous section. This leads to the representation of equation (16) as two matrix equations in respect to two groups of the entries of x,

(17)

and

$$\begin{array}{ccc}& & -\left[\begin{array}{ccccc}{\sigma }_1{x}_n& {\sigma }_2{x}_{m+n-2}& {\sigma }_3{x}_{2m+n-5}& \dots & {\sigma }_n{x}_{(n-1)(2m-n)/2+1}\\ {\sigma }_1{x}_{n+1}& {\sigma }_2{x}_{m+n-1}& {\sigma }_3{x}_{2m+n-4}& \dots & {\sigma }_n{x}_{(n-1)(2m-n)/2+2}\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ {\sigma }_1{x}_{m-1}& {\sigma }_2{x}_{2m-3}& {\sigma }_3{x}_{3m-6}& \dots & {\sigma }_n{x}_p\end{array}\right]\hfill \\ & & =-\left[\begin{array}{ccccc}\delta f_{n+1,1}& \delta f_{n+1,2}& \delta f_{n+1,3}& \dots & \delta f_{n+1,n}\\ \delta f_{n+2,1}& \delta f_{n+2,2}& \delta f_{n+2,3}& \dots & \delta f_{n+2,n}\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ \delta f_{m1}& \delta f_{m2}& \delta f_{m3}& \dots & \delta f_{mn}\end{array}\right]+{\Delta }_2,\hfill \end{array}$$

(18)

where

$$\begin{array}{ccc}\hfill {\Delta }_1& =& \delta U_1^T\delta U_1{\Sigma }_n-\delta U_1^{T}A\delta V-U_1^T\delta A\delta V-\delta U_1^T\delta AV-\delta U_1^T\delta A\delta V,\hfill \end{array}$$

(19)

$$\begin{array}{ccc}\hfill {\Delta }_2& =& \delta U_2^T\delta U_1{\Sigma }_n-\delta U_2^{T}A\delta V-U_2^T\delta A\delta V-\delta U_2^T\delta AV-\delta U_2^T\delta A\delta V,\hfill \end{array}$$

(20)

$$\Delta :=\left[\begin{array}{c}{\Delta }_1\\ {\Delta }_2\end{array}\right],{\Delta }_1\in {\mathbb{R}}^{n\times n},{\Delta }_2\in {\mathbb{R}}^{(m-n)\times n}.$$

We note that the estimation of ${\Delta }_2$ requires to know an estimate of $\delta U_2$ which is undetermined.

Equations (17) and (18) are the basic equation of the SVD perturbation analysis. They can be used to obtain asymptotic as well as global perturbation bounds on the elements of the vectors x and y.

Let us introduce the vectors

$$x_{\left(1\right)}={\Omega }_{1}x,and{x}_{\left(2\right)}={\Omega }_{2}x,$$

where

$$\begin{array}{ccc}\hfill {\Omega }_1& :=& \left[\mathrm{diag}({\omega }_1,{\omega }_2,\dots ,{\omega }_n)\right]\in {\mathbb{R}}^{q\times p},\hfill \\ \hfill {\omega }_i& :=& \left[I_{n-i},0_{(n-i)|\times (m-n)}\right]\in {\mathbb{R}}^{(n-i)\times (m-i)},i=1,2,\dots ,n,\hfill \\ & & {\Omega }_1{\Omega }_1^T=I_q,{\parallel {\Omega }_1\parallel }_2=1,\hfill \\ \hfill {\Omega }_2& :=& \left[\mathrm{diag}({\omega }_1,{\omega }_2,\dots ,{\omega }_n)\right]\in {\mathbb{R}}^{(m-n)n\times p},\hfill \\ \hfill {\omega }_i& :=& \left[0_{(m-n)\times (n-i)},I_{m-n}\right]\in {\mathbb{R}}^{(m-n)\times (m-i)},i=1,2,\dots ,n,\hfill \\ & & {\Omega }_2{\Omega }_2^T=I_{(m-n)n},{\parallel {\Omega }_2\parallel }_2=1.\hfill \end{array}$$

The vector $x_{\left(1\right)}$ contains the elements of the unknown vector x participating in (17), while $x_{\left(2\right)}$ contains the elements of x participating in (18). It is easy to prove that

$$x={\Omega }_1^T{x}_{\left(1\right)}+{\Omega }_2^T{x}_{\left(2\right)}.$$

Taking into account that

$$\begin{array}{ccc}\hfill \mathrm{Low}((\delta {D_U}_1-\delta {N_U}_1){\Sigma }_n-{\Sigma }_n(\delta D_V-\delta N_V))& =& 0,\hfill \\ \hfill \mathrm{Up}(\delta {D_U}_1{\Sigma }_n-{\Sigma }_n\delta D_V)& =& 0,\hfill \end{array}$$

the strictly lower part of (17) can be represented columnwise as the system of linear equations in respect to the unknown vectors $x_{\left(1\right)}$ and y,

$$-S_1{x}_{\left(1\right)}+S_{2}y=-f+\mathrm{vec}\left(\mathrm{Low}\left({\Delta }_1\right)\right),$$

(21)

where

$$\begin{array}{ccc}\hfill S_1& =& \mathrm{diag}(\underset{n-1}{\underbrace{{\sigma }_1,{\sigma }_1,\dots ,{\sigma }_1}},\underset{n-2}{\underbrace{{\sigma }_2,\dots ,{\sigma }_2}},\dots ,\underset{1}{\underbrace{{\sigma }_{n-1}}}),\hfill \end{array}$$

(22)

$$\begin{array}{ccc}\hfill S_2& =& \mathrm{diag}(\underset{n-1}{\underbrace{{\sigma }_2,{\sigma }_3,\dots ,{\sigma }_n}},\underset{n-2}{\underbrace{{\sigma }_3,\dots ,{\sigma }_n}},\dots ,\underset{1}{\underbrace{{\sigma }_n}}),\hfill \\ & & S_i\in {\mathbb{R}}^{q\times q}\hfill \end{array}$$

(23)

and

$$\begin{array}{ccc}\hfill f& =& \mathrm{vec}\left(\mathrm{Low}\left(\delta F_1\right)\right)={\Omega }_3\mathrm{vec}\left(\delta F_1\right)\in {\mathbb{R}}^q,\hfill \\ \hfill {\Omega }_3& :=& \left[\mathrm{diag}({\omega }_1,{\omega }_2,\dots ,{\omega }_n)\right]\in {\mathbb{R}}^{q\times n^2},\hfill \\ \hfill {\omega }_i& :=& \left[0_{(n-i)\times i},I_{n-i}\right]\in {\mathbb{R}}^{(n-i)\times n},i=1,2,\dots ,n,\hfill \\ & & {\Omega }_3{\Omega }_3^T=I_q,{\parallel {\Omega }_3\parallel }_2=1.\hfill \end{array}$$

Similarly, the strictly upper part of (17) is represented rowwise as the system of equations

$$S_2{x}_{\left(1\right)}-S_{1}y=-g-\mathrm{vec}\left({\left(\mathrm{Up}(\delta {N_U}_1{\Sigma }_n-{\Sigma }_n\delta N_V-{\Delta }_1)\right)}^T\right),$$

(24)

where

$$\begin{array}{ccc}\hfill g& =& \mathrm{vec}\left({\left(\mathrm{Up}\left(F_1\right)\right)}^T\right)={\Omega }_4\mathrm{vec}\left(\delta F_1\right)\in {\mathbb{R}}^q,\hfill \\ \hfill {\Omega }_4& :=& \left[\begin{array}{cc}{0}_{(n-1)\times n},\hfill & {\omega }_1,\\ 0_{(n-2)\times 2n},\hfill & {\omega }_2,\\ ⋮\hfill & ⋮\\ 0_{1\times (n-1)n},\hfill & {\omega }_{n-1}\end{array}\right]\in {\mathbb{R}}^{q\times n^2},\hfill \\ \hfill {\omega }_i& :=& \left[\begin{array}{cccc}[0_{1\times (i-1)},1,0_{1\times (n-i)}]& 0_{1\times n}& \dots & 0_{1\times n}\\ 0_{1\times n}& [0_{1\times (i-1)},1,0_{1\times (n-i)}]& \dots & 0_{1\times n}\\ ⋮& ⋮& \ddots & ⋮\\ 0_{1\times n}& 0_{1\times n}& \dots & [0_{1\times (i-1)},1,0_{1\times (n-i)}]\end{array}\right]\hfill \\ & & \in {\mathbb{R}}^{(n-i)\times (n-i).n},i=1,2,\dots ,n-1,\hfill \\ & & {\Omega }_4{\Omega }_4^T=I_q,{\parallel {\Omega }_4\parallel }_2=1.\hfill \end{array}$$

It should be noted that the operators $\mathrm{Low}$ and $\mathrm{Up}$ used in (21), (24), take only the entries of the strict lower and strict upper part, respectively, of the corresponding matrix which are placed by the operator $\mathrm{vec}$ column by column excluding the zeros above or under the diagonal, respectively. For instance, if $n=4$ the elements of the vectors f and g satisfy

In this way, the solution of (17) reduces to the solution of the two coupled equations (21) and (24) with diagonal matrices of size $q\times q$. The equation (18) can be solved independently yielding

$$x_{\left(2\right)}=\mathrm{vec}\left((\delta F_2-{\Delta }_2){\Sigma }_n^{-1}\right).$$

(25)

Note that the elements of $x_{\left(1\right)}$ depend on the elements of y and vice versa, while $x_{\left(2\right)}$ does not depend on y.

Thanks to the diagonal structure of the matrices $S_1,S_2$ and $\Sigma$, the equations (21) - (25) can be solved efficiently with high accuracy.

5. Asymptotic Bounds of the Perturbation Parameters

Equations (21) and (24) can be used to determine asymptotic approximations of the vectors $x_{\left(1\right)}$ and y. The exact solution of these equations satisfies

$$\begin{array}{ccc}\hfill x_{\left(1\right)}& =& S_{xf}(-f+\mathrm{vec}\left(\mathrm{Low}\left({\Delta }_1\right)\right)\hfill \\ & & +S_{xg}(-g-\mathrm{vec}\left({\left(\mathrm{Up}(\delta {N_U}_1{\Sigma }_n-{\Sigma }_n\delta N_V-{\Delta }_1)\right)}^T\right),\hfill \end{array}$$

(26)

$$\begin{array}{ccc}\hfill y& =& S_{yf}(-f+\mathrm{vec}\left(\mathrm{Low}\left({\Delta }_1\right)\right)\hfill \\ & & +S_{yg}(-g-\mathrm{vec}\left({\left(\mathrm{Up}(\delta {N_U}_1{\Sigma }_n-{\Sigma }_n\delta N_V-{\Delta }_1)\right)}^T\right),\hfill \end{array}$$

(27)

where taking into account that $S_1$ and $S_2$ commute, we have that

$$\begin{array}{ccc}\hfill S_{xf}& =& {(S_2-S_1{S}_2^{-1}{S}_1)}^{-1}{S}_1{S}_2^{-1}={(S_2^2-S_1^2)}^{-1}{S}_1\in {\mathbb{R}}^{q\times q},\hfill \\ \hfill S_{xg}& =& {(S_2-S_1{S}_2^{-1}{S}_1)}^{-1}={(S_2^2-S_1^2)}^{-1}{S}_2\in {\mathbb{R}}^{q\times q},\hfill \\ \hfill S_{yf}& =& S_{xg},\hfill \\ \hfill S_{yg}& =& S_{xf}.\hfill \end{array}$$

Exploiting these expressions, it is possible to show that

$$\begin{array}{ccc}\hfill Sxf& =& \mathrm{diag}\left(\underset{n-1}{\underbrace{\frac{{\sigma }_1}{{\sigma }_2^2-{\sigma }_1^2},\frac{{\sigma }_1}{{\sigma }_3^2-{\sigma }_1^2},\dots ,\frac{{\sigma }_1}{{\sigma }_n^2-{\sigma }_1^2}}},\underset{n-2}{\underbrace{\frac{{\sigma }_2}{{\sigma }_2^2-{\sigma }_1^2},\dots ,\frac{{\sigma }_2}{{\sigma }_n^2-{\sigma }_2^2}}},\dots ,\underset{1}{\underbrace{\frac{{\sigma }_{n-1}}{{\sigma }_n^2-{\sigma }_2^2}}}\right),\hfill \\ \hfill Sxg& =& \mathrm{diag}\left(\underset{n-1}{\underbrace{\frac{{\sigma }_2}{{\sigma }_2^2-{\sigma }_1^2},\frac{{\sigma }_3}{{\sigma }_3^2-{\sigma }_1^2},\dots ,\frac{{\sigma }_n}{{\sigma }_n^2-{\sigma }_1^2}}},\underset{n-2}{\underbrace{\frac{{\sigma }_3}{{\sigma }_3^2-{\sigma }_2^2},\dots ,\frac{{\sigma }_n}{{\sigma }_n^2-{\sigma }_2^2}}},\dots ,\underset{1}{\underbrace{\frac{{\sigma }_n}{{\sigma }_n^2-{\sigma }_{n-1}^2}}}\right).\hfill \end{array}$$

Let us consider the conditions for existence of a solution of equations (21), (24). These equations have an unique solution for $x_{\left(1\right)}$ and y, if and only if the matrix

$$\left[\begin{array}{cc}\hfill -S_1& \hfill S_2\\ \hfill S_2& \hfill -S_1\end{array}\right]$$

is nonsingular, or equivalently, the matrices $S_1$, $S_2$ and $S_2^2-S_1^2$ are nonsingular. The matrices $S_1$ and $S_2$ are nonsingular, since the matrix A has nonzero singular values. In turn, a condition for nonsingularity of the matrix $S_2^2-S_1^2$ can be found taking into account the structure of the matrices $Sxf$ and $Sxg$ shown above. Clearly, the denominators of the first $n-1$ diagonal entries of $Sxf$ and $Sxg$ will be different from zero, if ${\sigma }_1$ is distinct from ${\sigma }_2,{\sigma }_3,\dots ,{\sigma }_n$. Similarly, the denominators of the next group of $n-2$ diagonal entries will be different from zero if ${\sigma }_2$ is distinct from ${\sigma }_3,\dots ,{\sigma }_n$ and so on. Finally, ${\sigma }_{n-1}$ should be different from ${\sigma }_n$. Thus we come to the conclusion that the matrices $Sxf$ and $Sxg$ will exist and the equations (21), (24) will have an unique solution, if and only if the singular values of A are distinct. This conclusion should not come to surprise, since U is the matrix of the transformation of $A{A}^T$ to Schur (diagonal) form $U\Sigma {\Sigma }^T{U}^T$ and V is the matrix of the transformation of $A^{T}A$ to diagonal form $V{\Sigma }^T\Sigma V^T$. On the other hand, the perturbation problem for the Schur form is well posed only when the matrix eigenvalues (the diagonal elements of $\Sigma {\Sigma }^T$ or ${\Sigma }^T\Sigma$) are distinct.

Neglecting the higher order terms in (26), (27) and approximating each element of f and g by the perturbation norm ${\parallel \delta A\parallel }_2$, we obtain the linear estimates

$$\begin{array}{ccc}\hfill x_{\left(1\right)}^{lin}& =& \left(\right|S_{xf}|+|S_{xg}\left|\right)h,\hfill \end{array}$$

(28)

$$\begin{array}{ccc}\hfill y^{lin}& =& \left(\right|S_{yf}|+|S_{yg}\left|\right)h,\hfill \end{array}$$

(29)

where

$$h={\underset{q}{\underbrace{[1,1,\dots ,1]}}}^T\times {\parallel \delta A\parallel }_2.$$

Clearly, if the matrices $|S_{xf}|,|S_{xg}|,|S_{yf}|,|S_{yg}|$ have large diagonal elements, then the estimates of the perturbation parameters will be large. Using the expressions for $Sxf$ and $Sxg$, we may show that

$$\begin{array}{ccc}\hfill {\parallel Sxf\parallel }_2& =& \underset{i,j}{max}\frac{{\sigma }_i}{|{\sigma }_j^2-{\sigma }_i^2|},i=1,2,\dots ,n-1,j=i+1,i+2,\dots ,n,\hfill \\ \hfill {\parallel Sxg\parallel }_2& =& \underset{i,j}{max}\frac{{\sigma }_j}{|{\sigma }_j^2-{\sigma }_i^2|}.\hfill \end{array}$$

Note that the norms of $Sxf$ and $Sxg$ can be considered as condition numbers of the vectors $x_1$ and y with respect to the changes of $\delta A$.

An asymptotic estimate of the vector $x_{\left(2\right)}$ is obtained neglecting the higher order term ${\Delta }_2$ and approximating its elements according to (25) as

$$\begin{array}{ccc}\hfill x_{\left(2\right)}^{lin}& =& \mathrm{vec}\left(Z\right),\hfill \\ \hfill Z& =& {\parallel \delta A\parallel }_2\times \left[\begin{array}{cccc}1/{\sigma }_1& 1/{\sigma }_2& \dots & 1/{\sigma }_n\\ 1/{\sigma }_1& 1/{\sigma }_2& \dots & 1/{\sigma }_n\\ ⋮& ⋮& ⋮& ⋮\\ 1/{\sigma }_1& 1/{\sigma }_2& \dots & 1/{\sigma }_n\end{array}\right].\hfill \end{array}$$

(30)

Equation (30) shows that a group of n elements of ${x_{\left(2\right)}}^{lin}$ will be large, if the singular value participating in the corresponding column of Z is small. The presence of large elements in the vector x leads to large entries in $\delta W_U$ and consequently in the estimate of $\delta U$. This observation is in accordance with the well known fact that the sensitivity of a singular subspace is inversely proportional to the smallest singular value associated with this subspace.

As a result of solving the linear systems (28) - (30), we obtain an asymptotic approximation of the vector x as

$$x^{lin}={\Omega }_1^T{x}_{\left(1\right)}^{lin}+{\Omega }_2^T{x}_{\left(2\right)}^{lin}.$$

(31)

Example 1.

Consider the $6\times 4$ matrix

$$A=\left[\begin{array}{cccc}\hfill 3& \hfill 3& \hfill -3& \hfill -6\\ \hfill -3& \hfill -1& \hfill 1& \hfill 8\\ \hfill 3& \hfill 1& \hfill 0& \hfill -9\\ \hfill -3& \hfill 1& \hfill -2& \hfill 11.1\\ \hfill 6& \hfill 4& \hfill -6& \hfill -11.9\\ \hfill 3& \hfill 1& \hfill 1& \hfill -10.1\end{array}\right]$$

and assume that it is perturbed by

$$\begin{array}{ccc}\hfill \delta A& =& {10}^c\times A_0,\hfill \\ \hfill A_0& =& \left[\begin{array}{cccc}\hfill -5& \hfill 3& \hfill 1& \hfill -2\\ \hfill 7& \hfill -4& \hfill -9& \hfill -5\\ \hfill -3& \hfill 8& \hfill 5& \hfill 4\\ \hfill 2& \hfill -6& \hfill -3& \hfill -7\\ \hfill -5& \hfill 3& \hfill 1& \hfill 9\\ \hfill -4& \hfill -8& \hfill -2& \hfill 6\end{array}\right],\hfill \\ & & {\parallel \delta A\parallel }_2=2.081981\times {10}^c,\hfill \end{array}$$

where c is a varying parameter. (For convenience, the entries of the matrix $A_0$ are taken as integers.)

The singular value decompositions of the matrices A and $A+\delta A$ are computed by the functionsvdof MATLAB^®[7]. The singular values of A are

$$\begin{array}{ccc}\hfill {\sigma }_1& =& 25.460186918120350,\hfill \\ \hfill {\sigma }_2& =& 7.684342946021752,\hfill \\ \hfill {\sigma }_3& =& 1.248776186923002,\hfill \\ \hfill {\sigma }_4& =& 0.017709242587950.\hfill \end{array}$$

In the given case the matrices $S_1$ and $S_2$ in equations (22), (23), are

$$\begin{array}{ccc}\hfill S_1& =& \mathrm{diag}(25.4601869,25.4601869,25.4601869,\hfill \\ & & 7.6843429,7.6843429,1.2487762),\hfill \\ \hfill S_2& =& \mathrm{diag}(7.6843429,1.2487762,0.0177092,\hfill \\ & & 1.2487762,0.0177092,0.0177092),\hfill \end{array}$$

and the matrices participating in (26), (27) are

$$\begin{array}{ccc}\hfill S_{xf}& =& \mathrm{diag}(-0.0432135,-0.0393717,-0.0392770,\hfill \\ & & -0.1336647,-0.1301354,-0.8009451),\hfill \\ \hfill S_{xg}& =& \mathrm{diag}(-0.0130426,-0.0019311,-0.0000273,\hfill \\ & & -0.0217217,-0.0002999,-0.0113584)\hfill \\ \hfill S_{yf}& =& S_{xg},\hfill \\ \hfill S_{yg}& =& S_{xf}.\hfill \end{array}$$

The matrix

$$\left[\begin{array}{cc}\hfill -S_1& \hfill S_2\\ \hfill S_2& \hfill -S_1\end{array}\right],$$

which determines the solution for $x_{\left(1\right)}$ and y, has a condition number equal to $26.9234179$. The exact parameters $x_k$ and their linear approximates $x_k^{lin}$ computed by using (28) and (30), are shown to eight decimal digits for two perturbation sizes in Table 1 (the elements $x_k$ obtained from the equations (30) are highlighted in red boxes). The differences between the values of $x_k^{lin}$ and $x_k$ are due to the bounding of the elements of the vectors f and g by the value of ${\parallel \delta A\parallel }_2$ and taking the terms in (28) - (30) with positive signs. Both approximations are necessary to ensure that $x⪯x^{lin}$ for arbitrary small size perturbation.

Similarly, in Table 2 we show for the same perturbations of A the exact perturbation parameters $y_{\ell }$ and their linear approximations obtained from (29).

6. Bounding the Perturbations of the Singular Values

Equation (17) can also be used to determine linear and nonlinear estimates of the perturbations of the singular values. Considering the diagonal elements of this equation (highlighted in green boxes) and taking into account that $\mathrm{diag}\left(\delta {N_U}_1\right)=0_n,\mathrm{diag}\left(\delta N_V\right)=0_n$, we obtain

$$\delta {\Sigma }_n=\mathrm{diag}\left(\delta F_1\right)-\mathrm{diag}(\delta {D_U}_1{\Sigma }_n-{\Sigma }_n\delta D_V+{\Delta }_1)$$

or

$$\delta {\sigma }_i=\delta f_{ii}-{(\delta {D_U}_1{\Sigma }_n-{\Sigma }_n\delta D_V+{\Delta }_1)}_{ii},i=1,2,\dots ,n,$$

(32)

where ${(.)}_{ii}$ denotes the ith diagonal element of $(.)$. Neglecting the higher order terms, we determine the componentwise asymptotic estimate

$$\delta {\sigma }_i^{lin}=\left|\delta f_{ii}\right|,i=1,2,\dots ,n.$$

(33)

Bounding each diagonal element $|\delta f_{ii}|$ by ${\parallel \delta A\parallel }_2$, we find the normwise estimate of $\delta {\sigma }_i$,

$$\delta {\widehat{\sigma }}_i\le {\parallel \delta A\parallel }_2,$$

which is in accordance with the Weyl’s theorem (see (3)).

From (32) we also have that

$$\sqrt{\sum _{i=1}^n\delta {\sigma }_i^2}\le \parallel \delta F_1{\parallel }_F+\parallel \delta {D_U}_1{\Sigma }_n-{\Sigma }_n\delta D_V{\parallel }_F+{\parallel {\Delta }_1\parallel }_F.$$

In Table 3 we show the exact perturbations of the singular values of the matrix A of Example 1 along with the normwise bound $\delta {\widehat{\sigma }}_i={\parallel \delta A\parallel }_2$ and the asymptotic estimate $\delta {\sigma }_i^{lin}$ obtained from (33) under the assumption that $\delta A$ is known. The exact perturbations and their linear bounds are very close.

7. Asymptotic Bounds on the Perturbations of $U_1$ and V

Having componentwise estimates for the elements of x and y, it is possible to find bounds on the entries of the matrices $\delta U_1$ and $\delta V$. An asymptotic bound $\delta W_U^{lin}=\left|\delta Q_u\right|$ on the absolute value of the matrix $\delta W_U=U^T\delta U_1$ is given by

$$\delta W_U^{lin}=\left[\begin{array}{ccccc}0& |x_1|& |x_2|& \dots & |x_{n-1}|\\ |x_1|& 0& |x_m|& \dots & |x_{m+n-3}|\\ |x_2|& |x_m|& 0& ⋮& |x_{2m+n-6}|\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ |x_{n-1}|& |x_{m+n-3}|& |x_{m+2n-6}|& \dots & 0\\ |x_n|& |x_{m+n-2}|& |x_{2m+n-5}|& \dots & |x_{(n-1)(2m-n)/2+1}|\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ |x_{m-1}|& |x_{2m-3}|& |x_{3m-6}|& \dots & |x_p|\end{array}\right]\in {\mathbb{R}}^{m\times n}.$$

From (6), a linear approximation of the matrix $|\delta U_1|$ is determined as

$$|\delta U_1|⪯\delta U_1^{lin}=\left|U\right||U^T\delta U_1|=|U\left|\right|\delta W_U|.$$

(34)

The matrix $\delta U_1^{lin}$ gives bounds on the perturbations of the individual elements of the orthogonal transformation matrix U.

In a similar way we have that the linear approximation of the matrix $\delta W_V$ is given by

$$\delta W_V^{lin}=\left|\delta Q_V\right|=\left[\begin{array}{ccccc}0& |y_1|& |y_2|& \dots & |y_{n-1}|\\ |y_1|& 0& |y_n|& \dots & |y_{2n-3}|\\ |y_2|& |y_n|& 0& ⋮& |y_{3n-6}|\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ |y_{n-1}|& |y_{2n-3}|& |y_{3n-6}|& \dots & 0\end{array}\right]\in {\mathbb{R}}^{n\times n}.$$

From (7) we obtain that

$$\left|\delta V\right|⪯\delta V^{lin}=\left|V\right||\delta W_V|$$

(35)

Hence, the entries of the matrix $\delta V^{lin}$ give asymptotic estimates of the perturbations of the entries of V.

For the matrix A of Example 1, we obtain that the absolute values of the exact changes of the entries of the matrix $U_1$ for the perturbation $\delta A={10}^{-5}{A}_0$ satisfy

$$|\delta U_1|={10}^{-3}\times \left[\begin{array}{cccc}0.0020130& 0.0013238& 0.0330402& 0.2407168\\ 0.0009548& 0.0070093& 0.0356216& 0.9425346\\ 0.0006562& 0.0004107& 0.0576633& 0.1185295\\ 0.0000996& 0.0028839& 0.0020370& 0.1517173\\ 0.0004584& 0.0024637& 0.0316821& 0.3649600\\ 0.0004974& 0.0014734& 0.0630623& 0.1877213\end{array}\right]$$

and their asymptotic componentwise estimates found by using (34), are

$$\delta U_1^{lin}={10}^{-2}\times \left[\begin{array}{cccc}0.0018494& 0.0052303& 0.0215558& 0.9532765\\ 0.0012398& 0.0044262& 0.0221403& 1.3988805\\ 0.0014125& 0.0045112& 0.0219550& 0.5264197\\ 0.0017421& 0.0041695& 0.0159661& 0.5316375\\ 0.0015062& 0.0033947& 0.0117806& 0.6365127\\ 0.0016843& 0.0050286& 0.0180145& 0.9544088\end{array}\right].$$

Also, for the same $\delta A$ we have that

$$\left|\delta V\right|={10}^{-4}\times \left[\begin{array}{cccc}0.0188135& 0.0293793& 0.2621967& 0.0365476\\ 0.0348050& 0.0032323& 0.0065525& 0.2572914\\ 0.0153096& 0.0121733& 0.0836556& 0.0994065\\ 0.0036567& 0.0106224& 0.0446413& 0.0002219\end{array}\right]$$

and, according to (35),

$$\delta V^{lin}={10}^{-3}\times \left[\begin{array}{cccc}0.0110821& 0.0341290& 0.1616097& 0.0370981\\ 0.0130707& 0.0300485& 0.0179106& 0.1608764\\ 0.0161344& 0.0250801& 0.0803532& 0.1020436\\ 0.0056685& 0.0191880& 0.0672176& 0.0164314\end{array}\right].$$

It is seen that the magnitude of the entries of $\delta U^{lin}$ and $\delta V^{lin}$ reflect correctly the magnitude of the corresponding entries of $\left|\delta U\right|$ and $\left|\delta V\right|$, respectively. Note that the perturbations of the columns of U and V tend to increase with increasing of the column number.

8. Sensitivity of Singular Subspaces

The sensitivity of a left

$${\mathcal{U}}_r=\mathrm{span}(u_1,u_2,\dots ,u_r),r\le min\{m-1,n\}$$

or right

$${\mathcal{V}}_r=\mathrm{span}(v_1,v_2,\dots ,v_r),r\le n-1$$

singular subspace of dimension r is measured by the angles between the corresponding unperturbed and perturbed subspaces [14, Ch. V].

Let the unperturbed left singular subspace corresponding to the first r singular values is denoted by ${\mathcal{U}}_r$ and its perturbed counterpart as ${\tilde{\mathcal{U}}}_r$ and let $U_{\left(r\right)}$ and ${\tilde{U}}_{\left(r\right)}$ be the orthonormal bases for ${\mathcal{U}}_r$ and ${\tilde{\mathcal{U}}}_r$, respectively. Then the maximum angle between $\mathcal{U}$ and $\tilde{\mathcal{U}}$ is determined from

$$cos(\Phi max({\mathcal{U}}_r,{\tilde{\mathcal{U}}}_r)={\sigma }_{min}\left(U_{\left(r\right)}^T{\tilde{U}}_{\left(r\right)}\right).$$

(36)

The expression (36) has the disadvantage that if $\Phi max$ is small, then $cos\left(\Phi max\right)\approx 1$ and the angle $\Phi max$ is not well determined. To avoid this difficulty, instead of

$cos(\Phi max({\mathcal{U}}_r,{\tilde{\mathcal{U}}}_r))$ it is preferable to work with $sin(\Phi max({\mathcal{U}}_r,{\tilde{\mathcal{U}}}_r))$. It is possible to show that [2]

$$sin(\Phi max({\mathcal{U}}_r,{\tilde{\mathcal{U}}}_r))={\sigma }_{max}\left({U_{\left(r\right)}^{\perp }}^T{\tilde{U}}_{\left(r\right)}\right)$$

(37)

where $U_{\left(r\right)}^{\perp }$ is the orthogonal complement of $U_{\left(r\right)}$, ${U_{\left(r\right)}^{\perp }}^T{U}_{\left(r\right)}=0$. Since

$${\tilde{U}}_{\left(r\right)}=U_{\left(r\right)}+\delta U_{\left(r\right)},$$

we have that

$$sin\left(\Phi max({\mathcal{U}}_r,{\tilde{\mathcal{U}}}_r)\right)={\sigma }_{max}\left({U_{\left(r\right)}^{\perp }}^T\delta U_{\left(r\right)}\right).$$

(38)

Equation (38) shows that the sensitivity of the left singular subspace ${\mathcal{U}}_r$ is connected to the values of the perturbation parameters $x_k=u_i^T\delta u_j,k=i+(j-1)m-\frac{j(j+1)}{2},i>r,j=1,2,\dots ,r$. In particular, for $r=1$ the sensitivity of the first column of U (the left singular vector, corresponding to ${\sigma }_1$) is determined as

$$sin(\Phi max({\mathcal{U}}_1,{\tilde{\mathcal{U}}}_1))={\sigma }_{max}\left(\delta {W_U}_{2:m,1}\right),$$

for $r=2$ one has

$$sin\left(\Phi max({\mathcal{U}}_2,{\tilde{\mathcal{U}}}_2)\right)={\sigma }_{max}\left(\delta {W_U}_{3:m,1:2}\right)$$

and so on, see Figure 1, where the matrices ${U_{\left(r\right)}^{\perp }}^T\delta U_{\left(r\right)}=\delta {W_U}_{r+1:m,1:r}$ for different values of r are highlighted in boxes and (*) denotes entries which are not used.

In a similar way, utilizing the matrix $\delta W_V$, it is possible to find the sine of the maximum angle between the perturbed ${\tilde{\mathcal{V}}}_r$ and unperturbed right singular subspace ${\mathcal{V}}_r$,

$$sin\left(\Theta max({\mathcal{V}}_r,\widetilde{{\mathcal{V}}_r})\right)={\sigma }_{max}\left(\delta {W_V}_{r+1:n,1:r}\right)$$

(see Figure 2). Hence, if the perturbation parameters are determined, it is possible to find at once sensitivity estimates for the nested singular subspaces

$$\begin{array}{ccc}\hfill {\mathcal{U}}_1=\mathrm{span}\left(U_{\left(1\right)}\right)& =& \mathrm{span}\left(u_1\right),\hfill \\ \hfill {\mathcal{U}}_2=\mathrm{span}\left(U_{\left(2\right)}\right)& =& \mathrm{span}(u_1,u_2),\hfill \\ \hfill ⋮& & ⋮\hfill \\ \hfill {\mathcal{U}}_r=\mathrm{span}\left(U_{\left(r\right)}\right)& =& \mathrm{span}(u_1,u_2,\dots ,u_r),r=min\{m-1,n\}\hfill \end{array}$$

and

$$\begin{array}{ccc}\hfill {\mathcal{V}}_1=\mathrm{span}\left(V_{\left(1\right)}\right)& =& \mathrm{span}\left(v_1\right),\hfill \\ \hfill {\mathcal{V}}_2=\mathrm{span}\left(V_{\left(2\right)}\right)& =& \mathrm{span}(v_1,v_2),\hfill \\ \hfill ⋮& & ⋮\hfill \\ \hfill {\mathcal{V}}_{n-1}=\mathrm{span}\left(V_{\left(n\right)}\right)& =& \mathrm{span}(v_1,v_2,\dots ,v_{n-1}).\hfill \end{array}$$

Specifically, as

$$\delta W_U=\left[\begin{array}{ccccc}*& *& *& \dots & *\\ x_1& *& *& \dots & *\\ x_2& x_m& *& \dots & *\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ x_{n-1}& x_{m+n-3}& x_{2m+n-6}& \dots & *\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ x_{m-1}& x_{2m-3}& x_{3m-6}& \dots & x_p\end{array}\right]\in {\mathbb{R}}^{m\times n},$$

we have that the exact maximum angle between the perturbed and unperturbed left singular subspace of dimension r is given by

$$\Phi max({\mathcal{U}}_r,{\tilde{\mathcal{U}}}_r)=arcsin\left({\sigma }_{max}\left(\delta {W_U}_{r+1:m,1:r}\right)\right).$$

(39)

Similarly, the maximum angle between the perturbed and unperturbed right singular subspace of dimension r is

$$\Theta max({\mathcal{V}}_r,{\tilde{\mathcal{V}}}_r)=arcsin\left({\sigma }_{max}\left(\delta {W_V}_{r+1:m,1:r}\right)\right).$$

(40)

To find an asymptotic estimate of $\Phi max({\mathcal{U}}_r,{\tilde{\mathcal{U}}}_r)$, in the expression for the matrix $\delta W_U$ the elements $x_k$ are replaced by their linear approximations (31). Let

$$\delta W_U^{lin}=\left[\begin{array}{ccccc}*& *& *& \dots & *\\ x_1^{lin}& *& *& \dots & *\\ x_2^{lin}& x_m^{lin}& *& \dots & *\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ x_{n-1}^{lin}& x_{m+n-3}^{lin}& x_{2m+n-6}^{lin}& \dots & *\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ x_{m-1}^{lin}& x_{2m-3}^{lin}& x_{3m-6}^{lin}& \dots & x_p^{lin}\end{array}\right]\in {\mathbb{R}}^{m\times n}$$

Then the following asymptotic estimate holds,

$${\Phi max}^{lin}({\mathcal{U}}_r,{\tilde{\mathcal{U}}}_r)=arcsin(\parallel \delta {W_U}_{r+1:m,1:r}^{lin}{\parallel }_2)$$

(41)

In particular, for the sensitivity of the range $\mathcal{R}\left(A\right)$ of A we obtain that

$$sin\left({\Phi max}^{lin}({\mathcal{U}}_n,{\tilde{\mathcal{U}}}_n)\right)={\parallel \delta {W_U}_{n+1:m,1:n}^{lin}\parallel }_2.$$

Similarly, for the angles between the unperturbed and perturbed right singular subspaces we obtain the linear estimates

$${\Theta max}^{lin}({\mathcal{V}}_r,{\tilde{\mathcal{V}}}_r)=arcsin(\parallel \delta {W_V}_{r+1:n,1:r}^{lin}\parallel ).$$

(42)

We note that the use of separate x and y parameters decouples the SVD perturbation problem and makes it possible to determine the sensitivity estimates of the left and right singular subspaces independently. This is important in cases when the left or right subspace in a pair of singular subspaces is much more sensitive than its counterpart.

Consider the same perturbed matrix A as in Example 1. Computing the matrices $\delta W_U^{lin}$ and $\delta W_V^{lin}$, it is possible to estimate at once the sensitivity of all four pairs of singular subspaces of dimensions $r=1,2,3,4$ corresponding to the chosen ordering of the singular values. In Table 4 we show the actual values of the left and right singular subspaces sensitivity and the computed asymptotic estimates (41) of this sensitivity. To determine the sensitivity of other singular subspaces, it is necessary to reorder the singular values in the initial decomposition so that the desired subspace appears in the set of nested singular subspaces. Also, the computations related to the determining of the linear estimates should be done again.

9. Global Perturbation Bounds

Since analytical expressions for the global perturbation bounds of the singular value decompositions are not known up to this moment, we present an iterative procedure for finding estimates of these bounds based on the asymptotic analysis presented above. This procedure is similar to the corresponding iterative schemes proposed in [8,9], but is more complicated since the determining of bounds on the parameter vectors x and y must be done simultaneously due to fact that the equations for these parameters are coupled.

9.1. Perturbation Bounds of the Entries of $U_2$

The main difficulty in determining global bounds of x and y is to find an appropriate approximation of the high order term ${\Delta }_2$ in (18). As is seen from (20), the determining of such estimate requires to know the perturbation $\delta U_2$ which is not well determined since it contains the columns of the matrix $\delta U_2={\tilde{U}}_2-U_2$. This perturbation satisfies the equations

$$\begin{array}{ccc}\hfill {(U_1+\delta U_1)}^T(U_2+\delta U_2)& =& 0,\hfill \end{array}$$

(43)

$$\begin{array}{ccc}\hfill {(U_2+\delta U_2)}^T(U_2+\delta U_2)& =& I_{m-n}\hfill \end{array}$$

(44)

which follow from the orthogonality of the matrix

$$\tilde{U}=[U_1+\delta U_1,U_2+\delta U_2].$$

An estimate of $\delta U_2$ can be found based on a suitable approximation of

$$X=U^T\delta U_2:=\left[\begin{array}{c}{X}_1\\ X_2\end{array}\right].$$

As shown in [9], a first order approximation of the matrix X can be determined using the estimates

$$\begin{array}{ccc}\hfill X_1^{appr}& =& -{(I_n+\delta W_1^T)}^{-1}\delta W_2^T\in {\mathbb{R}}^{n\times (m-n)},\hfill \end{array}$$

(45)

$$\begin{array}{ccc}\hfill X_2^{appr}& =& -{X_1^{appr}}^T{X}_1^{appr}/2\in {\mathbb{R}}^{(m-n)\times (m-n)},\hfill \end{array}$$

(46)

where $\delta W_1=U_1^T\delta U_1$, $\delta W_2=U_2^T\delta U_1$ and for sufficiently small perturbation $\delta U_1$ the matrix $I_n+\delta W_1^T$ is nonsingular. (Note that $\delta W_U={\left[\delta W_1^T\delta W_2^T\right]}^T$ is already estimated.) Thus, we have that

$$\delta U_2^{appr}=U{X}^{appr}.$$

(47)

9.2. Iterative Procedure for Finding Global Bounds of x and y

Global componentwise perturbation bounds of the matrices U and V can be found using nonlinear estimates of the matrices $\delta W_U$ and $\delta W_V$, determined by (10) and (11), respectively. Such estimates are found correcting the linear estimates of the perturbation parameters $x_k=u_i^T\delta u_j$ and $y_{\ell }=v_i^T\delta v_j$ on each iteration step in a way similar to the one presented in [8] and [9].

Consider the case of estimating the matrix $\delta W_U$. It is convenient to substitute the terms containing the perturbations $\delta u_j$ in (10) by the quantities

$$\delta {w_u}_j=U^T\delta u_j,j=1,2,\dots ,n,$$

which have the same magnitude as $\delta u_j$. Since

$$\delta u_i^T\delta u_j=\delta u_i^{T}U{U}^T\delta u_j=\delta {w_u}_i^T\delta {w_u}_j,$$

the absolute value of the matrix $\delta W_U\in {\mathbb{R}}^{m\times n}$ (10) can be bounded as

$$\begin{array}{ccc}\hfill |\delta W_U^{nonl}|& =& |U^T\delta U_1|=[|\delta {w_u}_1|,\delta {w_u}_2|,\dots ,|\delta {w_u}_n\left|\right]\hfill \\ & ⪯& |\delta Q_U|+|\delta D_U|+|\delta N_U|,\hfill \end{array}$$

(48)

where

$$|\delta Q_U|:=\left[\begin{array}{c}\delta {Q_U}_1\\ \delta {Q_U}_2\end{array}\right]=\left[\begin{array}{ccccc}0& |x_1|& |x_2|& \dots & |x_{n-1}|\\ |x_1|& 0& |x_m|& \dots & |x_{m+n-3}|\\ |x_2|& |x_m|& 0& \dots & |x_{2m+n-6}|\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ |x_{n-1}|& |x_{m+n-3}|& |x_{2m+n-6}|& \dots & 0\\ |x_n|& |x_{m+n-2}|& |x_{2m+n-5}|& \dots & |x_{(n-1)(2m-n)/2+1}|\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ |x_{m-1}|& |x_{2m-3}|& |x_{3m-6}|& \dots & |x_p|\end{array}\right]\in {\mathbb{R}}^{m\times n},$$

$$|\delta D_U|:=\left[\begin{array}{c}\delta {D_U}_1\\ \delta {D_U}_2\end{array}\right]=\left[\begin{array}{ccccc}|{d_u}_{11}|& 0& 0& \dots & 0\\ 0& |{d_u}_{22}|& 0& \dots & 0\\ 0& 0& |{d_u}_{33}|& \dots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \dots & |{d_u}_{nn}|\\ 0& 0& 0& \dots & 0\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ 0& 0& 0& \dots & 0\end{array}\right]\in {\mathbb{R}}^{m\times n},$$

$$|\delta N_U|:=\left[\begin{array}{c}\delta {N_U}_1\\ \delta {N_U}_2\end{array}\right]=\left[\begin{array}{ccccc}0& |\delta {w_u}_1^T\left|\right|\delta {w_u}_2|& |\delta {w_u}_1^T\left|\right|\delta {w_u}_3|& \dots & |\delta {w_u}_1^T\left|\right|\delta {w_u}_n|\\ 0& 0& |\delta {w_u}_2^T\left|\right|\delta {w_u}_3|& \dots & |\delta {w_u}_2^T\left|\right|\delta {w_u}_n|\\ 0& 0& 0& \dots & |\delta {w_u}_3^T\left|\right|\delta {w_u}_n|\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \dots & |\delta {w_u}_{n-1}^T\left|\right|\delta {w_u}_n|\\ 0& 0& 0& \dots & 0\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ 0& 0& 0& \dots & 0\end{array}\right]\in {\mathbb{R}}^{m\times n}.$$

Since the unknown column estimates $|\delta w_j|$ participate in both sides of (48), it is possible to obtain them as follows. The first column of $\delta W_U$ is determined from

$$|\delta {w_u}_1|=|\delta {q_u}_1|+|\delta {d_u}_1|,$$

where $|\delta {q_u}_1|$, $|\delta d_1|$ are the first columns of $|\delta Q_U|$, $|\delta D_U|$, respectively. Then the next column estimates $|\delta w_j|,j=2,3,\dots n$ can be determined recursively from

$$|\delta {w_u}_j|=|\delta {q_u}_j|+|\delta {d_u}_j|+\left[\begin{array}{c}|\delta {w_u}_1^T\left|\right|\delta {w_u}_j|\\ |\delta {w_u}_2^T\left|\right|\delta {w_u}_j|\\ ⋮\\ |\delta {w_u}_{j-1}^T\left|\right|\delta {w_u}_j|\\ 0\\ ⋮\\ 0\end{array}\right],$$

which is equivalent to solving the linear system

$$|{S_U}_j\left|\right|\delta {w_u}_j|=|\delta {q_u}_j|+|\delta {d_u}_j|,$$

where

$$\left|{S_U}_j\right|=\left[\begin{array}{c}{e}_1^T-\left|\delta {w_u}_1^T\right|\\ e_2^T-\left|\delta {w_u}_2^T\right|\\ ⋮\\ e_{j-1}^T-\left|\delta {w_u}_{j-1}^T\right|\\ e_j^T\\ ⋮\\ e_n^T\end{array}\right]\in {\mathbb{R}}^{m\times m}$$

and $e_j$ is the jth column of $I_m$. The matrix $|{S_U}_j|$ is upper triangular with unit diagonal and if $|\delta {w_u}_i|,i=1,2,\dots ,j-1$ have small norms, then the matrix $|{S_U}_j|$ is diagonally dominant. Hence, it is very well conditioned with condition number close to 1.

As a result we obtain that

$$|\delta {w_u}_j|⪯|{S_U}_j{|}^{-1}\left(\right|\delta {q_u}_j|+|\delta {d_u}_j\left|\right)$$

(49)

which produces the jth column of $|\delta W_U^{nonl}|$.

A similar recursive procedure can be used to determine the quantities $\left|\delta {w_v}_j\right|=\left|v_i^T\delta v_j\right|$. In this case for each j it is necessary to solve the nth order linear system

$$\left|{S_V}_j\right|\left|\delta {w_v}_j\right|=\left|\delta {q_v}_j\right|+\left|\delta {d_v}_j\right|.$$

The estimates of $|\delta {w_u}_j|,j=1,2,\dots ,m$ and $|\delta {w_v}_j|,j=1,2,\dots ,n$ thus obtained, are used to bound the absolute values of the nonlinear elements ${\Delta }_1$ and ${\Delta }_2$ given in (19) and (20), respectively. Utilizing the approximation of $U^T\delta U_2$, it is possible to find an approximation of the matrix $U^T\delta U$ as

$$Z=[\delta W_U^{nonl},\delta U_2^{appr}]:=[Z_1,Z_2],Z_1\in {\mathbb{R}}^{m\times n},Z_2\in {\mathbb{R}}^{m\times (m-n)}$$

where $\delta U_2^{appr}=U{X}^{appr}$,

$$X^{appr}=\left[\begin{array}{c}{X}_1^{appr}\\ X_2^{appr}\end{array}\right]$$

and $X_1^{appr},X_2^{appr}$ are given by (45), (46). Then the elements of ${\Delta }_1$, ${\Delta }_2$ are bounded according to (19) and (20) as

$$\begin{array}{ccc}\hfill |{{\Delta }_1}_{ij}^{nonl}|& \le & {\sigma }_j\left|{Z_1}_i^T{Z_1}_j\right|+|{Z_1}_i^T\Sigma Y_j|+\hfill \\ & & (\parallel {Z_1}_i{\parallel }_2+\parallel Y_j{\parallel }_2+\parallel {Z_1}_i{\parallel }_2\parallel Y_j{\parallel }_2{)\parallel \delta A\parallel }_2,\hfill \\ & & i=1,2,\dots ,n,j=1,2,\dots ,n,\hfill \end{array}$$

(50)

$$\begin{array}{ccc}\hfill |{{\Delta }_2}_{ij}^{nonl}|& \le & {\sigma }_j|{Z_2}_i^T{Z_1}_j|+|{Z_2}_i^T\Sigma Y_j|+\hfill \\ & & (\parallel {Z_2}_i{\parallel }_2+\parallel Y_j{\parallel }_2+\parallel {Z_2}_i{\parallel }_2\parallel Y_j{\parallel }_2{)\parallel \delta A\parallel }_2.\hfill \\ & & i=1,2,\dots ,m-n,j=1,2,\dots ,n.\hfill \end{array}$$

(51)

Utilizing (26) and (27), the nonlinear corrections of the vectors $x_{\left(1\right)}$ and y can be determined from

$$\begin{array}{ccc}\hfill \delta x_{\left(1\right)}& =& |S_{xf}\left|\mathrm{vec}\right(\mathrm{Low}\left(\right|{\Delta }_1\left|\right)+|S_{xg}\left|\mathrm{vec}\right(\left(\mathrm{Up}\right(|\delta {N_U}_1|{\Sigma }_n\hfill \end{array}$$

$$\begin{array}{ccc}& & +{\Sigma }_n|\delta N_V|+|{\Delta }_1^{nonl}{\left|\right))}^T),\hfill \\ \hfill \delta y& =& |S_{yf}\left|\mathrm{vec}\right(\mathrm{Low}\left(\right|{\Delta }_1\left|\right)+|S_{yg}\left|\mathrm{vec}\right(\left(\mathrm{Up}\right(|\delta {N_U}_1|{\Sigma }_n\hfill \end{array}$$

(52)

$$\begin{array}{ccc}& & +{\Sigma }_n|\delta N_V|+|{\Delta }_1{|}^{nonl}{\left)\right)}^T),\hfill \end{array}$$

(53)

where $|\delta {N_U}_1|$ is estimated by using the corresponding expression (48) and $|\delta N_V|$ - by a similar expression.

The nonlinear correction of $x_{\left(2\right)}$ is found from

$$\begin{array}{ccc}\hfill \delta x_{\left(2\right)}& =& \mathrm{vec}\left(Z\right),\hfill \\ \hfill Z& =& \parallel {\Delta }_2^{nonl}{\parallel }_2\times \left[\begin{array}{cccc}1/{\sigma }_1& 1/{\sigma }_2& \dots & 1/{\sigma }_n\\ 1/{\sigma }_1& 1/{\sigma }_2& \dots & 1/{\sigma }_n\\ ⋮& ⋮& ⋮& ⋮\\ 1/{\sigma }_1& 1/{\sigma }_2& \dots & 1/{\sigma }_n\end{array}\right].\hfill \end{array}$$

(54)

and the total correction vector is determined from

$$\delta x={\Omega }_1^T\delta x_{\left(1\right)}+{\Omega }_2^T\delta x_{\left(2\right)}.$$

(55)

Now, the nonlinear estimates of the vectors x and y are found from

$$\begin{array}{ccc}\hfill x^{nonl}& =& x^{lin}+\delta x,\hfill \end{array}$$

(56)

$$\begin{array}{ccc}\hfill y^{nonl}& =& y^{lin}+\delta y.\hfill \end{array}$$

(57)

In this way we obtain an iterative scheme for finding simultaneously nonlinear estimates of the coupled perturbation parameter vectors x and y involving the equations (48) - (51), (52) - (57). In the numerical experiments presented below, the initial conditions are chosen as $x_0^{nonl}=eps{[1,1,\dots ,1]}^T$ and $y_0^{nonl}=eps{[1,1,\dots ,1]}^T$, where $eps$ is the MATLAB^®function eps, $eps=2^{-52}$. The stopping criteria for x- and y-iterations are taken as

$$\begin{array}{ccc}\hfill er{r}_x& =& \parallel x_s^{nonl}-x_{s-1}^{nonl}{\parallel }_2/{\parallel x_{s-1}^{nonl}\parallel }_2<tol,\hfill \\ \hfill er{r}_y& =& \parallel y_s^{nonl}-y_{s-1}^{nonl}{\parallel }_2/{\parallel y_{s-1}^{nonl}\parallel }_2<tol,\hfill \end{array}$$

where $tol=10eps$. The scheme converges for perturbations $\delta A$ of restricted size. It is possible that y converges while x does not converge.

The nonlinear estimate of the higher term ${\Delta }_1$ can be used to obtain nonlinear corrections of the singular value perturbations. Based on (32), a nonlinear correction of each singular value can be determined as

$$\delta {\sigma }_i^{corr}=\left(\right|\delta {D_U}_1|{\Sigma }_n+{\Sigma }_n|\delta D_V|+|{\Delta }_1{\left|\right)}_{ii},$$

(58)

so that the corresponding singular value perturbation is estimated as

$$\delta {\sigma }_i^{nonl}={\sigma }_i^{lin}+\delta {\sigma }_i^{corr}.$$

Note that ${\sigma }_i^{lin}=\left|\delta f_{ii}\right|$ is known only when the entries of the perturbation $\delta A$ are known and usually this is not fulfilled in practice. Nevertheless, the nonlinear correction (58) can be useful in estimating the sensitivity of a given singular value.

In Table 5, we present the number of iterations necessary to find the global bound ${\Delta }^{nonl}$ for the problem considered in Example 1 with perturbations $\delta A={10}^c\times A_0,c=-10,-9,\dots ,-3$. In the last two columns of the table we give the norm of the exact higher order term $\Delta$ and its approximation $\parallel {\Delta }^{nonl}{\parallel }_2$ computed according to (50), (51) (the approximation is given for the last iteration). In particular, for the perturbation $\delta A={10}^{-5}{A}_0$, the exact higher order term $\Delta$, found using (19) and (20), is

$$\begin{array}{ccc}\hfill |\Delta |& =& {10}^{-7}\times \left[\begin{array}{cccc}0.0045563& 0.0006839& 0.0146195& 0.0257942\\ 0.0012737& 0.0008692& 0.0080543& 0.0003708\\ 0.0000993& 0.0015257& 0.0086067& 0.0085598\\ 0.0019107& 0.0007767& 0.0258012& 0.0170781\\ 0.0032494& 0.0006701& 0.0423789& 0.2046100\\ 0.0064032& 0.0004567& 0.0384765& 0.0071094\end{array}\right],\hfill \\ \hfill {\parallel \Delta \parallel }_2& =& 2.12746\times {10}^{-8}.\hfill \end{array}$$

Implementing the described iterative procedure, after 10 iterations we obtain the nonlinear bound

$$\begin{array}{ccc}\hfill {\Delta }^{nonl}& =& {10}^{-5}\times \left[\begin{array}{cccc}0.0019716& 0.0021343& 0.0049111& 0.0047911\\ 0.0041273& 0.0053422& 0.0069174& 0.0069309\\ 0.0189857& 0.0186412& 0.0222135& 0.0179503\\ 0.8671798& 0.8680878& 0.8664829& 0.8607932\\ 0.5146400& 0.5162827& 0.5207565& 0.2618101\\ 0.5146400& 0.5162827& 0.5207565& 0.2618101\end{array}\right],\hfill \\ \hfill \parallel {\Delta }^{nonl}{\parallel }_2& =& 2.16338\times {10}^{-5}\hfill \end{array}$$

computed according to (50) and (51) on the base of the nonlinear bound $x^{nonl}$.

The global bounds $\parallel x^{nonl}{\parallel }_2$ and $\parallel y^{nonl}{\parallel }_2$, found for different perturbations along with the values of ${\parallel x\parallel }_2$ and ${\parallel y\parallel }_2$, are shown in Table 6. The results confirm that the global estimates of x and y are close to the corresponding asymptotic estimates.

In Figure 2 and Figure 3 we show the convergence of the relative errors

$$\parallel x_s^{nonl}-x_{s-1}^{nonl}\parallel /\parallel x_{s-1}^{nonl}\parallel$$

and

$$\parallel y_s^{nonl}-y_{s-1}^{nonl}\parallel /\parallel y_{s-1}^{nonl}\parallel ,$$

respectively, at step s of the iterative process for different perturbations $\delta A={10}^c\times A_0$. For the given example the iteration behaviours of x and y are close.

As it is seen from the figures, with the increasing of the perturbation size the convergence worsens and for $c=-3$$(\parallel \delta A{\parallel }_2=2.08198\times {10}^{-2})$ the iteration diverges. This demonstrates the restricted usefulness of the nonlinear estimates which are valid only for limited perturbation magnitudes.

In Table 7 we give normwise perturbation bounds of the singular values along with the actual singular value perturbations and their global bounds found for two perturbations of A under the assumption that the linear bounds of all singular values are known. As it can be seen from the table, the nonlinear estimates of the singular values are very tight.

9.3. Global Perturbation Bounds of $\delta U_1$ and $\delta V$

Having nonlinear bounds of x, y, $|\delta W_U|$ and $|\delta W_V|$, we may find nonlinear bounds on the perturbations of the entries of $U_1$ and V according to the relationships

$$\begin{array}{c}\hfill \delta U_1^{nonl}=\left|U\right||\delta W_U^{nonl}|,\end{array}$$

(59)

$$\begin{array}{c}\hfill \delta V^{nonl}=\left|V\right||\delta W_V^{nonl}|.\end{array}$$

(60)

For the perturbations of the orthogonal matrices of Example 1 we obtain the nonlinear componentwise bounds

$$\delta U_1^{nonl}={10}^{-2}\times \left[\begin{array}{cccc}0.0018632& 0.0052760& 0.0300043& 0.9595752\\ 0.0012438& 0.0044395& 0.0270841& 1.4007065\\ 0.0014395& 0.0046009& 0.0279076& 0.5379336\\ 0.0017550& 0.0042123& 0.0215531& 0.5375519\\ 0.0015075& 0.0033988& 0.0146480& 0.6375027\\ 0.0016913& 0.0050519& 0.0262884& 0.9579537\end{array}\right]$$

and

$$\delta V^{nonl}={10}^{-3}\times \left[\begin{array}{cccc}0.0110870& 0.0341508& 0.1619665& 0.0371843\\ 0.0130758& 0.0300693& 0.0179405& 0.1612236\\ 0.0161397& 0.0250977& 0.0805035& 0.1022483\\ 0.0056705& 0.0191966& 0.0673461& 0.0164515\end{array}\right].$$

These bounds are close to the obtained in sect. Section 7 linear estimates $\delta U_1^{lin}$ and $\delta V^{lin}$, respectively.

Based on (39), (40), global estimates of the maximum angles between the unperturbed and perturbed singular subspaces of dimension r can be obtained using the nonlinear bounds $\delta W_U^{nonl}$ and $\delta W_V^{nonl}$ of the matrices $\delta W_U$ and $\delta W_V$, respectively. For the pair of left and right singular subspaces we obtain that

$$\begin{array}{ccc}\hfill {\Phi max}^{nonl}({\mathcal{U}}_r,{\tilde{\mathcal{U}}}_r)& =& arcsin(\parallel {\delta W_U}_{r+1:m,1:r}^{nonl}{\parallel }_2),\hfill \\ & & r\le min\{m-1,n\},\hfill \end{array}$$

(61)

$$\begin{array}{ccc}\hfill {\Theta max}^{nonl}({\mathcal{V}}_r,{\tilde{\mathcal{V}}}_r)& =& arcsin(\parallel {\delta W_V}_{r+1:n,1:r}^{nonl}{\parallel }_2),\hfill \\ & & r\le n-1.\hfill \end{array}$$

(62)

In Table 8 we give the exact angles between the perturbed and unperturbed left and singular subspaces of different dimensions and their nonlinear bounds computed using (61) and (62) for the matrix A from Example 1 and two perturbations $\delta A={10}^c\times A_0,c=-10,-5$. The comparison with the corresponding linear bounds given in Table 4 shows that the two types of bounds produce close results. As in the estimation of the other elements of the singular value decomposition, the global perturbation bounds are slightly larger than the corresponding asymptotic estimates but give guaranteed bounds on the changes of the respective elements although for limited size of $\parallel \delta A\parallel$.

10. Two Higher Order Examples

In this section, we present the result of two numerical experiments with higher order matrices to illustrate the properties of the asymptotic and global estimates obtained in the paper.

Example 2.

Consider a $50\times 20$ matrix A, taken as

$$A=U_0\left[\begin{array}{c}{\Sigma }_0\\ 0\end{array}\right]V_0^T,$$

where ${\Sigma }_0=\mathrm{diag}(1,1,\dots ,1)$, the matrices $U_0$ and $V_0$ are constructed as proposed in [1],

$$\begin{array}{ccc}\hfill U_0& =& M_2{S}_U{M}_1,\hfill \\ \hfill M_1& =& I_m-2e{e}^T/m,M_2=I_m-2f{f}^T/m,\hfill \\ \hfill e& =& {[1,1,1,\dots ,1]}^T,f={[1,-1,1,\dots ,{(-1)}^{m-1}]}^T,\hfill \\ \hfill S_U& =& \mathrm{diag}(1,\sigma ,{\sigma }^2,\dots ,{\sigma }^{m-1}),\hfill \\ & & \\ \hfill V_0& =& N_2{S}_V{N}_1,\hfill \\ \hfill N_1& =& I_n-2u{u}^T/n,N_2=I_n-2v{v}^T/n,\hfill \\ \hfill g& =& {[1,1,1,\dots ,1]}^T,h={[1,-1,1,\dots ,{(-1)}^{n-1}]}^T,\hfill \\ \hfill S_V& =& \mathrm{diag}(1,\tau ,{\tau }^2,\dots ,{\tau }^{n-1}),\hfill \end{array}$$

and the orthogonal and symmetric matrices $M_1,M_2,N_1,N_2$ are Householder reflections. The condition numbers of $U_0$ and $V_0$ with respect to the inversion are controlled by the variables σ and τ and are equal to ${\sigma }^{m-1}$ and ${\tau }^{n-1}$, respectively. In the given case, $\sigma =1.2$, $\tau =1.3$ and $\mathrm{cond}\left(U_0\right)=7583.7,\mathrm{cond}\left(V_0\right)=146.2$. The minimum singular value of the matrix A is ${\sigma }_{min}\left(A\right)=0.114319461998077$. The perturbation of A is taken as $\delta A={10}^c\times A_0$, where c is a negative number and $A_0$ is a matrix with random entries generated by the MATLAB^®function rand .

In Figure 5, Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10 we show several results related to the perturbations of the singular value decomposition of A for 30 values of c between $-12$ and $-3$. As particular examples, in Figure 5 we display the perturbations of the entry $U_{20,15}$, which is an element of the matrix $U_1$ and in Figure 6 - the perturbations of the entry $V_{20,15}$, both as functions of ${\parallel \delta A\parallel }_2$. The componentwise linear bound reflect correctly the behavior of the actual perturbations and are valid for wide changes of the perturbation size. Note that this holds for all elements of U and V. The global (nonlinear) bounds practically coincide with the linear bounds but do not exist for perturbations whose size is larger than $2\times {10}^{-3}$.

In Figure 7 and Figure 8 we show the angles between the perturbed and unperturbed left

$${\mathcal{U}}_{15}=\mathrm{span}(u_1,u_2,\dots ,u_{15})$$

and right

$${\mathcal{V}}_{15}=\mathrm{span}(v_1,v_2,\dots ,v_{15})$$

singular subspaces of dimension 15. Again, the linear bounds on the angles $|\Phi max({\mathcal{U}}_{15},{\tilde{\mathcal{U}}}_{15}\left)\right|$ and $|\Theta max({\mathcal{V}}_{15},{\tilde{\mathcal{V}}}_{15}\left)\right|$ are valid for perturbation magnitudes from ${\parallel \delta A\parallel }_2={10}^{-13}$ to ${\parallel \delta A\parallel }_2={10}^{-1}$ and this also holds for singular subspaces of other dimensions.

In Figure 9 and Figure 10 we show the perturbations of the singular values and their nonlinear bounds for perturbations with ${\parallel \delta A\parallel }_2=2.7282784\times {10}^{-10}$ and ${\parallel \delta A\parallel }_2=4.3932430\times {10}^{-4}$. While in the first case the nonlinear bound $\delta {\sigma }_i^{nonl}$ coincides with the actual change $\delta {\sigma }_i$ of the singular values, in the second case the bound becomes significantly greater than the actual change due to the overestimating of the higher order term ${\Delta }_1$.

Example 3.

Consider a $150\times 100$ matrix A, constructed as in the previous example for $\sigma =1.05$ and $\tau =1.0$. The perturbation of A is taken as $\delta A={10}^{-9}\times A_0$, where $A_0$ is a matrix with random entries.

In Figure 11 and Figure 12 we show the entries of the matrices $|\delta U_1|$ and $\left|\delta V\right|$, respectively, along with the corresponding componentwise estimates $\delta U^{lin}$ and $\delta V^{lin}$. The absolute values of the exact changes in all 100 singular values along with the bounds $|\delta {\sigma }_i^{nonl}|$ and ${\parallel \delta A\parallel }_2$ are shown in Figure 13. The nonlinear bounds of $\delta U_1$ and $\delta V$ are found only for 7 iterations and are visually indistinguishable from the corresponding linear bounds.

11. Conclusions

The SVD perturbation analysis presented in this paper makes possible to determine componentwise perturbation bounds of the orthogonal matrices, singular values and singular subspaces of a full rank matrix. The analysis performed has some peculiarities which make it a challenging problem. On one hand, the SVD analysis is simpler than some other problems, like the QR decomposition perturbation problem. This is due to the diagonal form of the decomposed matrix which, among the other, allows to solve easily the equations for the perturbation parameters avoiding the use of the Kronecker product. On the other hand, the presence of two matrices in the decomposition requires the introduction of two different parameter vectors which are mutually dependent due to the relationship between the perturbations of the two orthogonal matrices. This makes necessary to solve a coupled system of equations about the parameter vectors which complicates the analysis.

The analysis presented in the paper can be extended with minor complications to the case of complex matrices.