Assignment 2

This lab explores two alternative decompositions to the LU and QR decompositions: the UL Decomposition and LQ decompositions.

UL Decomposition

An LU decomposition is constructed by multiplying by a lower triangular matrix that introduces zeros column-by-column starting with the first column.

Instead, we will investigate the UL decomposition: we will introduce zeros column-by-column starting with the last column.

Exercise 1(a) Given an $n \times n$ matrix A, what upper triangular matrix U_j introduces zeros in the $j$th column:

(U_j*A)[1:j-1,j] == zeros(1:j-1)?

Complete the following function that takes in a square Matrix A and an integer j and returns an upper triangular Matrix U_j that satisfies the above property.

In [11]:
# INPUT:  A square matrix A and integer j
# OUTPUT: A matrix U_j that is upper triangular

function upper_reduce_col(A,j)
    n=size(A,1)
    U=eye(n)

    for k=1:j-1
        U[k,j]=-A[k,j]/A[j,j]
    end

    U
end

A=rand(5,5)
U=upper_reduce_col(A,4)
U*A
Out[11]:
5x5 Array{Float64,2}:
 0.079882  -0.0327976  -0.312991   1.11022e-16  -0.0893706
 0.817018   0.235419    0.0803426  0.0           0.508066 
 0.165188   0.340396    0.713421   0.0           0.326151 
 0.456192   0.702257    0.656478   0.892522      0.481106 
 0.843693   0.24602     0.341197   0.703776      0.0503664

Exercise 1(c) Check that

n=100
A=rand(n,n)
U_n=upper_reduce_col(A,n)
norm((U_n*A)[1:n-1,n])

is small.

In [12]:
n=100
A=rand(n,n)
U_n=upper_reduce_col(A,n)
norm((U_n*A)[1:n-1,n])
Out[12]:
3.764949453935611e-16

Exercise 1(d) Complete the following function that returns a UL Decomposition.

In [14]:
# INPUT:  A square matrix A
# OUTPUT: A tuple (U,L) where U is upper triangular and L is lower triangular

function ul(A)
    n=size(A,1)

    Ui=eye(n)
    for j=n:-1:2
        Uj=upper_reduce_col(A,j)
        Ui=Uj*Ui
        A=Uj*A
    end

    inv(Ui),A
end
Out[14]:
ul (generic function with 1 method)

Exercise 1(e) Check that the following is small:

In [16]:
n=100
A=rand(n,n)
U,L=ul(A)

norm(U*L-A)
Out[16]:
4.935716917863025e-10

Badly conditioned PLU matrix

Recall from lectures the matrix that breaks the PLU Decomposition

$$A_n \triangleq \begin{pmatrix} 1 & && & 1 \cr -1 & 1 && & 1\cr -1 & -1 & 1 && 1 \cr \vdots & \vdots & \vdots & \ddots & \vdots \cr -1 & -1 & -1 & \cdots & 1 \end{pmatrix}$$

Exercise 2(a) Write a routine bad_plu_mat(n) that takes in an integer n and returns an n x n Matrix $A_n$ as defined above. Use a different approach to construct the matrix than what was done in lectures.

In [20]:
# INPUT:  An Integer n
# OUTPUT: A Matrix A_n that is n x n

function bad_plu_mat(n)
    A=2eye(n)-tril(ones(n,n)) 
    A[1:n-1,end]=ones(n-1)
    A
end

bad_plu_mat(3),bad_plu_mat(4)
Out[20]:
(
3x3 Array{Float64,2}:
  1.0   0.0  1.0
 -1.0   1.0  1.0
 -1.0  -1.0  1.0,

4x4 Array{Float64,2}:
  1.0   0.0   0.0  1.0
 -1.0   1.0   0.0  1.0
 -1.0  -1.0   1.0  1.0
 -1.0  -1.0  -1.0  1.0)

Exercise 2(b) Plot the condition numbers ${\rm cond}(L)*{\rm cond}(U)$ for $A_n = P L U$ as a function of n, for n=1:100, using semilogy.

In [21]:
using PyPlot

cnd=Float64[(A=bad_plu_mat(n);
    (L,U,p)=lu(A); 
    cond(L)*cond(U)) for n=1:100]

semilogy(cnd)
Out[21]:
1-element Array{Any,1}:
 PyObject <matplotlib.lines.Line2D object at 0x31329bf10>

Exercise 2(c) Plot the condition numbers ${\rm cond}(U)*{\rm cond}(L)$ for $A_n = UL$ as a function of n, for n=1:100, using semilogy. What do you conclude about the stability of the UL decomposition for this example?

In [22]:
cnd=Float64[(A=bad_plu_mat(n);
    (U,L)=ul(A); 
    cond(U)*cond(L)) for n=1:100]

semilogy(cnd);

LQ Decomposition

This part investigates the calculation of $A=LQ$, which is an alternative to $A=QR$.

Exercise 2(a) Write a function house(x,j) that takes in a Vector x and returns an orthogonal Matrix $Q$ so that

$$x^\top Q = [x_1,\ldots,x_{j-1},-{\rm sign}\, x_j \|x_{j:end}\|,0,\ldots,0]$$
In [29]:
# INPUT:  A Vector x and Integer j
# OUTPUT: A Matrix Q that is m x m and orthogonal

function house(x,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]
    # for     n=size(A,1)
    n=length(x)
    x=[zeros(j-1);x[j:end]]
    α=-sign(x[j])
    v=x-α*norm(x)*[zeros(j-1);1;zeros(n-j)]
    v=v/norm(v)
    I-2v*v'    
end

x=rand(5)
Q=house(x,2)
x'*Q
Out[29]:
1x5 Array{Float64,2}:
 0.906057  -1.39361  -5.55112e-17  -1.73472e-16  -2.22045e-16

Exercise 2(b) Check numerically that

$$\|x^\top Q - [x_1,\ldots,x_{j-1},-{\rm sign}\, x_j \|x_{j:end}\|,0,\ldots,0]\|$$

is small for a random vector of length 100 and j=3.

In [33]:
x=rand(100)
Q=house(x,3)
norm(Q*x-[x[1];x[2];-sign(x[3])*norm(x[3:end]);zeros(length(x)-3)])
Out[33]:
2.479007607181183e-15

Exercise 2(c) Use part 2(a) to write a function lq(A) that calculates the RQ Decomposition: it takes in a n x m Matrix A and returns a tuple of matrices L and Q so that L is n x m and lower triangular, Q is m x m and orthogonal, and L*Q is (up to rounding error) equal to A.

In [38]:
# INPUT:  An n x m matrix A
# OUTPUT: A tuple (L,Q)  where L is n x m lower triangular and Q is m x m and orthogonal

function lq(A)
    n,m=size(A)
    Qt=eye(m)
    for k=1:min(n,m)
        Q=house(A[k,:],k)
        A=A*Q
        Qt=Qt*Q
    end

    A,Qt'
end

A=rand(5,4)
L,Q=lq(A)
norm(L*Q-A),norm(Q'*Q-I)

A=rand(4,5)
L,Q=lq(A)
norm(L*Q-A),norm(Q'*Q-I)

A=rand(4,4)
L,Q=lq(A)
norm(L*Q-A),norm(Q'*Q-I)
Out[38]:
(2.5104854386371074e-15,1.677820222479436e-15)

Exercise 2(d) Demonstrate numerically that your lq routine satisfies $\|A-LQ\|$ is small for random matrices for the cases $n > m$, $n=m$ and $n < m$.

In [39]:
A=rand(10,5)
L,Q=lq(A)
norm(A-L*Q)
Out[39]:
1.2216993950229764e-15
In [40]:
A=rand(5,10)
L,Q=lq(A)
norm(A-L*Q)
Out[40]:
8.420114109796008e-16
In [41]:
A=rand(10,10)
L,Q=lq(A)
norm(A-L*Q)
Out[41]:
2.7239005432278633e-15
In [ ]: