Some Procrustes-type Problems

Procrustes in the Greek mythology fits his guests to his bed. The mathematical problem after his name has a similar flavour. Suppose we have two set of vectors $X$, $Y$, where columns represent the vectors, and we want to find a linear transform $Q$ that takes $X$ to $Y$, $$ QX = Y $$ meaning we assume the column-to-column correspondence. If there is no constraint on what kind of transform we would like to have, then the answer can be $$ Q = YX^{\dagger} $$ where $X^{\dagger} = X^{T}(XX^{T})^{-1}$ assuming $XX^{T}$ invertible. The question is more instereting when we constrain $Q$ to be a rotation with possible reflection (i.e. orthonormal) $$ Q^{T}Q = I $$ which is of practical interests. Now its clear we don't always have $QX=Y$, so it's natural to ask for the solution of $$ \arg\min_{Q^TQ=I} \|QX-Y\|_{F}^2 $$ where $\|\cdot\|_{F}$ denotes the Frobenius norm. In other words, we want the transformed vectors $QX$ to be closest to $Y$ in a least square sense. Below the fold we will describe two ways to solve this problem with closed form solution.

Two approaches to the orthonormal Procrustes Problem

The first way is through standard constrained optimization approach. We consider the Lagrangian $$ L(Q, \Lambda) = \|QX-Y\|_{F}^2 + \sum_{i,j} \Lambda \odot (Q^TQ - I) $$ where $\odot$ denotes the entries-wise multiplication of matrices, sometimes called Hardamard product; the summation is over all indices. To find a critical point of $L$ and hence the minimum of the convex objective, we would need to differentiate $L(Q,\Lambda)$ with respect to entries of $Q$. To make our life easier, we will need to simplify the expression, by using the following identites for $A$, $B$ of size $m\times n$: $$ \begin{aligned} \|A\|_{F}^2 = \text{tr}(A^{T}A) \\ \text{tr}(A^TB) = \text{tr}(BA^T) = \sum_{i,j} A \odot B \\ \frac{\partial}{\partial X_{ij}}\text{tr}(X^TAX) = ((A^T + A)X)_{ij} \\ \frac{\partial}{\partial X_{ij}}\text{tr}(AX) = A_{ji} \end{aligned} $$ First, we write $$ L(Q,\Lambda) = \text{tr}(QXX^TQ^T) - 2\text{tr}(XY^TQ) + \text{tr}(Y^TY) + \text{tr}(\Lambda^T(Q^TQ-I)). $$ Differentiate in $Q$ gives $$ \frac{\partial L}{\partial Q} = 2QXX^T - 2YX^T + Q(\Lambda^T + \Lambda) $$ Setting this to zero, we have $$ (\Lambda^T + \Lambda)/2 + XX^T = Q^TYX^T $$ Now here is the crux. This means $Q^TXY^T$ must be symmetric $$ Q^TYX^T = XY^TQ $$ Multiplying both sides by $Q$, we get $$ YX^T = QXY^TQ $$ Squaring both sides, $$ YX^TXY^T = QXY^TYX^TQ^T $$ Suppose the eigen decopositions are $$ YX^TXY^T = VDV^T, XY^TYX^T = UDU^T $$ since we know they have the same eigenvalues, we must have $$ V = QU $$ or $Q = VU^T$. We can compute the matrices $U$ and $V$ by computing the singular value decomposition of $XY^T$ $$ XY^T = U\Sigma V^T $$ In the second approach we will effectively prove $Q = UV^T$ is the solution. We write $$ \|QX-Y\|_{F}^2 = \|QX\|_{F}^2 - 2\text{tr}(XY^TQ) + \|Y\|_{F}^2 $$ The last term has nothing to do with $Q$, so can be safely ignored. Now we invoke a lemma.

Lemma 1 If $Q$ is an orthonormal matrix, then $\|QX\|_{F}^2 = \|X\|_{F}^2$ for any matrix $X$.

Hence in fact the first term does not matter as well. We now want to maximize $$ \text{tr}(XY^TQ) $$ We proceed by taking the singluar value decomposition of $XY^T= U\Sigma V^T$, to get $$ \begin{aligned} \text{tr}(XY^TQ) &= \text{tr}(U\Sigma V^TQ) \\ &= \text{tr}(\Sigma V^TQU) \end{aligned} $$ Note that $\Sigma$ is diagonal and recall that they are all positive. So we are instereted in maximizing the diagonal entries of $V^TQU$.

Lemma 2 If $u, v$ are unit vectors, then $\langle u, v\rangle$ is at most $1$.

Viewing $V^TQ$ together, the first diagonal entry in the dot product between the fisrt row of $V^TQ$ and the first column of $U$, and it is at most $1$. It is then easy to see that the diagonal entries are maximized when $$ V^TQU = I $$ We again arrive at the solution $Q = VU^T$.

Relation to As-Rigid-As-Possible (ARAP) deformation

In ARAP, an important subproblem to solve is the following $$ \arg\min_{Q^TQ = I} \|\nabla f - Q\|_F^2 \,\,\,(1) $$ where $f$ is a mapping between two triangulated surfaes and we usually assume $\det \nabla f$ to be positive. Since locally on a triangle, $\nabla f$ maps edges to edges linearly. In the language of Procrustes problem, we are actually looking for $$ \arg\min_{Q^TQ = I} \|YX^{-1}- Q\|_F^2 $$ where we assume $X, Y$ are square and invertible. It is curious to note that this problem is not equivalent to the original Procrustes problem $$ \arg\min_{Q^TQ=I} \|QX-Y\|_{F}^2 \,\,\,(2) $$ which is actually equivalent to $$ \arg\min_{Q^TQ=I} \|YX^T-Q\|_{F}^2 $$ It turns out we can relate Problem (1) to Problem (2) by reweighing the quadratic energy column-wisely with respect to the geometry of the domain vectors, namely the cotangent weight. We refer to a previous post for details on that formulation.

 A more general problem 

Now we focus on Problem (1). Suppose we allow $Q$ to be different from a orthonormal a bit. To be precise, we set bounds to the singular values $\sigma_1, \sigma_2$ of $Q$, namely $$ k \leq \sigma_2 < \sigma_1 \leq K $$ We again ask what would be the best $Q$ that maps $X$ to $Y$ in a least square sense. It's tempting to suggest that the answer is $Q = U \Sigma^{*}V^T$, where $\Sigma^{*}$ is obtained by truncating the sigular values of $\nabla f=U \Sigma V^T$ to be inside the range $[k, K]$, i.e. $\sigma_i^* = \min(\max(\sigma_i, k), K)$. This is indeed true. To see it, we write $$ \begin{aligned} \|\nabla f - Q\|_F^2 & = \|U \Sigma V^T - Q\|_F^2 \\ & = \| \Sigma - U^T Q V\|_F^2 \end{aligned} $$ Let $Q = S\Sigma' P^T$ be the sigular value decomposition of $Q$. The above expression is $$ \| \Sigma - U^T S\Sigma' P^T V\|_F^2 =: \| \Sigma - U'^T\Sigma'V'\|_F^2 $$ which is equal to $$ \text{tr}(\Sigma^T \Sigma) - 2\text{tr}(\Sigma^TV'^T\Sigma'U') + \text{tr}(\Sigma'^T\Sigma') \,\,\,\,(3) $$ Note that only the middle term has to do with the rotation matrices. The problem now turns into maximizing $$ \text{tr}(\Sigma^TU'^T\Sigma'V') $$ over orthonormal matrices $U', V'$. Suppose we have ordered the diagonal entries of $\Sigma$ from large to small. Since they are all positive, we would like to maximize the diagonal entries of $U'^T\Sigma'V'$ and order them from large to small.

Proposition 3 Let $U',V'$ be orthonormal matrices, $\Sigma'$ be a matrix of zeros but non-negative on the diagonal. Suppose without loss of generality that the diagonal of $\Sigma'$ is ordered from large to small $\sigma_1' \geq \sigma_2' \geq \cdots$. Then we have $$ \text{diag}(U'^T\Sigma' V') \leq \text{diag}(\Sigma') $$ where we use $\leq$ as entry-wise comparison.

To prove this, observe that the $k$-th diagonal entry of $U'^T\Sigma' V'$ is $$ u_{k}'^T \Sigma' v_{k}' $$ We state some simple lemmas following the above notation, and the proposition will follow immediately.

Lemma 4 If $u, v$ are unit vectors, then $u^T \Sigma' v$ is at most $\sigma_1'$.

Lemma 5 If $u, v$ are unit vectors such that $u^T \Sigma' v = \sigma_1'$, then for any unit vectors $u', v'$ orthogonal to $u, v$ repectively, $u'^T \Sigma' v'$ is at most $\sigma_2'$.

These lemmas are in fact equivalent to the min-max characterization for singular values. Continue with the expression (3), we see when $U', V'$ are optimal, it is equal to $$ \text{tr}(\Sigma^T \Sigma) - 2\text{tr}(\Sigma^T\Sigma') + \text{tr}(\Sigma'^T\Sigma') $$ which is minimized when $\Sigma'$ is thresholded as described.

MATH 3310 Note 11: Gradient and conjugate gradient methods

1. The steepest descent

Solving a linear system

$ \displaystyle B{\bf x}={\bf c} $

can be recast as a optimization problem

$ \displaystyle \min_{{\bf x}\in\mathbb{R}^{n}}f({\bf x})=\|B{\bf x}-{\bf c}\|^{2}. \ \ \ \ \ (1)$

One easily sees that a necessary condition for the optimum is

$ \displaystyle B^{T}B{\bf x}=B^{T}{\bf c}. $

This is a normal equation, and is equivalent to the original equation if $ {B}$ is square and invertible (although the condition number of $ {B^{T}B}$ becomes worse, and it requires twice the time to compute the matrix vector product). It is also possible for $ {B}$ to only have linear independent columns (so there are more rows than columns), giving the least square solution we have seen.
So we are interested in linear systems of the form

$ \displaystyle A{\bf x}={\bf b} $

where $ {A\in\mathbb{R}^{n\times n}}$ is symmetric positive definite (SPD).
One simple way to optimize (1) is a greedy algorithm. Suppose we are at $ {{\bf x}^{(k)}}$. We want to go for the best direction $ {{\bf p}^{(k)}}$ using only the local information of the objective function $ {f}$. This will require some smoothness assumption on $ {f}$, so that for small fixed step sizes $ {\alpha}$, the choice $ {{\bf p}^{(k)}}$ always win over all other directions (unless we are already at the optimum). Fortunately, this direction is quite easy to find.

Exercise 1 (Steepest descent direction) Suppose ${f:\mathbb{R}^{n}\rightarrow\mathbb{R}}$ is differentiable and let $ {\varphi(t)=f({\bf x}+t{\bf y})}$. Show that

$ \displaystyle -\nabla f(x)=\arg\min_{\underset{\|y\|=1}{{\bf y}\in\mathbb{R}^{n}}}\varphi'(0)=\nabla f(x)^{T}{\bf {\bf y}}. $

We are thus led to the update formula

$ \displaystyle {\bf x}^{(k+1)}={\bf x}^{(k)}+\alpha{\bf p}^{(k)},\quad\alpha>0\text{ small} $

This is known as the gradient descent algorithm. The step size requires a choice, and it is natural to find $ {\alpha^{*}}$ that gives the best result. This is a one dimensional problem, and when the objective function is in simple form such as (1), we have a closed form formula. Otherwise, it is find by a search algorithm.

Exercise 2 (Optimal step size) Show that

$ \displaystyle \alpha_{k}^{*}=-\frac{(A{\bf x}^{(k)}-{\bf b})^{T}{\bf p}_{k}}{{\bf p}^{(k)T}A{\bf p}^{(k)}} $

where

$ \displaystyle \alpha_{k}^{*}=\arg\min_{\alpha>0}f({\bf x}^{(k)}+\alpha{\bf p}^{(k)}). $

Gradient descent and its variants can be applied to quite general objective functions (even differentiablity can be somewhat relaxed), and is thus widely used in practice. But when the objective function is quadratic as in (1), there is a lot of room for improvement. We shall next turn to the most widely used method for solving SPD linear system.

2. The conjugate gradient algorithm

The conjugate gradient method for a positive definite quadratic optimization problem (equivalently solving normal equations) is very special -- it hits the exact solution in at most $ {n}$ steps! (Only in theory, though, in practice there is the devil of round-off errors). This dramatically distinguishes itself from other iterative schemes for linear system, which almost never give the exact solution but very good approximation of it. Conjugate gradient is a semi-direct method, in that it is almost a direct method like Gaussian elimination which yields the exact solution in finite number of steps, but more importantly one can also stop earlier and get an approximated solution as well.

Let's go to an example to motivate the ideas, where things are quite direct. We restrict our attention on the SPD tridiagonal system first

$ \displaystyle T{\bf y}={\bf f} $

where

$ \displaystyle T=\begin{pmatrix}d_{1} & u_{1}\\ u_{1} & d_{2} & u_{2}\\ & \ddots & \ddots & \ddots\\ & & u_{n-2} & d_{n-1} & u_{n-1}\\ & & & u_{n-1} & d_{n} \end{pmatrix},\quad{\bf f}=\begin{pmatrix}\|{\bf b}\|\\ 0\\ \vdots\\ \\ 0 \end{pmatrix}. $

The choice of $ {{\bf f}}$ is not important.
We could solve for the $ {{\bf y}}$ directly, but that misses the point of conjugate gradient. The idea now is that in each step, we only use a principle submatrix of $ {T}$ to solve for the iterate $ {{\bf y}^{(k)}}$. In the first step, we have

$ \displaystyle d_{1}y_{1}^{(1)}=\|{\bf b}\|,\quad{\bf y}^{(1)}=\begin{pmatrix}y_{1}^{(1)}\\ 0\\ \vdots\\ 0 \end{pmatrix} $

in the second step,

$ \displaystyle \begin{pmatrix}d_{1} & u_{1}\\ u_{1} & d_{2} \end{pmatrix}\begin{pmatrix}y_{1}^{(2)}\\ y_{2}^{(2)} \end{pmatrix}=\begin{pmatrix}\|{\bf b}\|\\ 0 \end{pmatrix},\quad{\bf y}^{(2)}=\begin{pmatrix}y_{1}^{(2)}\\ y_{2}^{(2)}\\ 0\\ \vdots\\ 0 \end{pmatrix} $

and in the third step,

$ \displaystyle \begin{pmatrix}d_{1} & u_{1} & 0\\ u_{1} & d_{2} & u_{2}\\ 0 & u_{2} & d_{3} \end{pmatrix}\begin{pmatrix}y_{1}^{(3)}\\ y_{2}^{(3)}\\ y_{3}^{(3)} \end{pmatrix}=\begin{pmatrix}\|{\bf b}\|\\ 0\\ 0 \end{pmatrix},\quad{\bf y}^{(3)}=\begin{pmatrix}y_{1}^{(3)}\\ y_{2}^{(3)}\\ y_{3}^{(3)}\\ 0\\ \vdots\\ 0 \end{pmatrix} $

and so on. We would like to relate $ {{\bf y}^{(k+1)}}$ to the previous $ {{\bf y}^{(k)}}$, and it is the basic step of the conjugate gradient method.

Remark 1 Several remarks are in order.

1. By definition, the exact solution $ {{\bf y}={\bf y}^{(n)}}$ is reached in at most $ {n}$ steps.
2. Each $ {{\bf y}^{(k)}}$ is the "least square'' choice among all vectors whose last $ {n-k}$ components are zero. In other words,
   $ \displaystyle {\bf y}^{(k)}=\arg\min_{\underset{y_{i}=0\,\forall i>k}{{\bf y}\in\mathbb{R}^{n}}}\frac{1}{2}{\bf y}^{T}T{\bf y}-{\bf f}^{T}{\bf y}. $
3. The direction of change $ {{\bf y}^{(k+1)}-{\bf y}^{(k)}}$ is "$ {T}$-conjugate'' to all previous and later directions:
   $ \displaystyle ({\bf y}^{(k+1)}-{\bf y}^{(k)})^{T}T({\bf y}^{(j+1)}-{\bf y}^{(j)})=0\quad j\neq k. $
4. The residuals $ {{\bf r}^{(k)}=T{\bf y}^{(k)}-{\bf f}}$ from each step are orthogonal to each other
   $ \displaystyle {\bf r}^{(k)T}{\bf r}^{(j)}=0\quad j\neq k. $

The third remark is not trivial and requires a proof. By symmetry of $ {T}$, it suffices to prove the statement for $ {j>k}$. We compare the entries of the vectors $ {{\bf y}^{(k+1)}-{\bf y}^{(k)}}$ and $ {T({\bf y}^{(j+1)}-{\bf y}^{(j)})}$. The vector $ {{\bf y}^{(k+1)}-{\bf y}^{(k)}}$ is nonzero only at the first $ {k+1}$ entries. As for the other vector, we observe that

$ \displaystyle T{\bf y}^{(j)}=\begin{pmatrix}\|{\bf b}\|\\ 0\\ \vdots\\ 0\\ f_{j+1}^{(j)}\\ 0\\ \vdots \end{pmatrix}, $

and so

$ \displaystyle T({\bf y}^{(j+1)}-{\bf y}^{(j)})=\begin{pmatrix}\|{\bf b}\|\\ 0\\ \vdots\\ 0\\ 0\\ f_{j+2}^{(j+1)}\\ \vdots \end{pmatrix}-\begin{pmatrix}\|{\bf b}\|\\ 0\\ \vdots\\ 0\\ f_{j+1}^{(j)}\\ 0\\ \vdots \end{pmatrix}=\begin{pmatrix}0\\ 0\\ \vdots\\ 0\\* \\* \\ 0\\ \vdots \end{pmatrix} $

is nonzero only at the $ {(j+1)}$ and $ {(j+2)}$-th entries. We see that the nonzero entries of the two vectors don't overlap, hence the remark is proved. By the way, we also see that the nonzero entries of $ {T{\bf y}^{(k)}-{\bf f}}$ and $ {T{\bf y}^{(j)}-{\bf f}}$ does not overlap, so they are orthogonal as well, proving the last remark.
It turns out the above method can be extended to all SPD matrices as well, by tridiagonalizing the SPD matrices. There is a general technique that applies to all such matrices. We won't need to go into the full theory here. But the point is that the technique involves orthogonalizing the set of vectors produced by the matrix powers

$ \displaystyle \begin{pmatrix}{\bf b} & A{\bf b} & A^{2}{\bf b} & \cdots & A^{n-1}{\bf b}\end{pmatrix}. $

This leads us to the Krylov subspace method.

Definition 1 A sequence of nested subspaces ${\mathcal{K}_{0}\subset\mathcal{K}_{1}\subset\cdots}$ is called the Krylov subspaces associated to the equation $ {A{\bf x}={\bf b}}$ if

$ \displaystyle \mathcal{K}_{0}=\{0\},\quad\mathcal{K}_{k}=\text{span}\{b,Ab,\dots A^{k-1}b\}. $

Note that if $ {\mathcal{K}_{k+1}=\mathcal{K}_{k}}$ for some ${k}$, then $ {\mathcal{K}_{j}=\mathcal{K}_{k}}$ for all $ {j>k}$. The sequence of vectors

$ \displaystyle {\bf x}^{(k)}=\arg\min_{{\bf x}\in\mathcal{K}_{k}}\frac{1}{2}{\bf x}^{T}A{\bf x}-{\bf b}^{T}{\bf x} $

is called a Krylov sequence.

Conjugate gradient is a recursive method for computing the Krylov sequence. One should recognize that the vectors $ {{\bf y}^{(k)}}$ in the earlier example is a Krylov sequence for the tridiagonal system $ {T{\bf y}={\bf f}}$, which follows from the second remark above. The corresponding properties still hold true for a general Krylov sequence. We summerise it here.

Proposition 2 Using the same notation from above, the following statements hold

1. The residual $ {{\bf r}^{(k)}}$ is such that $ {{\bf r}^{(k)}\in\mathcal{K}_{k}^{\perp}}$ and $ {{\bf r}^{(k)}\in\mathcal{K}_{k+1}}$. Hence the residuals form an orthogonal basis of the corresponding Krylov subspaces.
2. The steps $ {{\bf d}^{(k)}={\bf x}^{(k)}-{\bf x}^{(k-1)}}$ are "$ {A}$-conjugate'' to each other
   $ \displaystyle {\bf d}^{(k)T}A{\bf d}^{(j)}=0,\quad j\neq k. $
   Hence the steps form an "$ {A}$-conjugate'' basis for the corresponding Krylov subspaces.
3. We have
   $ \displaystyle {\bf d}^{(k)T}A{\bf d}^{(k)}={\bf d}^{(k)T}{\bf r}^{(k-1)},\quad\forall k. $

From the first and second statements above, we can see that $ {{\bf d}^{(k)}}$ can be expressed in terms of $ {{\bf d}^{(1)},\dots,{\bf d}^{(k-1)}}$ and $ {{\bf r}^{(k-1)}}$, since

$ \displaystyle {\bf d}^{(k)}\in\mathcal{K}_{k}=\text{span}\{{\bf d}^{(1)},\dots,{\bf d}^{(k-1)},{\bf r}^{(k-1)}\}. $

In fact, if we write

$ \displaystyle {\bf p}^{(k)}=\frac{{\bf r}^{(k-1)T}{\bf r}^{(k-1)}}{{\bf d}^{(k)T}{\bf r}^{(k-1)}}{\bf d}^{(k)}, $

then

$ \displaystyle {\bf p}^{(k)}=-{\bf r}^{(k-1)}-\frac{{\bf p}^{(k-1)T}A{\bf r}^{(k-1)}}{{\bf p}^{(k-1)T}A{\bf p}^{(k-1)}}{\bf p}^{(k-1)}, $

which is the most important step of updating the direction in the conjugate gradient algorithm.
As a final remark, we mention that conjugate gradient converges fast if the matrix is well-conditioned, that means the eigenvalues of $ {A}$ cluster in the small region. In particular, repeated eigenvalues leads to exact solution in less than $ {n}$ steps. Hence it is often the case in practice that the matrix $ {A}$ has to be preconditioned. For more about this, you can refer to the notes by Prof. Vandenbergh here.

MATH3310 Note 10: Numerical methods for linear eigenvalue problems

1. The power method: finding the eigenvalues and eigenvectors

The power method for eigenvalue problem takes the advantage of the discrepancy between the eigenvalues. Strictly speaking, it is the difference in modulus that matters. Suppose an vector $ {{\bf u}}$ can be written as linear combination of the eigenvectors $ {\{{\bf v}_{i}\}}$ of a matrix $ {A}$

$ \displaystyle {\bf u}=u_{1}{\bf v}_{1}+u_{2}{\bf v}_{2}+\cdots+u_{n}{\bf v}_{n},\quad u_{i}\neq0\ \forall i. \ \ \ \ \ (1)$

Then, by linearity

$ \displaystyle A^{k}{\bf u}=\lambda_{1}^{k}u_{1}{\bf v}_{1}+\lambda_{2}^{k}u_{2}{\bf v}_{2}+\cdots\lambda_{n}^{k}u_{n}{\bf v}_{n}. $

If, say, $ {\lambda_{2}}$ is the only largest one in modulus (let's call it dominant), then the component $ {\lambda_{2}^{k}u_{2}{\bf v}_{2}}$ should stand out more significantly than any other component as $ {k}$ becomes larger. Note that we have no knowledge about what the eigenvectors or eigenvalues are, except that we assume

• the vector $ {{\bf u}}$ can be written in the form of (1).
• there exists one and only one largest eigenvalue.

Of course, in general if we simply iterate $ {A^{k}{\bf u}}$, we don't have convergence in the usual sense (it can blow-up or vanish). We must introduce a rescaling, or normalization, namely to divide the vector $ {A^{k}u}$ by a suitable scalar. There are many choices for the scalar, one example is the largest component in absolute value $ {\|A^{k}u\|_{\infty}}$. In general, it can be any norm of the vector. Then one goes on again multiplying the matrix $ {A}$ and rescaling the vector. This simple algorithm is called the power iteration. In the limit we recover that dominant component.
There are a number of variants of the power iteration, designed to find eigenvalues other than the dominant one, and to accelerate convergence.

• The inverse power iteration. Here we iterate $ {{\bf u}^{(k)}=A^{-1}{\bf u}^{(k-1)}}$ to find the least eigenvalue of $ {A}$ in modulus, of course assuming $ {A}$ is invertible. Instead of find the inverse matrix $ {A^{-1}}$, we solve the equation
  $ \displaystyle A{\bf u}^{k}={\bf u}^{(k-1)}. $
• Power iteration with shift. Here the matrix $ {A}$ is replaced by $ {A-\mu I}$. By doing this, we see a shift in the eigenvalues of $ {A}$. If we have a good estimation $ {\mu}$ of the eigenvalue $ {\lambda_{i}}$ of interest, then we can simply subtract $ {\mu}$ on the diagonal to get the matrix $ {A-\mu I}$, so that the shifted eigenvalue of interest $ {\lambda_{i}-\mu}$ is the least eigenvalue in modulus. And thus we can apply the inverse power iteration to find $ {\lambda_{i}-\mu}$.
  The shift parameter $ {\mu}$ does not have to be the same for each iteration step. One way to update it is to use the Rayleign Quotient
  $ \displaystyle \frac{{\bf u}^{(k)T}A{\bf u}^{(k)}}{{\bf u}^{(k)T}{\bf u}^{(k)}}. $
  It is a good choice of estimating $ {\lambda_{i}}$ when $ {{\bf u}^{(k)}}$ is a good approximation of the eigenvector corresponding to $ {\lambda_{i}}$. Note that Rayleigh Quotient method needs stronger condition than the power iteration: the angle between the iteration vector and the true eigenvector has to be small, not merely containing a non-trivial component. Otherwise it may converge to a different eigenvalue/eigenvector than expected.
• Simultaneous iteration. Here we need the matrix to be normal. Instead of interating a single vector $ {{\bf u}^{(k)}}$, here we iteration several vectors $ {U^{(k)}=\begin{pmatrix}{\bf u}_{1}^{(k)} & {\bf u}_{2}^{(k)} & \cdots & {\bf u}_{p}^{(k)}\end{pmatrix}}$. How can we get different eigenvalues from such an iteration? The answer is to do something like a Gramm-Schmidt orthogonalization procedure on the iterated matrix $ {U^{(k)}=Q^{(k)}R^{(k)}}$, and keep only the orthogonal matrix part $ {Q^{(k)}}$, and start iterating again. This scheme in a way can lead to the QR iteration, which we shall see next. Because of the orthogonalization, the iteration is numerically very stable (as the condition number of an orthogonal matrix is $ {1}$, so round off errors cannot amplify). 

2. The QR algorithm: finding the eigenvalues and the invariant subspaces

2.1. The QR decomposition. The basic component of the QR iteration is the QR factorisation, which in theory can be achieved by the Gramm-Schmidt orthogonalization. The standard QR decomposition of a matrix $ {A\in\mathbb{R}^{m\times n}}$ has the form

$ \displaystyle A=QR,\quad Q\in\mathbb{R}^{m\times n},R\in\mathbb{R}^{n\times n}, $

where $ {Q}$ has orthonormal columns, and $ {R}$ is a upper triangular matrix. But computing the matrix $ {Q}$ following the Gramm-Schmidt recipe can be numerically unstable (because subtracting away the all components at once can leave us non-orthogonal results due to numerical error). There is a modified Gramm-Schmidt procedure to the rescue, using a one-dimensional procedure each time. However, here we want to talk about something more preferable in practice.

The idea is to transform the matrix $ {A}$ to some upper triangular matrix $ {\hat{R}}$ of the same size, using elements in the orthogonal group $ {O(m)}$. In other words, we want to find a finite sequence of orthogonal matrices $ {H_{1},H_{2}\cdots H_{p}}$ such that

$ \displaystyle H_{p}\cdots H_{2}H_{1}A=\hat{R}. $

Then we obtain a variant of QR factorisation of $ {A}$ by

$ \displaystyle A=H_{1}^{T}H_{2}^{T}\cdots H_{p}^{T}\hat{R}=\hat{Q}\hat{R}, $

where $ {\hat{Q}\in\mathbb{R}^{m\times m}}$ and $ {\hat{R}\in\mathbb{R}^{m\times n}}$. Compared to the original standard QR decomposition, the $ {Q}$ is completed to an orthogonal matrix $ {\hat{Q}}$, and zeros are attached to the additional rows (if there are any) of $ {R}$. Since the matrices $ {H_{i}}$ are orthogonal, the procedure is much more stable.

We haven't spoken about the choice of these matrices. A particularly useful and simple choice is the Householder reflection. As the name suggests, it is a reflection in the space $ {\mathbb{R}^{m}}$, with respect to certain $ {m-1}$ dimensional hyperplane containing the orgin. This hyperplane is defined by a unit normal vector $ {{\bf u}}$ to it. The Householder reflection associated to the unit vector $ {{\bf u}}$ is defined to be

$ \displaystyle H=I-2{\bf u}{\bf u}^{T}. $

As a first check, note that

$ \displaystyle H{\bf u}=-{\bf u}. $

And if $ {{\bf x}}$ is such that $ {{\bf u}^{T}{\bf x}=0}$, then $ {H{\bf x}={\bf x}}$. Hence $ {H}$ is a reflection. Also, $ {H^{T}=H}$.

Exercise 1 Show directly that the above matrix $ {H}$ is orthogonal.

Why Householder matrices are useful in our context? Suppose

$ \displaystyle A=\begin{pmatrix}\\ {\bf a}_{1} & {\bf a}_{2} & \cdots & {\bf a}_{n}\\ \\ \end{pmatrix}. $

We can tranform the first column of $ {A}$, using a Householder reflection, to $ {\|{\bf a}_{1}\|{\bf e}_{1}}$, namely

$ \displaystyle H_{1}{\bf a}_{1}=\begin{pmatrix}\|{\bf a}_{1}\|\\ 0\\ \vdots\\ 0 \end{pmatrix}={\bf r}_{1}. $

The matrix $ {H_{1}}$ can be found by noting that the vector $ {{\bf a_{1}}-{\bf r}_{1}}$ should be perpendicular to the hyperplane, and so the defining unit vector of $ {H_{1}}$ should be

$ \displaystyle \frac{{\bf a_{1}}-{\bf r}_{1}}{\|{\bf a_{1}}-{\bf r}_{1}\|}. $

Exercise 2 Let $ {{\bf u=}\frac{{\bf a_{1}}-{\bf r}_{1}}{\|{\bf a_{1}}-{\bf r}_{1}\|}}$, ${H=I-2{\bf u}{\bf u}^{T}}$. Show that

$ \displaystyle H{\bf a}_{1}={\bf r}_{1}. $

Now, to proceed further, we want to find a matrix $ {H_{2}}$ such that the second column of $ {H_{2}H_{1}A}$ is from a upper triangular matrix. The choice for the normal unit vector is now

$ \displaystyle \frac{{\bf a}-{\bf r}}{\|{\bf a}-{\bf r}\|} $

where

$ \displaystyle {\bf a}^{(2)}=\begin{pmatrix}0\\ H_{1}{\bf a}_{2}(2:end) \end{pmatrix},\quad{\bf r}^{(2)}=\begin{pmatrix}0\\ \|{\bf a}^{(2)}\|\\ 0\\ \vdots\\ 0 \end{pmatrix}. $

You should check that the matrix $ {H_{2}}$ actually does not operate on the first entry in $ {H_{1}{\bf a}_{2}}$. In fact the first dimension remain fixed under the action of $ {H_{2}}$. And so the first column of $ {H_{1}{\bf a}_{1}}$ is not changed. We can keep going in this manner to define $ {H_{3},\dots H_{m}}$, and obtain the variant of QR decomposition we mentioned earlier

$ \displaystyle A=H_{1}H_{2}\cdots H_{m}\hat{R}=\hat{Q}\hat{R}. $

2.2. Solving least square problem using the QR decomposition. The QR decomposition of a "slim matrix'' can be useful for solving the least square problem

$ \displaystyle \min_{{\bf x}\in\mathbb{R}^{n}}\|A{\bf x}-{\bf b}\|^{2} $

where $ {A\in\mathbb{R}^{m\times n}}$, $ {m>n}$ and $ {\|\cdot\|}$ is the usual $ {\ell^{2}}$ norm. This formulation is one way to solve the over-determined system $ {A{\bf x}={\bf b}}$.

Using the standard QR decomposition of $ {A=QR}$, we find

$ \displaystyle QR{\bf x}={\bf b}. $

Now multiplying two sides by $ {Q^{*}}$, and noting that $ {Q^{*}Q=I_{n\times n}}$, we get

$ \displaystyle R{\bf x}=Q^{*}{\bf b}. $

Hence if $ {A}$ has full rank, so $ {R}$ is invertible and we can obtain a unique solution $ {{\bf x}^{*}}$of the above transformed system.

We claim that

$ \displaystyle {\bf x}^{*}=\arg\min_{{\bf x}\in\mathbb{R}^{n}}\|A{\bf x}-{\bf b}\|^{2}. $

It is a quadratic distance problem and we have the Pythagoras theorem at our hand. It amounts to show that the error $ {A{\bf x}^{*}-{\bf b}}$ is orthogonal to the column space $ {\text{col}(A)}$ of $ {A}$. Now $ {A=QR}$, where colums of $ {Q}$ form an orthonormal basis of $ {\text{col}(A)}$. So by multiplying $ {Q^{*}}$ with $ {A{\bf x}^{*}-{\bf b}}$ on the right, we orthogonally project the error to $ {\text{col}(A)}$, and we get

$ \displaystyle Q^{*}(A{\bf x}^{*}-{\bf b})=R{\bf x}^{*}-Q^{*}{\bf b}=0 $

by definition of $ {{\bf x}^{*}}$. Hence the error is orthogonal to $ {\text{col}(A)}$ and the claim is thus proved.

2.3. The QR algorithm.
What does $ {QR}$ decomposition have anything to do with the eignvalues? It is indeed a surprise that a fast and stable method of computing the eigenvalues is based on this decomposition. We describe the algorithm first.

Given a square matrix $ {A^{(0)}=A\in\mathbb{R}^{n\times n}}$, the QR iteration does the following: for $ {k=1,2,\dots}$

1. Do the QR factorisation $ {A^{(k-1)}=Q^{(k)}R^{(k)}}$.
2. Update $ {A^{(k)}=R^{(k)}Q^{(k)}}$.

Under the certain assumptions, the matrix $ {A^{(k)}}$ should eventually approach a triangular matrix, where on the diagonal there lie the eigenvalues of $ {A}$! The QR iteration actually computes the Schur decomposition of $ {A}$:

$ \displaystyle A=Q^{(\infty)*}A^{(\infty)}Q^{(\infty)}, $

where $ {Q^{(\infty)},A^{(\infty)}}$ are the limits of the QR algorithm, being orthogonal and upper triangular, respectively. If $ {A}$ is real symmetric, then $ {A^{(\infty)}}$ is in fact diagonal and $ {Q^{(\infty)}}$ contains the corresponding eigenvectors.

We note that each $ {A^{(k)}}$ is actually similar to $ {A}$, since

$ \displaystyle Q^{(k)*}A^{(k-1)}Q^{(k)}=R^{(k)}Q^{(k)}=A^{(k)} $

so inductively

$ \displaystyle A^{(k)}=Q^{(k)*}\cdots Q^{(1)*}AQ^{(1)}\cdots Q^{(k)}. $

So each $ {A^{(k)}}$ has the same eigenvalues with $ {A}$. Furthermore, we have

$ \displaystyle A^{2}=Q^{(1)}R^{(1)}Q^{(1)}R^{(1)}=Q^{(1)}Q^{(2)}R^{(2)}R^{(1)}, $

so inductively

$ \displaystyle A^{k}=Q^{(1)}\cdots Q^{(k)}R^{(k)}R^{(k-1)}\cdots R^{(1)}. $

We now have a QR decomposition of $ {A^{k}}$! Now $ {Q^{(1)}\cdots Q^{(k)}}$ can be seen to be equivalent to applying the simutaneous iteration to the standard basis $ {I=\begin{pmatrix}{\bf e}_{1} & {\bf e}_{2} & \cdots & {\bf e}_{n}\end{pmatrix}}$. More explanation about its relation to the power iteration and convergence in a special case can be found here.

We make some final remarks about the QR algorithm as stated. It is too inefficient to be effective. There are two refinements that make it more competitive:

1. Reduction of the matrix $ {A}$ to Hessenberg form.
2. Use of shift to accelerate convergence.

We won't pursue these topics here.

MATH3310 Note 9: Iterative methods for linear systems, continued

1. Successive over-relaxation

When the classical methods of Jacobi and Gauss-Seidel converge, they converge monotonically. It is thus natural to go beyond a little in the descending direction, a method known as over-relaxation. One can introduce this technique to a iterative scheme of the following general form

$ \displaystyle Mx^{k+1}=(M-A)x^{k}+b $

as follows. First, define an auxilary sequence $ {\hat{x}^{k}}$ by

$ \displaystyle M\hat{x}^{k+1}=(M-A)x^{k}+b. $

And then update $ {x^{k+1}}$ according to

$ \displaystyle x^{k+1}=x^{k}+\omega(\hat{x}^{k+1}-x^{k}),\quad\omega>0 $

When $ {\omega=1}$, this iteration is the same with the original one. And the idea of the term $ {\omega(\hat{x}^{k+1}-x^{k})}$ is to adjust the correction given by the auxilary iteration. One usually refers to the case $ {\omega>1}$ as over-relaxation, since one enlarges the correction compared to the original iteration.

The successive over-relaxation (SOR) method is a modification of the Gauss-Seidel iteration along this line. Recall the notation

$ \displaystyle A=L+D+U $

where $ {L,D,U}$ are the strict lower triangular, diagonal and strict upper triangular part of $ {A}$, respectively. Then, following the above, we have

$ \displaystyle (D+L)\hat{x}^{k+1}=-Ux^{k}+b $

and

$ \displaystyle x^{k+1}=x^{k}+\omega(\hat{x}^{k+1}-x^{k}). $

Exercise 1 Show that the above scheme can be simplified to

$ \displaystyle (D+\omega L)x^{k+1}=(-\omega U+(1-\omega)D)x^{k}+\omega b. $

This means we have a splitting of $ {A=M+(A-M)}$ where

$ \displaystyle M=L+\frac{1}{\omega}D. $

The benefit of such a scheme is that in the right situation, the convergence rate can be significantly faster. Let's look at an example.

$ \displaystyle A=\begin{bmatrix}2 & -1\\ -1 & 2 \end{bmatrix},\quad M=\begin{bmatrix}2/\omega & 0\\ -1 & 2/\omega \end{bmatrix} $

The iteration matrix is

$ \displaystyle B=\begin{bmatrix}1-\omega & 1/2\omega\\ \omega(1-\omega)/2 & 1-\omega+\omega^{2}/4 \end{bmatrix}. $

We can optimize the eigenvalues in $ {\omega}$, it turns out that the optimal choice is when the two eigenvalues equal to $ {1-\omega^{*}}$, where

$ \displaystyle \omega^{*}=4(2-\sqrt{3})\approx1.07. $

In other words the spectral radius of the matrix $ {B}$ is around $ {\rho(B)=0.07}$. Compare this to the case of Jacobi iteration ($ {\rho(B)=\frac{1}{2}}$) and Gauss-Seidel $ {(\rho(B)=\frac{1}{4}}$). A rule of thumb is that the convergence is fast when all the eigenvalues of the iteration matrix stays about the same level.

Block relaxation scheme. It is easy to extend the previously described iterations to handle linear system in blocks. For example, the block Jacobi iteration for the block matrix $ {A=(A_{ij})_{1\leq i,j\leq p}}$ will be as follows. Define

$ \displaystyle D=\begin{bmatrix}A_{11}\\ & A_{22}\\ & & \ddots\\ & & & A_{pp} \end{bmatrix},\quad L=\begin{bmatrix}0\\ A_{21}\\ \vdots & \ddots\\ A_{p1} & \cdots & A_{p,p-1} & 0 \end{bmatrix}, $

and

$ \displaystyle L=\begin{bmatrix}0 & A_{12} & \cdots & A_{1,p}\\ & & \ddots & \vdots\\ & & & A_{p-1,p}\\ & & & 0 \end{bmatrix}. $

The block Jacobi iteration is of the same form

$ \displaystyle Dx^{k+1}=-(L+U)x^{k}+b. $

Such a scheme can be preferable for a parallel computer. There are more general techniques , dealing with e.g. overlapping blocks, which belong to the family of domain decomposition methods.

Iterative method as preconditioning. A preconditioner for a square linear system $ {Ax=b}$ is a matrix $ {P}$, of the same size with $ {A}$, which is used to produce a equivalent system $ {PAx=Pb}$, with (hopefully) $ {PA}$ having smaller condition number (usually defined as the ratio between the modulus of the largest and smallest eigenvalue). The preconditioned system can be preferable in terms of stability and convergence rate.

We can view the Jacobi, Gauss-Seidel and SOR as fixed point iteration on a precondtioned system

$ \displaystyle M^{-1}A=M^{-1}b, \ \ \ \ \ (1)$

where $ {M}$ is same notation in $ {A=M+(A-M)}$ as before. To see this, recall that the general scheme is of the form

$ \displaystyle x^{k+1}=(I-M^{-1}A)x^{k}+M^{-1}b $

which is the fixed point iteration of

$ \displaystyle (I-I+M^{-1}A)x=M^{-1}b $

which is equivalent to (1). These are particular simple cases of the preconditioned iterative methods.

2. Consistenly ordered matrices

In class, we called a matrix $ {A}$ to be consistently ordered if eigenvalues of the matrices

$ \displaystyle \alpha L+\frac{1}{\alpha}U,\quad\alpha\neq0 $

are independent of $ {\alpha}$. We are interested in this property because the convergence analysis of SOR scheme for such matrices are particularly simple. Here, we would like to find a special class of matrices that satisfy this property.
We have in mind the tridiagonal matrices, which are of the form

$ \displaystyle A=\begin{bmatrix}d_{1} & t_{12}\\ t_{21} & d_{2} & t_{23}\\ & t_{32} & d_{3} & t_{34}\\ & & \ddots & \ddots & \ddots\\ & & & t_{N,N-1} & d_{N} \end{bmatrix}. $

Examples include the central difference matrix for $ {\frac{d^{2}}{dx^{2}}}$ with Dirichlet or Neumann boundary condition. We note that, the matrices

$ \displaystyle L+U=\begin{bmatrix}0 & t_{12}\\ t_{21} & 0 & t_{23}\\ & t_{32} & 0 & t_{34}\\ & & \ddots & \ddots & \ddots\\ & & & t_{N,N-1} & 0 \end{bmatrix} $

and

$ \displaystyle \alpha L+\frac{1}{\alpha}U=\begin{bmatrix}0 & \frac{1}{\alpha}t_{12}\\ \alpha t_{21} & 0 & \frac{1}{\alpha}t_{23}\\ & \alpha t_{32} & 0 & \frac{1}{\alpha}t_{34}\\ & & \ddots & \ddots & \ddots\\ & & & \alpha t_{N,N-1} & 0 \end{bmatrix} $

are related by a similarity transform:

$ \displaystyle L+U=X^{-1}\left(\alpha L+\frac{1}{\alpha}U\right)X $

where

$ \displaystyle X=\begin{bmatrix}1\\ & \alpha\\ & & \alpha^{2}\\ & & & \ddots\\ & & & & \alpha^{N} \end{bmatrix}. $

It then follows that the eigenvalues of $ {\alpha L+\frac{1}{\alpha}U}$ is the same with those of $ {L+U}$.

It turns out that this argument can be extended to block tridiagonal matrices whose diagonal blocks are diagonal matrices. Matrices with this property is called T-matricex, and can be further extended to similarity transforms of such matrices.

Now we recall the result we proved in class.

Theorem 1 Let $ {A}$ be given. Suppose

1. $ {0<\omega<2}$.
2. The Jacobi iteration converges, i.e. $ {\rho(D^{-1}(L+U))<1}$ and the matrix $D^{-1}(L+U)$ has real eigenvalues.
3. A is consistenly ordered.

Then the SOR method converges. The optimal $ {\omega}$ is given by

$ \displaystyle \omega^{*}=\frac{2}{1+\sqrt{1-\rho(D^{-1}(L+U))^{2}}}. $

So we have a very useful corollary.

Corollary 2 Let $ {A}$ be a tridiagonal matrix, such that $ {\rho(D^{-1}(L+U))<1}$ and $D^{-1}(L+U)$ has real eigenvalues. Then the SOR iteration converges if and only if $ {0<\omega<2}$.

MATH3310 Note 8: Iterative methods for linear systems, an introduction

The main problem of interest here is to solve a square linear system, which in matrix form reads

$ \displaystyle Ax=b,\quad A\in\mathbb{R}^{N\times N}. $

The type of algorithm that operating directly on the matrix $ {A}$ belongs to the family of direct methods. The simple Gaussian elimination (LU decomposition) is one of this kind, and can be used to solve such a system in $ {O(N^{3})}$ time, yielding an exact solution. But that is not always preferable, because it is typically time and power consuming when the matrix $ {A}$ is too large and no structural information of $ {A}$ is used in the algorithm. Various improvements of computation complexity in direct methods are possible, but will be a research topic. Here, we will take a look at some simple indirect methods, or iterative methods, which try to avoid computation with $ {A}$ as much as possible.

Generally speaking, one instead deals with the equivalent system

$ \displaystyle Mx=(M-A)x+b $

where $ {M}$ has to be designed based on $ {A}$. This can be seen as a splitting of $ {A}$

$ \displaystyle A=M+(A-M). $

And then we will take the fixed point iteration

$ \displaystyle Mx^{k+1}=(M-A)x^{k}+b $

since if $ {x^{*}=A^{-1}b}$, it is a fixed point of the above iteration

$ \displaystyle Mx^{*}=(M-A)x^{*}+b. $

Whether this iteration converges or not needs to be analysed. To do this, we consider the discrepancy of $ {x^{*}}$ from the $ {k}$-th iteration

$ \displaystyle \begin{array}{rcl} e^{k+1} & = & x^{*}-x^{k+1}\\ & = & M^{-1}(M-A)(x^{*}-x^{k})\\ & = & M^{-1}(M-A)e^{k}\\ & =: & Be^{k}. \end{array} $

This means the convergence is controled by the powers of the error matrix $ {B}$, and the iteration converges if and only if

$ \displaystyle B^{n}\rightarrow0\text{ as }n\rightarrow\infty. $

Exercise 1 Show that $ {B^{n}\rightarrow0}$ if and only if every eigenvalue $ {\lambda_{i}}$ of ${B}$ has modulus strictly less than one, i.e.

$ \displaystyle |\lambda_{i}|<1\quad\forall i. $

(Hint: use the Jordan normal form of $ {B}$.) The largest eigenvalue in modulus governs the rate of convergence, which is called the spectral radius, or $L^2$ operator norm of the matrix $ {B}$.

There are several criteria on what a "good'' $ {M}$ should be:

• $ {M}$ is easy to invert.
• $ {M-A}$ is easy to compute.
• The matrix $ {B:=M^{-1}(M-A)}$ has spectral radius stricly less than one.

Below the fold let us look at some simple concrete choices of $ {M}$. For this purpose we write

$ \displaystyle A=L+D+U $

where $ {L,D,U}$ are the strict lower triangular part, diagonal part and strict upper triangular part of $ {A}$ respectively.

1.1. The Jacobi iteration

This unfortunately is a rather uninformative name. We take

$ \displaystyle M=D. $

Hence in order for $ {M}$ to be invertible, the diagonal part of $ {A}$ must be non-zero. We will see that a sufficient condition for the iteration to converge is a condition on the diagonal part of $ {A}$.

Let us take a small concrete example to see what is going on. Take

$ \displaystyle A=\begin{bmatrix}2 & -1\\ -1 & 2 \end{bmatrix}. $

One Jacobi iteration looks like the following

$ \displaystyle \begin{bmatrix}2\\ & 2 \end{bmatrix}\begin{bmatrix}x_{0}^{1}\\ x_{1}^{1} \end{bmatrix}=\begin{bmatrix} & 1\\ 1 \end{bmatrix}\begin{bmatrix}x_{0}^{0}\\ x_{1}^{0} \end{bmatrix}+b. $

Note that one has to store the entire $ {x^{0}}$ in order to compute $ {x^{1}}$. This is what happens in general, since $ {M-A}$ is usually not triangular.

The error matrix is

$ \displaystyle B=D^{-1}(D-A), $

where we recognize $ {D-A}$ is the negative non-diagonal part of $ {A}$. So $ {B}$ is the negative non-diagonal part of $ {A}$, with each row divided by the corresponding diagonal entry, which was assumed to be non-zero. The larger the diagonal entries are, the closer to zero $ {B}$ becomes. To make this condition precise, we need the following useful fact.

Exercise 2 (Gershgorin disk) Let $ {A=(a_{ij})\in\mathbb{R}^{N\times N}}$. Then each eigenvalue of $ {A}$ (which can be complex) lies within at least one of the closed disks $ {\{z:|z-a_{ii}|\leq R_{i}\}}$ in $ {\mathbb{C}}$, where the radius $ {R_{i}}$ is defined by

$ \displaystyle R_{i}=\sum_{j\neq i}|a_{ij}|. $

Exercise 3 What are the Gershgorin disks of the matrix $ {B}$ in the Jacobi iteration? Hence show that Jacobi iteration is convergent if the matrix $ {A}$ is strictly diagonal dominant (SDD), i.e.

$ \displaystyle |a_{ii}|>\sum_{j\neq i}|a_{ij}|. $

Exercise 4 Show that a SDD matrix is invertible.

1.2. The Gauss-Seidel iteration

Again unfortunately this name is rather uninformative. I would rather call it method of successive displacement, for reason to be clear in a moment. We take

$ \displaystyle M=D+U\:\text{ or }\:M=D+L. $

Now since $ {M}$ is triangular, its eigenvalues lie on the diagonal. Hence in order for $ {M}$ to be invertible, the diagonal of $ {A}$ must be nonzero. The same SDD condition is also sufficient for the Gauss-Seidel to converge.

We take the same example to illustrate what is going on. With

$ \displaystyle A=\begin{bmatrix}2 & -1\\ -1 & 2 \end{bmatrix},\quad M=\begin{bmatrix}2\\ -1 & 2 \end{bmatrix} $

one Gauss-Seidel iteration looks like the following

$ \displaystyle \begin{bmatrix}2\\ -1 & 2 \end{bmatrix}\begin{bmatrix}x_{0}^{1}\\ x_{1}^{1} \end{bmatrix}=\begin{bmatrix} & 1\\ \ \end{bmatrix}\begin{bmatrix}x_{0}^{0}\\ x_{1}^{0} \end{bmatrix}+b. $

Note that to obtain $ {x_{1}^{1}}$, one needs $ {x_{0}^{1}}$ and $ {x_{1}^{0}}$ but not $ {x_{0}^{0}}$, due to our choice that $ {M}$ is the lower triangular part of $ {A}$, and cpnsequently $ {M-A}$ is strictly upper triangular. So, instead of storing the entire $ {x^{0}}$ before $ {x^{1}}$ is obtained, one can replace the entry $ {x_{0}^{0}}$ with the computed $ {x_{0}^{1}}$ immediately, hence saving half of the storage.

Exercise 5 Check the case $ {M=D+U}$ using the above $ {2\times2}$ example.

The SDD condition still guarantees the convergence here. Assume $ {M}$ is taken to be the lower triangular part (the other case is similar). Now consider the eigenvlaues $ {\lambda}$ of

$ \displaystyle B=M^{-1}(M-A)=-M^{-1}U. $

We have

$ \displaystyle \ker(-M^{-1}U-\lambda I)\neq\{0\}. $

Since $ {M}$ is invertible, this implies

$ \displaystyle \ker(U+\lambda M)\neq\{0\}. $

Note that $ {U+\lambda M=U+\lambda D+\lambda L}$. If $ {A}$ is SDD, then $ {U+\lambda D+\lambda L}$ is SDD only if $ {|\lambda|\geq1}$. But if $ {|\lambda|\geq1}$, $ {U+\lambda D+\lambda L}$ is invertible, by Exercise 4. This shows that eigenvalues of $ {B}$ in modulus cannot be greater than $ {1}$ if $ {A}$ is SDD. This guarantees the convergence by Exercise 1.


Another sufficient condition for convergence. Besides the SDD conditions, we also have

Theorem 1

1. The Gauss-Seidel iteration converges if $ {A}$ is symmetric and positive definite.
2. The Jacobi iteration converges if both $ {A}$ and $ {2D-A}$ are symmetric positive definite.

These are immediate consequences of

Theorem 2 (Householder-John) If $ {A,B\in\mathbb{R}^{N\times N}}$ and both $ {A}$ and $ {A-B-B^{T}}$ are symmetric positive definite, then the spectral radius of $ {(A-B)^{-1}B}$ is stricly less than one.

Proof: Let $ {\lambda}$,$ {v\neq0}$ be an eigenvalue-eigenvector pair of $ {-(A-B)^{-1}B}$, which may well be complex. Then

$ \displaystyle -Bv=\lambda(A-B)v. $

Hence $ {\lambda\neq1}$ since $ {A}$ is invertible. Multiplying both sides by $ {\bar{v}^{T}}$, we obtain

$ \displaystyle -\bar{v}^{T}Bv=\lambda\bar{v}^{T}Av+\lambda\bar{v}^{T}Bv. $

Thus

$ \displaystyle \bar{v}^{T}Bv=\frac{\lambda}{\lambda-1}\bar{v}^{T}Av. $

Now one easily check that the symmetric positive definiteness of a matrix $ {M}$ implies

$ \displaystyle \bar{x}^{T}Mx>0\quad,\forall x\in\mathbb{C}^{N}. $

Also,

$ \displaystyle \bar{v}^{T}B^{T}v=\left(\bar{v}^{T}Bv\right)^{*}=\frac{\bar{\lambda}}{\bar{\lambda}-1}\bar{v}^{T}A^{*}v=\frac{\bar{\lambda}}{\bar{\lambda}-1}\bar{v}^{T}Av. $

Apply the above observations to $ {M=A-B-B^{T}}$ and $ {M=A}$, with $ {x=v}$, we get

$ \displaystyle \begin{array}{rcl} 0 & < & \bar{v}^{T}(A-B-B^{T})v\\ & = & \bar{v}^{T}Av-\frac{\lambda}{\lambda-1}\bar{v}^{T}Av-\frac{\bar{\lambda}}{\bar{\lambda}-1}\bar{v}^{T}Av\\ & = & \frac{1-|\lambda|^{2}}{|1-\lambda|^{2}}\bar{v}^{T}Av \end{array} $

and

$ \displaystyle 0<\bar{v}^{T}Av. $

We deduce that $ {|\lambda|<1}$. $ \Box$

QR-Algorithm for Eigenvalue Problem

  I still remember Prof. Gilbert Strang in his linear algebra course once commented (which I learned from YouTube) that the Jordan canonical form of a matrix (see also my discussion here in a more general context) was once the culminating point in a linear algebra course long ago (50 years ago, perhaps). Indeed, the Jordan Canonical form $ {J}$ of a square matrix $ {A:V\rightarrow V}$ (non-sqaure matrix can be complete to sqaure matrix by adding zeros) gives us the precise information of its invariant subspaces, that is a linear subspace $ {W\subset V}$ satisfying

$ \displaystyle A\downharpoonright_{W}:W\rightarrow W. $

Abusing the notation, if we also write $ {W\in\mathbb{C}^{n\times r}}$ as a column arrangement of a basis of the subspace $ {W}$, then that the space $ {W}$ is invariant is characterized by

$ \displaystyle AW=W\Lambda $

for some $ {\Lambda\in\mathbb{C}^{r\times r}}$. Moreover, the Jordan block form tells us exactly how $ {A}$ acts on these invariant spaces, and once we know the actions on invariant subspaces, then the matrix is determined up to similarity transformations. One method to find the Jordan blocks is as follows: determine the eigenspaces and their bases (generalized eigenvectors), then put the basis vectors in columns in a particular order to form matrix $ {P}$ so that

$ \displaystyle P^{-1}AP=J. $

The theory side is clean. However, numerically to find such a decomposition is likely to be erraneous. After all, for matrices of very large size usually one cannot even evaluate the characteristic polynomial and obtain the eigenvectors directly. To remedy, there is a clever and primitve numerical method to obtain the eigenvectors, i.e. the power method, by iterating the sequence

$ \displaystyle x^{(m+1)}=Ax^{(m)}, $

where $ {A\in\mathbb{C}^{n\times n}}$ is our matrix and $ {x^{(m)}\in\mathbb{C}^{n}}$. With carefully chosen initial $ {x^{(0)}}$ and renormalization at each step, one can hope $ {x^{(m)}}$ will converge to an eigenvector that corresponds to the largest eigenvalue of $ {A}$. To do better, one can upgrage to

$ \displaystyle X^{(m+1)}=AX^{(m)}, $

where $ {X^{(m)}\in\mathbb{C}^{n\times r}}$ for some $ {1\leq r\leq n}$, and hope that after renormalization at each step, it should converge to a basis of a invariant subspace. So one can use this iterative method to approximate $ {P}$, in the form

$ \displaystyle \hat{P}^{-1}A\hat{P}=P_{k}^{-1}\cdots P_{1}^{-1}AP_{1}\cdots P_{k}. $

But the column vectors of $ {\hat{P}}$ may be ``too close'' to each other, and thus becomes ``numerically dependent''. This motivates us to re-orthonormalize $ {X^{(\cdot)}}$ at each iteration, using Gramm-Schmidt orthogonalization procedure. This has another cute name, called the QR factorization :

$ \displaystyle X^{(m)}=Q^{(m)}R^{(m)}, $

where $ {Q^{(m)}}$ is an isometric matrix (or unitary matrix if the numbers are complex), and $ {R^{(m)}}$ is upper triangular. Replacing $ {X^{(m)}}$ by $ {Q^{(m)}R^{(m)}}$, we get

$ \displaystyle X^{(m+1)}=AQ^{(m)}R^{(m)} $

Then taking the QR-factorization of $ {AQ^{(m)}=Q^{(m+1)}\hat{R}^{(m+1)}}$, we automatically get a QR factorization of $ {X^{(m+1)}=Q^{(m+1)}\hat{R}^{(m+1)}R^{(m)}=Q^{(m+1)}R^{(m+1)}}$. So the $ {R^{(m)}}$ term in fact becomes redundant in the iteration. Hence instead we can do

$ \displaystyle Y^{(m+1)} = AQ^{(m)}, Y^{(m+1)} = Q^{(m+1)}\hat{R}^{(m+1)}. (\star) $

The iteration should eventually take $ {Q^{(\cdot)}}$ to an orthonormal basis of its invariant subspace. In case $ {r=n}$ (i.e. $ {X^{(\cdot)}}$'s are $ {n\times n}$), this means there is an ismometric $ {Q}$ and upper triangular $ {\Lambda}$ such that

$ \displaystyle AQ=Q\Lambda. $

Multiplying $ {Q^{*}}$ from the right, we obtain the so-called Schur decomposition:

$ \displaystyle A=Q\Lambda Q^{*}. $

And since $ {\Lambda}$ is upper triangular, we get a nested sequence of invariant subspaces

$ \displaystyle V_{1}\subset V_{2}\subset\cdots V_{n}=\mathbb{C}^{n} $

and the 1st, 2nd, ..., and the $ {j}$-th columns spans the space $ {V_{j}}$, for each $ {j}$. Therefore, the iteration ($ {\star}$) can be considered as of getting an approximation of the Schur decomposition of $ {A}$ when $ {r=n}$. More precisely, we can take

$ \displaystyle A^{(m)}=Q^{(m)*}AQ^{(m)} $

to be an approximation of $ {\Lambda}$. So substituting this back into ($ {\star}$), we get

$ \displaystyle Q^{(m)*}Y^{(m+1)}=Q^{(m)*}Q^{(m+1)}\hat{R}^{(m+1)}=Q^{(m)*}AQ^{(m)}=A^{(m)}. $

Letting

$ \displaystyle Q^{(m)*}Q^{(m+1)}=\hat{Q}^{(m+1)}, $

we automatically get a QR-factorization of

$ \displaystyle A^{(m)}=\hat{Q}^{(m+1)}\hat{R}^{(m+1)}. $

Now, suppose we have found $ {A^{(m+1)}=Q^{(m+1)*}AQ^{(m+1)}}$. Then substituting $ {Q^{(m+1)}}$ by $ {Q^{(m)}\hat{Q}^{(m+1)}}$, we find

$ \displaystyle \begin{array}{rcl} A^{(m+1)} & = & \left(Q^{(m)}\hat{Q}^{(m+1)}\right)^{*}AQ^{(m)}\hat{Q}^{(m+1)}\\ & = & \hat{Q}^{(m+1)*}Q^{(m)*}AQ^{(m)}\hat{Q}^{(m+1)}\\ & = & \hat{Q}^{(m+1)*}A^{(m)}\hat{Q}^{(m+1)}\\ & = & \hat{Q}^{(m+1)*}\hat{Q}^{(m+1)}\hat{R}^{(m+1)}\hat{Q}^{(m+1)}=\hat{R}^{(m+1)}\hat{Q}^{(m+1)}. \end{array} $

Summerising, we obtained

$ \displaystyle \begin{array}{rcl} A^{(m)} & = & \hat{Q}^{(m+1)}\hat{R}^{(m+1)}\\ A^{(m+1)} & = & \hat{R}^{(m+1)}\hat{Q}^{(m+1)}. (\star \star) \end{array} $

This cute little algorithm is called QR iteration, standing at its elegance that the next iteration $ {A^{(m+1)}}$ is obtained by simply swapping the matrix multiplication.

Structure theorem for finitely generated modules over PIDs, Jordan canonical forms and Cayley Hamilton's theorem.

    Introduction    

The structure theorem, stated for finitely generated modules over PIDs, is a nice and elementary example of the usefulness of method of representation (in this case, that of a ring). Many important results in algebra follows immediately from this theorem: the structure theorem of finite abelian groups, rational canonical forms of matrices, also Jordan canonical forms, Cayley-Hamilton theorem. We will discuss the main theorem, with enough preparation until the proof seems trivial (thus omitted and refered to [A]: Artin., Algebra.) Several mentioned consequences will also be discussed. One of the slightly technical result is refered to a different exposition [C]: Conrad., Modules over a PID.