%% ================================================================
\section{Linear algebra}

In this paper I compare and contrast two techniques for computation of
determinants and inverses of square matrices: the more-familiar
Gaussian-elimination method, and the less-familiar Householder method.  I
primarily make use of geometric methods to do so.  This requires some
preliminaries from linear algebra, including geometric interpretations of
determinant, matrix norm, and error propagation.

%% ----------------------------------------------------------------
\subsection{Geometric meanings of determinant and matrix norm}
\label{sec:geom_det_and_norm}

Let $A: \R^n \to \R^n$ be a linear transformation (which I will always think of
as a matrix with respect to the standard basis).  The \bfidx{determinant} of
$A$ can be thought of as the \bfidx{signed volume} of the image of the
\bfidx{unit cube}:

\includegraphics{det1.eps}

The \bfidx{matrix norm}\index{norm} is, by definition, the \bfidx{maximum
extent}\index{extent} of the image of the \bfidx{unit ball}:

\includegraphics{norm1.eps}

\begin{defn}
\label{defn:matrix_norm}
We define the matrix norm of $A$ either by
$$
	\|A\| = \sup_{\|x\|=1} \|Ax\|
$$
or equivalently by
$$
	\|A\| = \sup_{\|x\|\ne 0} \frac{\|Ax\|}{\|x\|}.
$$
That is, either we consider vectors in the unit ball and find the maximum
extent of their images, or we consider \emph{all} non-zero vectors and find the
maximum amount by which they are scaled.  The equivalence of these two
characterizations is due to the linearity of $A$.
\end{defn}

A matrix with determinant $\pm 1$ preserves (unsigned) volume, but does not
necessarily preserve norm.  Of particular interest for this paper are three
kinds of matrices:

A $2 \times 2$ \bfidx{rotation matrix} is of the form
$$
	A=\left[\begin{array}{rr}
	\cos(t) & -\sin(t) \\
	\sin(t) & \cos(t)
	\end{array}\right],
$$
and has determinant 1:

\includegraphics{det3.eps}

An example of a $2 \times 2$ \bfidx{reflection matrix}, reflecting
about the $y$ axis, is
$$
	A=\left[\begin{array}{rr}
	-1 & 0 \\
	0  & 1 \\
	\end{array}\right],
$$
which has determinant $-1$:

\includegraphics{det4.eps}

Another example of a reflection is a \bfidx{permutation matrix}:
$$
	A=\left[\begin{array}{rr}
	0 & 1 \\
	1 & 0 \\
	\end{array}\right],
$$
which has determinant $-1$:

\includegraphics{det5.eps}

This reflection is about the $45^\circ$ line $x=y$.  (Construction of a
reflection matrix about an arbitrary axis is accomplished using Householder
transformations, as discussed in section \ref{sec:hh}.)

An example of a $2 \times 2$ \bfidx{shear matrix} is
$$
	A=\left[\begin{array}{rr}
	1 & a \\
	0 & 1 \\
	\end{array}\right],
$$
which has determinant $1$:

\includegraphics{det2.eps}

%% ----------------------------------------------------------------
\subsection{Computation of determinants}
\label{sec:comp_det}

In elementary linear algebra (see perhaps \cite{bib:fis}), we are first taught
to compute determinants\index{determinant} using \bfidx{cofactor expansion}.
This is fine for computing determinants of $2 \times 2$'s or $3 \times 3$'s.
However, it is ineffective for larger matrices.

\begin{itemize}
\item The determinant of a $2 \times 2$ matrix
$$
	\begin{bmatrix} a & b \\ c & d \\ \end{bmatrix}
$$
is $ad-bc$: two multiplies and an add.

\item To compute the determinant of a $3 \times 3$ using cofactor expansion, we
work down a row or column.  There are 3 determinants of $2 \times 2$'s.

\item To do a $4 \times 4$, there are 4 determinants of $3 \times 3$'s, each of
which takes (recursively) 3 determinants of $2 \times 2$'s.

\item Continuing this pattern, we see that there are on the order of $n!$
multiplies and adds for an $n \times n$.  For example, $25! \approx 10^{25}$.
Even at a billion operations per second, this requires $10^{16}$ seconds,
which is hundreds of millions of years.  We can do better.

\end{itemize}

\begin{center} * * * \end{center}

If a matrix is \bfidx{upper-triangular} or \bfidx{lower-triangular}, then its
determinant is the product of its diagonal entries.  There are $n$ of these and
we only need to multiply them, so there are $n$ multiplies.

The question is, how efficiently can we \emph{put} a given (square) matrix into
upper-triangular form, and how does this modification affect the determinant?

Upper-triangularization is the process of putting zeroes in certain elements of
a matrix, while modifying other entries.  Recall that when a matrix $Q$ acts by
premultiplication on a matrix $A$, we can think of $Q$ acting on each
\bfidx{column vector} of $A$.  That is, the $j$th column of $QA$ is simply $Q$
times the $j$th column of $A$.  And certainly we can transform column vectors
to put zeroes in various locations.

How does this affect the determinant?  If we apply $n$ transformations,
$M_1$ through $M_m$, then
$$
	\det(M_m \cdots M_2 M_1 A) = \det(M_m) \cdots \det(M_2) \det(M_1) \det(A)
$$
i.e.
$$
	\det(A) = \frac{\det(M_m \cdots M_2 M_1 A) }
		{\det(M_m) \cdots \det(M_2) \det(M_1)
	}
	= \frac{\textrm{product along diagonal of upper-triangular matrix}}
		{\det(M_m) \cdots \det(M_2) \det(M_1)
	}.
$$
That is, all we need to do is keep track of the determinants of our
transformations, multiply them up, and compute the determinant of the resulting
upper-triangular matrix.

How long does it take to compute and apply all these $M_i$'s, though?  It can
be shown (but I will not show in this paper) that the Gaussian-elimination and
Householder methods for upper-triangularization are on the order of $n^3$.
Thus, these methods are far more efficient than naive cofactor expansion.

%% ----------------------------------------------------------------
\subsection{Computation of matrix inverses}
\label{sec:comp_inv}

In elementary linear algebra, we are taught to compute \plainidx{inverses}
using \bfidx{cofactor expansion}.  This also can be shown to require on the
order of $n!$ operations.

\begin{center} * * * \end{center}

A more efficient method, which we are also taught in elementary linear
algebra, is to use an \bfidx{augmented matrix}.  That is, we paste our
$n \times n$
input matrix $A$ next to an $n \times n$ identity matrix:
$$
	[ \; A \; | \; I \; ]
$$
and put the augmented matrix into reduced row echelon form.  If it is
possible to get the identity matrix on the left-hand side, then the inverse
will be found in the right-hand side:
$$
	[ \; I \; | \; A^{-1} \; ].
$$

\begin{prop}
This works.
\begin{proof}
When we row-reduce the augmented matrix, we are applying a sequence $M_1,
\ldots, M_m$ of linear transformations to the augmented matrix.  Let their
product be $M$:
$$
	M = M_n \cdots M_1.
$$
Then, if the row-reduction works, i.e. if we get $I$ on the left-hand side,
then
$$
	M[\;A\;|\;I\;] = [\;I\;|\;*\;],
$$
i.e. $MA = I$.  But this forces $M=A^{-1}$.  Thus
$$
	* = MI = A^{-1}I = A^{-1}.
$$
\end{proof}
\end{prop}

Note that putting a matrix into reduced row echelon form is a two-step
process:
\begin{itemize}
\item Putting it into row echelon form.  This simply means putting it into
\bfidx{upper-triangular} form and diving each row by its leading coefficient:
$$
	\begin{bmatrix}
		* & * & * & * & * & * \\
		* & * & * & * & * & * \\
		* & * & * & * & * & * \\
	\end{bmatrix}
	\mapsto
	\begin{bmatrix}
		* & * & * & * & * & * \\
		0 & * & * & * & * & * \\
		0 & 0 & * & * & * & * \\
	\end{bmatrix}
	\mapsto
	\begin{bmatrix}
		1 & * & * & * & * & * \\
		0 & 1 & * & * & * & * \\
		0 & 0 & 1 & * & * & * \\
	\end{bmatrix}.
$$
\item Clearing the above-diagonal entries:
$$
	\begin{bmatrix}
		1 & * & * & * & * & * \\
		0 & 1 & * & * & * & * \\
		0 & 0 & 1 & * & * & * \\
	\end{bmatrix}
	\mapsto
	\begin{bmatrix}
		1 & 0 & 0 & * & * & * \\
		0 & 1 & 0 & * & * & * \\
		0 & 0 & 1 & * & * & * \\
	\end{bmatrix}.
$$
\end{itemize}

The remaining question (to be addressed below) is what particular sequence
of $M_i$'s to apply.

%% ----------------------------------------------------------------
\subsection{Error propagation}
\label{sec:err_prop}

Whenever we use numerical techniques, we need to ask about error.  A poorly
designed algorithm can permit error to grow ridiculously, producing nonsense
answers.  This paper is about well-designed algorithms.  The study of error is
\bfidx{numerical analysis} and is discussed, for example, in \cite{bib:bf} and
\cite{bib:nr}.  Here I simply make a few points.

Consider a true matrix $A$ at the $k$th step of an upper-triangularization
process, along with an actual matrix $A + \veps$ where $\veps$ is a matrix.
We want to know how error \bfidx{propagates}\index{error propagation} when
we apply a transformation $Q$.
Since we use linear transformations, we have
$$
	Q(A+\veps) = Q(A) + Q(\veps).
$$
Also, since the $j$th column of $QA$ is $Q$ times the $j$th column of $A$,
we have
$$
	Q(\vecv+\eta) = Q(\vecv) + Q(\eta)
$$
where $\vecv$ and $\eta$ are column vectors.  The \bfidx{vector
norm}\index{norm} of $Q(\eta)$ is related to the vector norm of $\eta$ by the
\bfidx{matrix norm} of $Q$:  since (from definition \ref{defn:matrix_norm})
$$
	\|Q\| \ge \frac{\|Q(\eta)\|}{\|\eta\|},
$$
we have
$$
	\|Q(\eta)\| \le \|Q\| \; \|\eta\|.
$$
In particular, if a transformation matrix $Q$ is norm-preserving, it
does not increase error.

Next we consider the norms of the matrices which are used in the
Gaussian-elimination and Householder methods.

%% ================================================================
\newpage
\section{Gaussian elimination}

In this section, Gaussian elimination is discussed from the viewpoint of
concatenated linear transformations, with an emphasis on geometrical
interpretation of those transformations.

%% ----------------------------------------------------------------
\subsection{Row reduction using Gaussian elimination}

For clarity, consider matrices of height 2.  Let's row-reduce an
example matrix:
$$
	\left[\begin{array}{rr} 0 & 3 \\ 1 & 2 \end{array} \right] \mapsto
	\left[\begin{array}{rr} 1 & 2 \\ 0 & 3 \end{array} \right] \mapsto
	\left[\begin{array}{rr} 1 & 2 \\ 0 & 1 \end{array} \right] \mapsto
	\left[\begin{array}{rr} 1 & 0 \\ 0 & 1 \end{array} \right].
$$

This example illustrates the three matrices which we use in Gaussian
elimination:

The \bfidx{row-swap matrix}
(a \plainidx{permutation matrix} from section \ref{sec:geom_det_and_norm})
has determinant $-1$ and norm 1:
$$
	\begin{bmatrix}
		0 & 1 \\
		1 & 0 \\
	\end{bmatrix}
$$
The \bfidx{row-scaling matrix} has determinant $m$ (in the example,
$m=1/3$) and norm $\max(1,m)$:
$$
	\begin{bmatrix}
		1 & 0 \\
		0 & m \\
	\end{bmatrix}
$$
The \bfidx{row-update matrix}
(a \plainidx{shear matrix} from section \ref{sec:geom_det_and_norm})
has determinant $1$ and a norm which we can most
easily visualize using the diagrams in section \ref{sec:geom_det_and_norm}:
$$
	\begin{bmatrix}
		1 & a \\
		0 & 1 \\
	\end{bmatrix}.
$$
(In the example, $a=-2$.)

It's clear that Gaussian elimination magnifies error when a row-scaling matrix
has large $|m|$, and/or when a row-update matrix has a large $|a|$.

\includegraphics{badshear.eps}

% [ 1 -2 ] [ 1 2 ] = [ 1 0 ]
% [ 0  1 ] [ 0 1 ]   [ 0 1 ]

%% ----------------------------------------------------------------
\subsection{Gaussian elimination without pivoting}

Here is an example of the error which can accumulate when we naively use
Gaussian elimination.  First let's solve the linear system
\begin{eqnarray*}
	y &=& 1 \\
	x-y &=& 0
\end{eqnarray*}
which is, in matrix form,
$$
	\left[\begin{array}{rr}
		0 & 1 \\
		1 & -1
	\end{array}\right]
	\left[\begin{array}{r}
		x \\ y
	\end{array}\right]
	=
	\left[\begin{array}{r}
		1 \\ 0
	\end{array}\right]
$$
or in augmented form
$$
	\left[\begin{array}{rr|r}
		0 &  1 & 1\\
		1 & -1 & 0
	\end{array}\right].
$$
Clearly, the solution is
$$
	\left[\begin{array}{r}
		x \\ y
	\end{array}\right]
	=
	\left[\begin{array}{r}
		1 \\ 1
	\end{array}\right].
$$

Now let's perturb it a bit, solving the system
$$
	\left[\begin{array}{rr|r}
		0.001 &  1 & 1\\
		1 & -1 & 0
	\end{array}\right].
$$
We expect a solution close to $(1,1)$.  Here is Gaussian elimination, rounded
to 4 and 3 decimal places, respectively, with no pivoting operations.  (3
decimal places means that, in particular, 1001 rounds to 1000.) At each step I
include the transformation matrix which takes us to the next step, along with
the norm of that transformation matrix.
$$
\begin{small}
\begin{array}{lrr|lrr}
&\textrm{4 places}&&
&\textrm{3 places}& \\
\hline &&&&&\\
	\textrm{(norm 1000)}
	&
	\left[\begin{array}{rrr}
		1000 & 0 \\
		0    & 1 \\
	\end{array}\right]
	&
	\left[\begin{array}{rr|r}
		0.001 &  1 & 1\\
		1 & -1 & 0\\
	\end{array}\right]
	&
	%%%%%
	\textrm{(norm 1000)}
	&
	\left[\begin{array}{rrr}
		1000 & 0 \\
		0    & 1 \\
	\end{array}\right]
	&
	\left[\begin{array}{rr|r}
		0.001 &  1 & 1\\
		1 & -1 & 0\\
	\end{array}\right]
\\&&&&&\\
%%%%%%%
	\textrm{(norm $\approx 1$)}
	&
	\left[\begin{array}{rrr}
		1   & 0 \\
		-1  & 1 \\
	\end{array}\right]
	&
	\left[\begin{array}{rr|r}
		1 &  1000 & 1000\\
		1 & -1 & 0\\
	\end{array}\right]
	&
	%%%%%
	\textrm{(norm $\approx 1$)}
	&
	\left[\begin{array}{rrr}
		1 & 0 \\
		-1 & 1 \\
	\end{array}\right]
	&
	\left[\begin{array}{rr|r}
		1 &  1000 & 1000\\
		1 & -1 & 0\\
	\end{array}\right]
\\&&&&&\\
%%%%%%%
	\textrm{(norm $= 1$)}
	&
	\left[\begin{array}{rrr}
		1   & 0 \\
		0  & -1/1001 \\
	\end{array}\right]
	&
	\left[\begin{array}{rr|r}
		1 &  1000 & 1000\\
		0 & -1001 & -1000\\
	\end{array}\right]
	&
	%%%%%
	\textrm{(norm $= 1$)}
	&
	\left[\begin{array}{rrr}
		1   & 0 \\
		0  & -1/1000 \\
	\end{array}\right]
	&
	\left[\begin{array}{rr|r}
		1 &  1000 & 1000\\
		0 & -1000 & -1000\\
	\end{array}\right]
\\&&&&&\\
%%%%%%%
	\textrm{(norm $= 1000$)}
	&
	\left[\begin{array}{rrr}
		1  & -1000 \\
		0  & 1\\
	\end{array}\right]
	&
	\left[\begin{array}{rr|r}
		1 &  1000 & 1000\\
		0 & 1 & 0.999\\
	\end{array}\right]
	&
	%%%%%
	\textrm{(norm $= 1000$)}
	&
	\left[\begin{array}{rrr}
		1  & -1000 \\
		0  & 1\\
	\end{array}\right]
	&
	\left[\begin{array}{rr|r}
		1 &  1000 & 1000\\
		0 & 1 & 1 \\
	\end{array}\right]
\\&&&&&\\
%%%%%%%
	&
	&
	\left[\begin{array}{rr|r}
		1 & 0 & 0.999\\
		0 & 1 & 0.999\\
	\end{array}\right]
	&
	%%%%%
	&
	&
	\left[\begin{array}{rr|r}
		1 & 0 & 0\\
		0 & 1 & 1\\
	\end{array}\right]
\\&&&&&
%%%%%%%
\end{array}
\end{small}
$$

Here, the roundoff error dominates:  the right-hand solution is very wrong.
You can check that the left-hand solution is nearly correct.  There were two
norm-1000 operations, contributing to the error.

%
%A=[0.001,1,1; 1,-1,0];
%M1=[1000,0;0,1]
%M2=[1,0;-1,1];
%M3=[1,0;0,-1/1001];
%M4=[1,-1000;0,1];
%
%A
%M1*A
%M2*M1*A
%M3*M2*M1*A
%M4*M3*M2*M1*A
%
%fprintf(1, '*************\n')
%


%% ----------------------------------------------------------------
\subsection{Gaussian elimination with pivoting}

Now let's do the same example, but with pivoting.  This means that
we use a permutation matrix to place the largest absolute-value column
head at the top of a column.

Again using 4 and 3 decimal places, respectively, we have the following.
$$
\begin{array}{lrr|lrr}
&\textrm{4 places}&&
&\textrm{3 places}& \\
\hline &&&&&\\
	\textrm{(norm 1)}
	&
	\left[\begin{array}{rrr}
		0 & 1 \\
		1 & 0 \\
	\end{array}\right]
	&
	\left[\begin{array}{rr|r}
		0.001 &  1 & 1\\
		1 & -1 & 0\\
	\end{array}\right]
	&
	%%%%%
	\textrm{(norm 1)}
	&
	\left[\begin{array}{rrr}
		0 & 1 \\
		1 & 0 \\
	\end{array}\right]
	&
	\left[\begin{array}{rr|r}
		0.001 &  1 & 1\\
		1 & -1 & 0\\
	\end{array}\right]
\\&&&&&\\
%%%%%%%
	\textrm{(norm $\approx 1$)}
	&
	\left[\begin{array}{rrr}
		1 & 0 \\
		-1/1000 & 1 \\
	\end{array}\right]
	&
	\left[\begin{array}{rr|r}
		1 & -1 & 0\\
		0.001 &  1 & 1\\
	\end{array}\right]
	&
	%%%%%
	\textrm{(norm $\approx 1$)}
	&
	\left[\begin{array}{rrr}
		1 & 0 \\
		-1/1000 & 1 \\
	\end{array}\right]
	&
	\left[\begin{array}{rr|r}
		1 & -1 & 0\\
		0.001 &  1 & 1\\
	\end{array}\right]
\\&&&&&\\
%%%%%%%
	\textrm{(norm $\approx 1$)}
	&
	\left[\begin{array}{rrr}
		1 & 1/1.001 \\
		0 & 1 \\
	\end{array}\right]
	&
	\left[\begin{array}{rr|r}
		1 & -1 & 0\\
		0 &  1.001 & 1\\
	\end{array}\right]
	&
	%%%%%
	\textrm{(norm $\approx 1$)}
	&
	\left[\begin{array}{rrr}
		1 & 1 \\
		0 & 1 \\
	\end{array}\right]
	&
	\left[\begin{array}{rr|r}
		1 & -1 & 0\\
		0 &  1 & 1\\
	\end{array}\right]
\\&&&&&\\
%%%%%%%
	\textrm{(norm $\approx 1$)}
	&
	\left[\begin{array}{rrr}
		1 & 0 \\
		0 & 1/1.001 \\
	\end{array}\right]
	&
	\left[\begin{array}{rr|r}
		1 & 0 & 0.999\\
		0 &  1.001 & 1\\
	\end{array}\right]
	&
	%%%%%
	&
	&
	\left[\begin{array}{rr|r}
		1 & 0 & 1\\
		0 & 1 & 1\\
	\end{array}\right]
\\&&&&&\\
%%%%%%%
	&
	&
	\left[\begin{array}{rr|r}
		1 & 0 & 0.999\\
		0 &  1& 0.999\\
	\end{array}\right]
	&
	%%%%%
	&
	&
\\&&&&&
%%%%%%%
\end{array}
$$

Here, both solutions are approximately equal and approximately correct.  All
the transformations have norm on the order of 1, which in turn is the case
because of the pivoting operation.

The point is that successful Gaussian elimination \emph{requires} pivoting,
which is a \bfidx{data-dependent decision}.  Also, we must store the
\bfidx{permutations}.  Neither of these is a serious detriment in a carefully
implemented software system.  We will next look at Householder transformations
not as a \emph{necessity}, but as an \emph{alternative}.  It will turn out that
they will permit somewhat simpler software implementations.

%% ================================================================
\newpage
\section{Householder transformations}
\label{sec:hh}

In this section the Householder transformation is derived, then put to use in
applications.

%% ----------------------------------------------------------------
\subsection{Geometric construction}
\label{sec:geometric_construction}

Here we construct a \bfidx{reflection matrix} $Q$ acting on $\R^n$ which sends
a chosen \bfidx{axis vector}, $\veca$, to its negative, and reflects all other
vectors through the hyperplane perpendicular to $\veca$:

\includegraphics{hh1.eps}

How do we do this?  We can \plainidx{decompose} a given matrix $\vecu$ into
its \plainidx{components} which are \plainidx{parallel} and \plainidx{perpendicular} to $\veca$:

\includegraphics{projpar.eps}

Then, if we can subtract off from $\vecu$ twice its parallel component
$\vecu_{||}$, then we'll have obtained the reflection:

\includegraphics{projpar2.eps}

That is, we want
$$
	Q \vecu = \vecu - 2 \vecu_{||}.
$$
Moreover, we would like a single matrix $Q$ which would transform all
$\vecu$'s.

What is the projection vector $\vecu_{||}$?  It must be in the
\bfidx{direction} of $\veca$, and it must have \bfidx{magnitude}
$\|\vecu\|\cos\theta$:
$$
	\vecu_{||} = \ahat \|\vecu\| \cos\theta
		= \frac{\veca}{\|\veca\|} \|\vecu\| \cos\theta.
$$
I use the notation
$$
	\ahat = \frac{\veca}{\|\veca\|}.
$$
It turns out (this is the central trick of the Householder transformation) that
we can reformulate this expression in terms of dot products, eliminating
$\cos\theta$.  To do this, recall that
$$
	\veca \cdot \vecu = \|\veca\| \, \|\vecu\| \cos\theta.
$$
So,
$$
	\vecu_{||} = \frac{\veca}{\|\veca\|^2} \|\veca\|\,\|\vecu\| \cos\theta
	= \frac{\veca}{\|\veca\|^2} \veca \cdot \vecu.
$$
Also recall that
$$
	\|\veca\|^2 = \veca \cdot \veca \qqand
	\|\vecu\|^2 = \vecu \cdot \vecu.
$$
Now we have
$$
	\vecu_{||} = \veca \frac{\veca \cdot \vecu}{\veca \cdot \veca}.
$$

Now we can use this to rewrite the reflection of $\vecu$:
\begin{eqnarray*}
	Q\vecu &=& \vecu - 2 \vecu_{||} \\
	&=& \vecu - 2 \veca \frac{\veca\cdot\vecu}{\veca\cdot\veca}.
\end{eqnarray*}
Now (here's the other trick of Householder) remember that we can write
the \plainidx{dot product}, or \bfidx{inner product}, as a matrix product:
$$
	\veca \cdot \veca = \veca^t \veca.
$$
Here, we think of a \plainidx{column vector} as an $n \times 1$ matrix, while
we think of a \plainidx{row vector} as a $1 \times n$ matrix which is the
\plainidx{transpose} of the column vector.  So,
\begin{eqnarray*}
	Q\vecu &=& \vecu - 2 \veca \frac{\veca^t\vecu}{\veca^t\veca}.
\end{eqnarray*}
But now that we've recast the dot products in terms of
matrix multiplication, we can use the \plainidx{associativity}
of matrix multiplication:
$$
	\veca (\veca^t \vecu) = (\veca \veca^t) \vecu.
$$
The product $\veca\veca^t$ is the \bfidx{outer product} of $\veca$
with itself.  It is an $n \times n$ matrix with $ij$th entry $a_ia_j$.
Now we have
\begin{eqnarray*}
	Q\vecu &=& \left(I - 2 \frac{\veca \veca^t}{\veca^t \veca}\right) \vecu
\end{eqnarray*}
from which
\begin{eqnarray*}
	Q &=& I - 2 \frac{\veca \veca^t}{\veca^t \veca}.
\end{eqnarray*}

%% ----------------------------------------------------------------
\subsection{Construction with specified source and destination}
\label{sec:Q_spec_uv}

In section \ref{sec:geometric_construction}, we saw how to construct a
Householder matrix when given an axis of reflection $\veca$.  This would send
any other vector $\vecu$ to its image $\vecv$ through the hyperplane
perpendicular to $\veca$.

However, what if we were given $\vecu$ and $\vecv$, and wanted to find the axis
of reflection $\veca$ through which $\vecv$ is $\vecu$'s mirror image?  First
note that since Householder transformations are norm-preserving, $\vecu$ and
$\vecv$ must have the same norm.  (Normalize first them if not.) Drawing a
picture, and using the head-to-tail method, it looks like $\veca=\vecu-\vecv$
ought to work:

\includegraphics{uva.eps}

We can check that this is correct.  If $\veca$ is indeed the reflection axis,
then it ought to get reflected to its own negative:
\begin{eqnarray*}
	Q(\veca) &=& Q(\vecu - \vecv) \\
	&=& Q(\vecu) - Q(\vecv) \\
	&=& \vecv - \vecu \\
	&=& -\veca.
\end{eqnarray*}

%% ----------------------------------------------------------------
\subsection{Properties of $Q$, obtained algebraically}
\label{sec:Q_props_alg}

Here we attempt to prove algebraically that the Householder matrix
has certain properties.  In the next section, we do the same proofs
geometrically.

\begin{itemize}
\item $Q$ is \bfidx{orthogonal}, i.e. $Q^t Q = I$:
Compose $(I - 2\veca\veca^t/\veca^t\veca)$ with itself and FOIL it out.

\item $Q$ is \bfidx{symmetric}, i.e. $Q = Q^t$.  This is easy since
$Q=I-2\veca\veca^t/\veca^t\veca$.  The identity is symmetric; $\veca\veca^t$
has $ij$th entry $a_i a_j = a_j a_i$ and so is symmetric as well.

\item $Q$ is \bfidx{involutary}, i.e. $Q^2=I$: Same as orthogonality,
due to symmetry, since $Q=Q^t$.

\item \bfidxtwo{Determinant}{determinant}:  I don't see a nice way to find this
algebraically.  It seems like it would be a mess.

\item \bfidxtwo{Matrix norm}{matrix norm}:  Likewise.
\end{itemize}

%% ----------------------------------------------------------------
\subsection{Properties of $Q$, obtained geometrically}
\label{sec:Q_props_geom}

In this section, we prove certain properties of the Householder matrix using
geometric methods.  Remember that dot products, via $\vecu \cdot \vecv =
\|\vecu\|\,\|\vecv\|\,\cos\theta$, give us norms as well as angles.

\begin{itemize}
\item $Q$ is \bfidx{involutary}, i.e. $Q^2=I$:  $Q$ reflects $\vecu$ to
its mirror image $\vecv$; a second application of $Q$ sends it back again.

\item $Q$ is \bfidx{orthogonal}, i.e. $(Q\vecu \cdot Q \vecv) = \vecu \cdot
\vecv$ for all vectors $\vecu, \vecv$: Since $Q$ is a reflection, it preserves
norms; also, from the picture, it's clear that it preserves angles:

\includegraphics{orthog.eps}

\item $Q$ is \bfidx{symmetric}, i.e. $(Q\vecu) \cdot \vecv =
\vecu \cdot (Q\vecv)$ for all vectors $\vecu, \vecv$:  This is the
same geometric argument as for orthogonality, making use of involutarity.

\item \bfidxtwo{Determinant}{determinant}:  Since $Q$ is a reflection about the
$\veca$ axis, leaving all the axes orthogonal to $\veca$ fixed, $Q$ must have
determinant -1.  That is, it turns the unit cube inside out along one axis.

\item \bfidxtwo{Matrix norm}{matrix norm}:  Since $Q$ is a reflection, it is
length-preserving.
\end{itemize}

%% ----------------------------------------------------------------
\subsection{Repeated Householders for upper-triangularization}
\label{sec:HHUT}

The goal of upper triangularization is to put a matrix of the form
$$
	\begin{bmatrix}
		* & * & * & * & * & *\\
		* & * & * & * & * & *\\
		* & * & * & * & * & *\\
	\end{bmatrix}
$$
into the form
$$
	\begin{bmatrix}
		* & * & * & * & * & *\\
		0 & * & * & * & * & *\\
		0 & 0 & * & * & * & *\\
	\end{bmatrix}.
$$
We can do this one step at a time:  from
$$
	\begin{bmatrix}
		* & * & * & * & * & *\\
		* & * & * & * & * & *\\
		* & * & * & * & * & *\\
	\end{bmatrix}
$$
to
$$
	\begin{bmatrix}
		* & * & * & * & * & *\\
		0 & * & * & * & * & *\\
		0 & * & * & * & * & *\\
	\end{bmatrix}
$$
to
$$
	\begin{bmatrix}
		* & * & * & * & * & *\\
		0 & * & * & * & * & *\\
		0 & 0 & * & * & * & *\\
	\end{bmatrix}.
$$
The first step (we will use a Householder transformation for this) is on all of
the $m \times n$ input matrix, with an operation which transforms the first
column.   The second step is on the $(m-1) \times (n-1)$ submatrix obtained by
omitting the top row and left column:
$$
	\left[\begin{array}{r|rrrrr}
		* & * & * & * & * & *\\
		\hline
		0 & * & * & * & * & *\\
		0 & * & * & * & * & *\\
	\end{array}\right].
$$
We can keep operating on submatrices
until there are no more of them.  So, the problem of upper-triangularization
reduces to that of modifying the first column of a matrix to have all non-zero
entries except the entry at the top of the column.

%% ----------------------------------------------------------------
\subsection{Householders for column-zeroing}
\label{sec:hh_col_zero}

The goal is to put a matrix
$$
	A = \begin{bmatrix}
		* & * & * & * & * & *\\
		* & * & * & * & * & *\\
		* & * & * & * & * & *\\
	\end{bmatrix}
$$
into the form
$$
	QA = \begin{bmatrix}
		* & * & * & * & * & *\\
		0 & * & * & * & * & *\\
		0 & * & * & * & * & *\\
	\end{bmatrix}
$$
The key insight is that when $Q$ acts on $A$ by left multiplication, the $j$th
column of $QA$ is the matrix-times-vector product of $Q$ times the $j$th column
of $A$.

Let $\vecu$ = column 0 of $A$, and $\vecv$ = column 0 of $QA$.  At this point
all we know is that we want $\vecv$ to have all zeroes except the first entry.
But if we're going to use Householder transformations, which are
norm-preserving, we know that $\|\vecv\| = \|\vecu\|$.  So, $v_1 = \pm
\|\vecu\|$.  (We'll see in the next paragraph how to choose the postive or
negative sign.) We want a Householder matrix $Q$ which sends $\vecu$ to
$\vecv$, and which can do what it will with the rest of the matrix as a side
effect.

What is the reflection axis $\veca$?   This is just as in section
\ref{sec:Q_spec_uv}:  $\veca = \vecu - \vecv$.  Above we had $v_1 = \pm
\|\vecu\|$.  Since what we're going to do with that is to compute
$\veca=\vecu-\vecv$, we can choose $v_1$ to have the opposite sign of $u_1$, to
avoid the cancellation error which can happen when we subtract two nearly equal
numbers.

Error analysis:  Since these $Q$'s are orthogonal, they're norm-preserving,
and so
$$
	Q(\vecu+\veceps) = Q(\vecu) + Q(\veceps)
$$
but
$$
	\|Q(\veceps)\| \le \|\veceps\|
$$
as discussed in section \ref{sec:err_prop}.

Also, there is no pivoting involved, and thus (other than the choice of the
sign of $v_1$) no \bfidxtwo{data-dependent decisions}{data-dependent decision}.

%% ----------------------------------------------------------------
\subsection{Computation of determinants}
\label{sec:HHdet}

Given a square matrix $A$, we can use repeated Householder transformations to
turn it into an upper-triangular matrix $U$.  As discussed in section
\ref{sec:comp_det}, $\det(A)$ is the product of diagonal entries of $U$,
divided by the product of the determinants of the transformation matrices.
However, as seen in sections \ref{sec:Q_props_alg} and \ref{sec:Q_props_geom},
Householder matrices have determinant $-1$.  So, we just have to remember
whether the number of Householders applied was odd or even.  But it takes $n-1$
Householders to upper-triangularize an $n \times n$ matrix:
$$
	\begin{bmatrix}
	* & * & * & * \\
	0 & * & * & * \\
	0 & 0 & * & * \\
	0 & 0 & 0 & * \\
	\end{bmatrix}
$$
So,
$$
	\det(A) = (-1)^{n-1} \det(U).
$$

%% ----------------------------------------------------------------
\subsection{Computation of matrix inverses}
\label{sec:HHinv}

Given an invertible (square) matrix $A$, we can make an augmented matrix just
as in section \ref{sec:comp_inv}.  We can use Householder
upper-triangularization to put the augmented matrix into upper-triangular form,
i.e. row echelon form.  The rest of the work is in clearing out the
above-diagonal entries of the left-hand side of the augmented matrix, for which
we can use scale and shear matrices as in section \ref{sec:geom_det_and_norm}.

%[Probably should put an example here.]

%% ----------------------------------------------------------------
\subsection{Rotation matrices}
\label{sec:HHrot}

$O(n)$ has two connected components: $SO(n)$ and $O^-(n)$.  The former are
(rigid) \bfidx{rotations} (determinant $+1$); the latter are
\bfidx{reflections} (determinant $-1$).  (See \cite{bib:frankel} or
\cite{bib:simon} for more information on $SO(n)$.) One might ask:  If
Householder transformations provide easily understandable, readily computable
elements of $O^-(n)$, then is there another, equally pleasant technique to
compute elements of $SO(n)$?  In fact there is.  One can consult the literature
for the Rodrigues formula, and modify that.

Alternatively, one can use the fact that since determinants multiply, the
product of two reflections, each having determinant $-1$, will have determinant
1.  Thus the product will be a rotation.  If the same reflection is composed
with itself, the product will be the identity, as seen in sections
\ref{sec:Q_props_alg} and \ref{sec:Q_props_geom}.  But the product of two
\emph{distinct} reflections will be a non-trivial rotation.

Let $\vecu$ and $\vecw$ be given.  Since rotations are norm-preserving, we can
only hope to rotate $\vecu$ into $\vecw$ if they both have the same norm.  We
can help this by computing $\uhat$ and $\what$.  Now note that $\uhat+\what$
will lie halfway between $\uhat$ and $\what$:

\includegraphics{rotuvw.eps}

Let $\vecv=\uhat+\what$, and compute $\vhat$.  Then let $Q_1$ reflect $\uhat$
to $\vhat$.  Just as in section \ref{sec:Q_spec_uv}, this is a reflection about
an axis $\veca_1 = \uhat - \vhat$.  Likewise for the second reflection.  We
then have
$$
	R = Q_2 Q_1 =
	\left(I - 2 \frac{\veca_2 \veca_2^t}{\veca_2^t \veca_2}\right)
	\left(I - 2 \frac{\veca_1 \veca_1^t}{\veca_1^t \veca_1}\right).
$$

%% ----------------------------------------------------------------
\subsection{Software}

The algorithms described in this document have been implemented in the Python
programming language in the files

\begin{itemize}
\item \texttt{http://math.arizona.edu/{\textasciitilde}kerl/python/sackmat\_m.py}
\item \texttt{http://math.arizona.edu/{\textasciitilde}kerl/python/pydet}
\item \texttt{http://math.arizona.edu/{\textasciitilde}kerl/python/pymatinv}
\item \texttt{http://math.arizona.edu/{\textasciitilde}kerl/python/pyhh}
\item \texttt{http://math.arizona.edu/{\textasciitilde}kerl/python/pyrefl}
\item \texttt{http://math.arizona.edu/{\textasciitilde}kerl/python/pyrot}
\end{itemize}

The first is library file containing most of the logic; the remaining four
handle I/O and command-line arguments, and call the library routines.

\begin{itemize}
\item The library routine \texttt{householder\_vector\_to\_Q}
	takes a vector $\veca$ and produces a matrix Q, as described
	in section \ref{sec:geometric_construction}.
\item The program \texttt{pyhh} is a direct interface to this routine.
\item The library routine \texttt{householder\_UT\_pass\_on\_submatrix}
	takes a matrix and a specified row/column number, applying a Householder
	transformation on a submatrix, as described in section \ref{sec:HHUT}.
\item The library routine \texttt{householder\_UT} puts a matrix into
	upper-triangular form, using repeated calls to
	\texttt{householder\_UT\_pass\_on\_submatrix}.
\item The program \texttt{pydet} computes determinants as described in
	section \ref{sec:HHdet}.
\item The program \texttt{pymatinv} computes matrix inverses as described in
	section \ref{sec:HHinv}.
\item The program \texttt{pyrefl} reads two vectors $\vecu$ and $\vecv$,
	normalizes them to $\uhat$ and $\vhat$, and then uses the technique
	described in section \ref{sec:Q_spec_uv} to compute a reflection matrix
	sending $\uhat$ to $\vhat$.
\item The program \texttt{pyrot} reads two vectors $\vecu$ and $\vecv$,
	normalizes them to $\uhat$ and $\vhat$, and then uses the technique
	described in section \ref{sec:HHrot} to compute a rotation matrix sending
	$\uhat$ to $\vhat$.
\end{itemize}

The Python programming language is clear and intuitive.  (Whether my coding
style is equally clear and intuitive is left to the judgment of the reader).
As a result, these algorithms should be readily translatable to, say, Matlab.

%% ----------------------------------------------------------------
\newpage
\section{Acknowledgments}

Thanks to Dr. Ross Martin of Lockheed Martin Corporation and Dr. John Palmer of
the University of Arizona for many productive conversations.
