Matrix decompositions

Some decompositions

The algorithms for the QR and LU decompositions are fairly self-explanatory from their definition.

The Cholesky decomposition can be computed for a nonnegative-definite matrix through a simultaneous LU decomposition from both sides. Starting with a matrix in the block form:

\[A = \left[ {\begin{array}{*{20}{c}}
  {{a_1}}&{{v^T}} \\
  v&{{A_2}}
\end{array}} \right]\]
We perform an LU decomposition first on the columns, from the left, then on the rows, from the right:

\[\begin{gathered}
  \left[ {\begin{array}{*{20}{c}}
  {{a_1}}&{{v^T}} \\
  v&{{A_2}}
\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
  {\sqrt {{a_1}} }&0 \\
  {v/\sqrt {{a_1}} }&1
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
  {\sqrt {{a_1}} }&{{v^T}/\sqrt {{a_1}} } \\
  0&{{A_2} - v{v^T}/{a_1}}
\end{array}} \right] \\
   = \left[ {\begin{array}{*{20}{c}}
  {\sqrt {{a_1}} }&0 \\
  {v/\sqrt {{a_1}} }&1
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
  1&0 \\
  0&{{A_2} - v{v^T}/{a_1}}
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
  {\sqrt {{a_1}} }&{{v^T}/\sqrt a } \\
  0&1
\end{array}} \right] \\
\end{gathered} \]
The process can then be continued to give an LU decomposition where the $L$ and $U$ are transpose, thus the Cholesky decomposition.

(Why does the matrix have to be nonnegative-definite?)

Eigenvalue algorithm

A lot of factorisations are fundamentally related to the "eigen"-stuff: this includes the change-of-basis factorizations: diagonalization, SVD, Schur -- and also the polar decomposition. Eigenstuff is fundamentally analytic (whatever that means), so terminating algorithms won't work, we need an asymptotic algorithm.

A simple such algorithm to calculate the eigenvectors and eigenvalues is power iteration. It is related to the flows of systems of linear differential equations, in which the eigenspace with the largest eigenvalue provides the only stable equilibrium of the system. Unless you choose an initial vector $v$ that happens to itself be an eigenvector, $A^nv$ will approach the projection of $v$ onto the principal eigenspace of $A$.

So that gives you the largest eigenvalue. There are many ways to get the remaining eigenvalues -- one inefficient way I thought of was to subtract off the part corresponding to the stretching of the primary eigenvector, but this is incredibly inefficient and numerically unstable.

An actually sensible approach is as follows: consider the matrix $(A-\mu I)^{-1}$, which has the same eigenspaces as $A$. The largest eigenvalue of this is the eigenvalue of $A$ closest to $\mu$. By varying $\mu$ across the real line, one can discover all the eigenvalues of $A$.

Schur algorithm

You can use the eigenvalue algorithm to calculate the Schur decomposition directly via its definition -- however, there is a more efficient method.

A triangular matrix is invariant under conjugation by the $Q$ in its own $QR$ decomposition (because the $Q$ is the identity matrix). So similar to power iteration, one may iteratively QR-factorize $A$ and replace it with $RQ$. This is called the QR algorithm.

Apparently this is expensive in its crude form, so you instead bring into a form called "upper Hessenberg form" (upper-triangular but with a subdiagonal) and that apparently makes the computations less expensive. The mechanism to bring a matrix into upper Hessenberg form is based on householder reflections, and involves doing householder reflections on the parts of the columns below the subdiagonal, so the corresponding householder reflections on the rows (since the Schur is a change-of-basis transformation) do not interfere with the column of interest and screw up all your zeroes. It's not very interesting.

Least-squares

Suppose you wanted to solve $Ax\approx b$, i.e. find the $x$ that minimizes $\|Ax-b\|$. This is relevant when $A$ is not surjective.

The idea is that $\|Ax-b\|$ is minimised when $Ax$ is the projection $b_A$ of $b$ onto the image of $A$. One can equivalently write $A^{T}(Ax-b)=0$, as $Ax-b$ is perpendicular to the column space of $A$. Thus it suffices to solve $A^TAx=A^Tb$.

If $A$ is full-rank (this does not mean surjective), we can actually have a unique $x$, given by $x=(A^TA)^{-1}A^Tb$, and $A^+=(A^TA)^{-1}A^T$ is called the "Moore-Penrose pseudoinverse" of $A$.

An alternative, more efficient algorithm is to use the QR factorization (the one with the non-square $Q$) to construct the projector as $QQ^T$ (which works, as $Q$ definitionally provides a basis for the image, and $QQ^T$ satisfies the defining property of a Hermitian projector). Then one can solve $QRx=QQ^Tb$, or equivalently $Rx=Q^Tb$. In the full-rank case, $R^{-1}Q^T$ is the Moore-Penrose pseudoinverse.

To add:

• SVD algorithm
• Algorithm complexity for each thing

Triangular matrices: Schur, QR, Cholesky and LU/LPU decompositions

Upper-triangular matrices seem like they should be something interesting. They seem to be far more "controllable" in their behavior than general matrices -- the first dimension is just stretched, and each successive dimension's behavior does not depend on the higher dimensions. And there are of course practical computational reasons why they're simpler in a certain sense.

It's sensible to wonder if a general matrix can be understood in this way -- yes, yes, of course you can with the Jordan normal form, but if there's something more interesting to be said.

Some properties of triangular matrices:

• The determinant of a triangular matrix is the product of its diagonal values: this is a generalization of the fact about the area of a parallelogram being its base times perpendicular height.
• The eigenvalues of a triangular matrix are its diagonal values, including (algebraic) multiplicity: Follows from the above since $A-\lambda I$ remains triangular. 
• Normal triangular matrices are diagonal: The $k+1$th eigenvector must be in the orthogonal complement of $\mathbb{R}^k$ in $\mathbb{R}^{k+1}$, which is one-dimensional and just the $x_{k+1}$-axis. 

I tend to think that all linear transformations already "geometrically" look like triangular matrices, though -- like you should be able to tilt your head some way so that the transformation. What this means is that we believe all linear transformations are unitarily triangularizable.

Is this so?

We can certainly identify the first vector in the change-of-basis matrix -- an actual eigenvector of the matrix. This corresponds to the eigenvalue represented by the top-left element of the triangular matrix.

Now we want a vector whose image is a linear combination of this eigenvector and itself. Here's an idea: take the projection of the matrix onto the orthogonal complement of the first eigenvector. Then this projection itself has an eigenvector, and the action of the entire matrix on this vector is it scaled plus some component in the direction orthogonal to the subspace, i.e. proportional to the first eigenvector. We can repeat this process to obtain the Schur decomposition of $A$:

$$A=UTU^{*}$$
For unitary $U$, triangular $T$.

(By the way, a sequence of embedded subspaces of consecutively increasing dimension is known as a flag. For this reason, the Schur decomposition is said to say that every linear transformation stabilizes a flag.)

---

The Gram-Schmidt process might remind you of triangular matrices in the nature of the transformations applied.

$$\begin{array}{*{20}{c}}
  {{v_1} = {x_1}}&{{u_1} = \frac{{{v_1}}}{{\left| {{v_1}} \right|}}} \\
  {{v_2} = {x_2} - ({x_2} \cdot {u_1}){u_1}}&{{u_2} = \frac{{{v_2}}}{{\left| {{v_2}} \right|}}} \\
   \vdots & \vdots  \\
  {{v_k} = {x_k} - \sum {({x_i} \cdot {u_i}){u_i}} }&{{u_k} = \frac{{{v_k}}}{{\left| {{v_k}} \right|}}} \\
   \vdots & \vdots
\end{array}$$
Indeed, one may check that this is equivalent to writing $A=QR$ for unitary $Q$, triangular $R$. This factorization is called the QR-factorization.

But there's another way to think about this, right? What $A=QR$ says is that the transformation $A$ arises from first performing a triangular transformation, then rotating it into the desired orientation.

Imagine this in three dimensions, which I couldn't be bothered to draw.

And one can easily check that the algorithm for doing so is precisely the Gram-Schmidt algorithm.

But there's yet another way to see this factorization. Instead of just calculating the orthonormal vectors at each step, one could actually try to transform our column vectors into those of the triangular matrix through a bunch of orthogonal transformations.

One specific way of doing this is through reflections, specifically reflections in $n-1$-dimensional planes, known as Householder reflections.

1. The first reflection is in the plane midway between $a_1$ and $e_1$, and brings $a_1$ to the $e_1$ axis.
2. The next reflection is in the plane that contains the $e_1$ axis and is midway between $a_2$ and the $(e_1,e_2)$ plane, and brings $a_2$ onto the $(e_1,e_2)$ plane.
3. The next reflection is in the plane that contains the $(e_1,e_2)$ plane and is midway between $x_3$ and the $(e_1,e_2,e_3)$ plane, and brings $a_3$ onto the $(e_1, e_2, e_3)$ plane.
4. ...

I.e. in step $i$ (between 1 and $n$), all the vectors of $A$ are reflected in the plane that contains the $(e_1,\dots e_{i-1})$ plane and is midway between $a_i$ and the $(e_1,\dots e_{i})$ plane. The vector is then never reflected again, as it is in all the future planes of reflection.

The reflection matrix $Q_i$ can be seen to be of the form (do you see why?):

$${Q_i} = \left[ {\begin{array}{*{20}{c}}
  I&0 \\
  0&{{F_i}}
\end{array}} \right]$$
Where $F_i$ is $(n-i)$-dimensional.

Here's how we can calculate $F_i$: we know that $F_i$ transforms the latter $n-k$ elements of $a_i$ like this:

\[{a^{i + 1:n}_i} \mapsto \left[ {\begin{array}{*{20}{c}}
  {\left| {{a^{i + 1:n}_i}} \right|} \\
  0 \\
   \vdots  \\
  0
\end{array}} \right]\]
Although this generally isn't enough to pin down a transformation, we already know that the transformation is a reflection in an $n-1$-dimensional plane. We know that the reflection matrix can be given as $1-2hh^T$ where $h$ is the unit normal to the plane of reflection. We can calculate this normal (up to normalization) as $F_ia^{i+1:n}_i-a^{i+1:n}_i$.

The QR decomposition is basically just another "all matrices are basically ____" theorem -- similar to the polar decomposition "all matrices are basically positive-definite" -- they're just a rotation away.

In particular, this gives us a new square root of a positive-definite matrix $M=AA^T$ -- much like the polar decomposition gives a unique positive-definite square root of a positive-definite matrix, the QR decomposition gives a unique triangular square root of a positive-definite matrix $M=RR^T$.

This is known as the Cholesky decomposition.

(BTW, we can QR-decompose an $m$ by $n$ matrix ($m\ge n$) too. There are two ways to do this

• Thinking in the Gram-Schmidt sense: we're transforming into an orthogonal basis is the image of $A$, so $Q$ is an $m$ by $n$ matrix of orthonormal vectors and $R$ is an $n$ by $n$ matrix.
• Thinking in the Householder sense, we're obliged to end up with an $m$ by $m$ $Q$ and an $m$ by $n$ $R$ (with the bottom extra rows being zeroes). Then the extra columns of $Q$ are just some arbitrary basis for the orthogonal complement of the image of $A$ (i.e. the left null space of $A$).) 

---

We've all solved linear equations $Ax=b$ through elimination, but let's get a consistent general algorithm (Gaussian elimination) to do so.

The idea is to get a $A$ down to a triangular form, isn't it? And we can do that by zeroing each column one-by-one. That's LU decomposition. But sometimes you have a zero in the row you don't want to zero that column of, so you can't multiply it with infinity, so just swap two rows ("pivoting") and continue. That's LPU decomposition, where $P$ is a permutation matrix (the $P$ can come before or after the $LU$, too -- the fact that you can factor out the permutation matrix this way is left as an exercise).