You are given a full column rank matrix $A_{m\times n}.$ You
want to premultiply it with an orthogonal $P$ to get $PA =
R,$ an upper triangular matrix. Then with $Q = P'$ you
get a $QR$ decomposition.
We plan to construct $P$ as a product $P_nP_{n-1}\cdots
P_1,$ where each $P_i$ is an orthogonal matrix such that
if $A$ is like this:
then $P_1A$ is like this (the shaded entries are 0's):
Then $P_2P_1A$ is like
and $P_3P_2P_1A$ is like
The question is how to get such $P_i$'s. Let's focus
on $P_1.$ If $\v a$ is the first column of $A,$
the $P\v a$ is the first column of $PA.$ So we want
$P\v a$ to have all entries (except the top one) to be zero.
If, for example, $\v a = \left[\begin{array}{ccccccccccc}1\\2\\3
\end{array}\right]$ we must
have $P\v a = \left[\begin{array}{ccccccccccc}x\\0\\0
\end{array}\right].$ Since $\|\v a\| = \|P
\v a\|,$ hence we must have $x = \pm \|\v a\| = \pm \sqrt{14}.$
Now $\left[\begin{array}{ccccccccccc}1\\2\\3
\end{array}\right]$
and $\left[\begin{array}{ccccccccccc}\sqrt{14}\\0\\0
\end{array}\right]$ are two vectors of equal
length in ${\mathbb R}^3,$ and hence we can rotate one to the
other. Householder's transform is just such an idea except that
it works with reflection rather than rotation.
Given two distinct vectors $\v a$ and $\v b$ of equal
length, Householder's transform is an orthogonal matrix $H$ such that
$H\v a = \v b$ and $H\v b = \v a.$
Intuitively it is a reflection along the hyperplane bisecting the
angle between $\v a$ and $\v b.$
Algebraically, define $\v u$ to be the unit vector along $\v a - \v b.$
Then find an $\alpha\in{\mathbb R}$ such that $H = I+ \alpha \v
u\v u'$ is an orthogonal matrix (just find $H'H$, and see
what value of $\alpha$ makes it equal to $I.$).
It is trivial to check that this $H$ works. This $H$ is
called a Householder matrix.
Notice that though $H$ is $n\times n$ it is completely
determined by $\v u$ (ie, by just $n$ numbers). We
shall call it $H(\v u).$
Also computing $H\v x$ is easy: $H(\v u)\v x = (I+ \alpha
\v u \v u')\v x = \v x +(\alpha \v u'\v x)\v u,$ which
requires $2n+1$ multiplications instead of $n^2$
multiplications (needed for multiplying an $n\times n$ matrix
and an $n\times1$ vector).
It is best explained via an example.
Let
$$A = \left[\begin{array}{ccccccccccc}1 & 2 & 3\\4 & 5 & 6\\7 & 8 & 1 \\ 2 & 3 &4
\end{array}\right]$$
After first step we want to make it
$$A = \left[\begin{array}{ccccccccccc}\sqrt{70} & * & *\\0 & * & *\\0 & * & * \\ 0 & *
&*
\end{array}\right].$$
For this we premultiply by $H(\v u_1)$ to
convert $\v a = \left[\begin{array}{ccccccccccc}1\\4\\7\\2
\end{array}\right]$ to $\v b
= \left[\begin{array}{ccccccccccc}\sqrt{70}\\0\\0\\0
\end{array}\right].$
Clearly $\v u_1 = $ unit vector along $\v a-\v b$ will do
the job.
After you compute the $*$ entries, repeat the same process
with
$$A = \left[\begin{array}{ccccccccccc}\sqrt{70} & * & *\\0 & a & *\\0 & b & * \\ 0 & c &*
\end{array}\right]$$
to get
$$A = \left[\begin{array}{ccccccccccc}\sqrt{70} & * & *\\0 & \sqrt{a^2+b^2+c^2} & *\\0 & 0
& * \\ 0 & 0 &*
\end{array}\right].$$
Here we need to multiply with
$$
\left[\begin{array}{ccccccccccc}
1 & \v 0'\\
\v 0 & H(\v u_2)
\end{array}\right],
$$
where $\v u_2$ is the $3\times 1$ unit vector along
$$
\left[\begin{array}{ccccccccccc}\sqrt{a^2+b^2+c^2}\\0\\0
\end{array}\right] - \left[\begin{array}{ccccccccccc}a\\b\\c
\end{array}\right].
$$
Similarly after the third (and final) step we obtain
$$\left[\begin{array}{ccccccccccc}* & * & *\\0 & * & *\\0 & 0 & * \\ 0 & 0 &0
\end{array}\right],$$
which is our required $R.$ In this algorithm e never
explicitly compute and store $H.$ Instead, we just store the
sequence $\v u_1,\v u_2$ and $\v u_3.$
You are to implement the above algorithm in C (preferred) or
R. If you are using C, then your program must have the following
interface:
First it should read the numbers of rows and columns of the
matrix.
Then it should read the entries of the matrix row by row.
The program should print the final $R$ and the sequence
of $\v u_i$'s.
If you are writing an R program, then you submit a function that
takes the matrix as ints only input. The function should return
(not print) the R matrix as well as the $\v u_i$ vectors.
1 marks bonus (above the regular 10 marks) for a program that
works even if $A$ is not full column rank.
