High-dimensional filtering with permutohedral lattice

  High-dimensional filtering refers to convolution with kernel depending on more than the spatial variable. Examples include the bilateral filtering,

$ \displaystyle I'_{i}=\sum_{j=1}^{N}k(f_{i}-f_{j})I_{j}, $

where the kernel $ {k(\cdot)}$ depends on a featrure variable including both the position and the image intensity

$ \displaystyle f_{i}=(p_{i},I_{i}); $

with a similar form there is also the more complicated non-local mean filtering. In a 2009 paper (referred to as the paper below), Adams, Baek and Davis proposed a fast, novel algorithm based on permutohedral lattice, whose memory usage also scales well with dimension. This was later utilized in the mean field inference of the fully connected conditional random field model of Krähenbühl and Koltun. In this post I wish to describe some details of the permutohedral lattice approach.

Generally speaking, there are two type of filters:

• High pass filters, which keep the high frequency information;
• Low pass filters, which keep the low frequency information.

The filter is usually locally supported in a window of fixed size, after applying which it calculates the local weighted average in that window. Conventional convolutions on a equally spaced rectangular grid can be represented as (block-)circulant matrices, where the weight is the same in each window thanks to the fact that filter only depends on the spatial variable as well as the regular grid structure. This is no longer true for high-dimensional filtering, or an unstructured grid that is the case for a geometric domain. Since filters now depends on additional variables other than the spatial variable. In bilateral filtering the filter depends on the intensities of the pixels inside the window. In a naive implementation filters need to be re-computed for each new window, yielding quadratic complexity in filter's radius. The situation is similar for an unstructured grid.

The approach here is to embed the features into a high-dimensional space (feature space) and interpolate the input signals to the regular grid point in the high-dimensional space. Efficient local averaging can then be done on the grid. And finally the values of the filtering are interpolated back to the original input feature position. These procedures are usually described vividly as splatting - blurring - slicing.

High-dimension filter with permutohedral lattice. Figure taken from the paper.

The choice of the regular grid affects the efficiency of this approach. If the standard grid $ {\mathbb{Z}^{d}}$ is chosen, say $ {d}$ is the dimension of the feature, then each input signal needs to be "splatted'' onto $ {2^{d}}$ nearby grid points. This is exponential complexity in dimension. In contrast, for the permutohedral lattice a input signal is "splatted'' onto only $ {d+1}$ nearby grid points, thus reducing the complexity to linear. Note that $ {d+1}$ is the number of vertices in a $ {d}$-dimensional simplex, meaning that this number cannot be improved further.

In the following we shall introduce the permutohedral lattice and layout some properties that will be useful for understanding the algorithm.

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}$.