Lab 6: QR Decomposition using Householder reflections

In this lab, we construct an algorithm for calculating the QR Decomposition of a matrix $A$ using Householder reflections. In this lab we always use the 2-norm.

Householder reflections

Given a vector $\mathbf v$ satisfying $\|\mathbf v\|=1$, the Householder reflection is the orthogonal matrix as

$$Q_{\mathbf v} \triangleq I - 2 \mathbf v \mathbf v^\top$$

Exercise 1(a) Demonstrate numerically that $Q_{\mathbf v}$ is orthogonal on a random vector of length 5 satisfying $\|\mathbf v\|=1$. Prove that this is always the case.

In [ ]:

Exercise 1(b) Given n, how would you construct the length n basis vector $\mathbf e_1$? How about $\mathbf e_k$? (Advanced: do this exercise in one line using [v;w] syntax to concatinate two or more vectors, along with zeros to create vectors of zeros.) Show that your approach works for n=20 and k=5.

In [ ]:

Exercise 1(c) Given a vector $\mathbf x$, define

$$\mathbf w = \mathbf x - \alpha \|{\mathbf x}\| \mathbf e_1, \qquad \alpha = \pm 1,$$

where the choice of $\alpha$ is arbitrary. Demonstate numerically that

$$\|\mathbf w\|^2=2(\|x\|^2 - \alpha \|x \| x_1)$$

where $x_1 = \mathbf e_1^\top \mathbf x$, for a random vector $\mathbf x$ of length 5. (Advanced: Prove this is true.)

In [ ]:

Exercise 1(c) With $\mathbf v = \mathbf w / \|\mathbf w\|$, demonstrate numerically that $Q_{\mathbf v}$ satsifies the property

$$Q_{\mathbf v}\mathbf x = \alpha \|{\mathbf x}\| \mathbf e_1.$$

(Advanced: Prove this is true.)

In [ ]:

Calculating the QR Decomposition

Exercise 2(a) Use 1(c) to design a $\mathbf v$ so that $Q_{\mathbf v} A$ has zeros in the $2:n$ rows of the first column, where $A$ is an $n \times m$ matrix, where $n \geq m$. Demonstrate this work on a random $10 \times 5$ matrix.

In [ ]:

Exercise 2(b) Given $k<n$ and a vector $\mathbf x$ of length $n$, design a $\mathbf v$ so that

$$Q_{\mathbf v} \mathbf x = \begin{pmatrix} x_1 \cr \vdots \cr x_{k-1} \cr \alpha \|x[k:n]\| \cr 0 \cr \vdots \cr 0\end{pmatrix}.$$

Demonstrate that this works on a random vector of length 10 for k=3;

In [ ]:

Exercise 2(c) How would you use 2(b) to construct to choose a $\mathbf v$ so that $Q_{\mathbf v} A$ has zeros in the $(j+1):n$ rows of the $j$th column, and but leaves the $1:(j-1)$ rows alone?

Exercise 2(d) How would you calculate the QR Decomposition using 2(c) with $m-1$ choices of $\mathbf v$ to, column-by-column, introduce zeros below the diagonal?

Implementation

Exercise 3(a) Complete the following functions that returns a Householder reflection corresponding to exercise 2(c). Use α=-sign(x[j]). (This choice will be investigated next lab.)

In [ ]:
function house_col(A,j)
    # TODO: Return a Householder reflection Q so that
    #    QA=Q*A satsifies QA[j+1:n,j]==zeros(n-j) and QA[1:j-1,j]==A[1:j-1,j]
    # where     n=size(A,1)

end

Exercise 3(b) Complete the following function:

In [ ]:
function house_qr(A)
    n,m=size(A)

    Qi=eye(n)
    for j=1:m-1  # introduce zeros in each column
        #TODO: Update Qi and A using house_col(A,j)
    end

    Qi',A
end

Exercise 3(d) Check that Q,R=house_qr(A) returns the correct result on a random $10 \times 5$ matrix.

In [ ]:
A=rand(10,5)
Q,R=house_qr(A)
norm(Q*R-A)

Complexity

Recall that the complexity of an algorithm gives the order of the number of operations. Remember, we saw in lectures that the number of FLOPs for backsubstution was $n^2 + 2n$, which was an order of complexity of $O(n^2)$.

We can estimate order of complexity by timing the algorithm. To do this, we use two commands: @time, which prints the time it takes to run a calculation and @elapsed, which returns the time it takes to run a computation.

Exercise 5(a) Explain the following block of code:

In [150]:
N=300
tms=zeros(N)

for k=1:N
    A=rand(k,k)
    tms[k]=@elapsed (qr(A))
end

plot(1:N,tms);

A loglog plot scales both the x and y axis logarithmically, and is a useful tool for studying complexity, as the slope of $cn^\beta$ is dictated by $\beta$:

In [148]:
r=1:10000

loglog(r,r)
loglog(r,2r)
loglog(r,r.^2)
loglog(r,2.5r.^2)
loglog(r,r.^3)
loglog(r,3r.^3);

It is also convenient to use a logarithmic choice of sample points, given by logspace(a,b,n), which returns n points between $10^a$ and $10^b$:

In [149]:
r=logspace(0,4,10)  # 10 points sampled logarithmically between 1 and 10^4

loglog(r,r)
loglog(r,2r)
loglog(r,r.^2)
loglog(r,2.5r.^2)
loglog(r,r.^3)
loglog(r,3r.^3);

Exercise 5(b) Construct a vector logtms of the time taken for QR on a $k \times k$ matrix, for k = floor(Int,j) for each j in logspace(0,3.2,10). Plot a loglog of the results.

In [ ]:

Exercise 5(c) Can you guess the complexity of qr? That is, it has complexity of $O(n^\beta)$ for what choice of $\beta$? Verify your guess by comparing a loglog plot of $0.00000001n^\beta$ to the plot of logtms.

In [ ]:

Exercise 5(d) Repeat 5(c) for house_qr, for logspace(1,2.5,10) (for speed reasons).

In [ ]: