Problem Statement

Any rectangular

m × n

{\displaystyle m\times n}

matrix

A

{\displaystyle {\mathit {A}}}

can be decomposed to a product of an orthogonal matrix

Q

{\displaystyle {\mathit {Q}}}

and an upper (right) triangular matrix

R

{\displaystyle {\mathit {R}}}

, as described in QR decomposition. Task Demonstrate the QR decomposition on the example matrix from the Wikipedia article: and the usage for linear least squares problems on the example from Polynomial regression. The method of Householder reflections should be used: Method Multiplying a given vector

a

{\displaystyle {\mathit {a}}}

, for example the first column of matrix

A

{\displaystyle {\mathit {A}}}

, with the Householder matrix

H

{\displaystyle {\mathit {H}}}

, which is given as reflects

a

{\displaystyle {\mathit {a}}}

about a plane given by its normal vector

u

{\displaystyle {\mathit {u}}}

. When the normal vector of the plane

u

{\displaystyle {\mathit {u}}}

is given as then the transformation reflects

a

{\displaystyle {\mathit {a}}}

onto the first standard basis vector which means that all entries but the first become zero. To avoid numerical cancellation errors, we should take the opposite sign of

a

1

{\displaystyle a_{1}}

: and normalize with respect to the first element: The equation for

H

{\displaystyle H}

thus becomes: or, in another form with Applying

H

{\displaystyle {\mathit {H}}}

on

a

{\displaystyle {\mathit {a}}}

then gives and applying

H

{\displaystyle {\mathit {H}}}

on the matrix

A

{\displaystyle {\mathit {A}}}

zeroes all subdiagonal elements of the first column: In the second step, the second column of

A

{\displaystyle {\mathit {A}}}

, we want to zero all elements but the first two, which means that we have to calculate

H

{\displaystyle {\mathit {H}}}

with the first column of the submatrix (denoted *), not on the whole second column of

A

{\displaystyle {\mathit {A}}}

. To get

H

2

{\displaystyle H_{2}}

, we then embed the new

H

{\displaystyle {\mathit {H}}}

into an

m × n

{\displaystyle m\times n}

identity: This is how we can, column by column, remove all subdiagonal elements of

A

{\displaystyle {\mathit {A}}}

and thus transform it into

R

{\displaystyle {\mathit {R}}}

. The product of all the Householder matrices

H

{\displaystyle {\mathit {H}}}

, for every column, in reverse order, will then yield the orthogonal matrix

Q

{\displaystyle {\mathit {Q}}}

. The QR decomposition should then be used to solve linear least squares (Multiple regression) problems

A

x

b

{\displaystyle {\mathit {A}}x=b}

by solving When

R

{\displaystyle {\mathit {R}}}

is not square, i.e.

m

n

{\displaystyle m>n}

we have to cut off the

m

− n

{\displaystyle {\mathit {m}}-n}

zero padded bottom rows. and the same for the RHS: Finally, solve the square upper triangular system by back substitution:

Let's start with the solution:

import ;
import ;
import ;

main :=
    let
        qrTest := [[12.0, -51.0,   4.0],
                   [ 6.0, 167.0, -68.0],
                   [-4.0,  24.0, -41.0]];
        
        qrResult := qr(qrTest);
        
        x := 1.0*(0 ... 10);
        y := 1.0*[1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321];
        
        regResult := polyfit(x, y, 2);
    in
        "q:\n" ++ delimit(delimit(floatToString(qrResult[1], 6), ','), '\n') ++ "\n\n" ++ 
        "r:\n" ++ delimit(delimit(floatToString(qrResult[2], 1), ','), '\n') ++ "\n\n" ++
        "polyfit:\n" ++ "[" ++ delimit(floatToString(regResult, 1), ',') ++ "]";

//---Polynomial Regression---

polyfit(x(1), y(1), n) :=
    let
        a[j] := x ^ j foreach j within 0 ... n;
    in  
        lsqr(transpose(a), transpose([y]));
    
lsqr(a(2), b(2)) :=
    let
        qrDecomp := qr(a);
        prod := mm(transpose(qrDecomp[1]), b);
    in
        solveUT(qrDecomp[2], prod);
        
solveUT(r(2), b(2)) := 
    let 
        n := size(r[1]);
    in
        solveUTHelper(r, b, n, duplicate(0.0, n)); 

solveUTHelper(r(2), b(2), k, x(1)) :=
    let
        n := size(r[1]);
        newX :=  setElementAt(x, k, (b[k][1] - sum(r[k][(k+1) ... n] * x[(k+1) ... n])) / r[k][k]);
    in
        x when k <= 0
    else
        solveUTHelper(r, b, k - 1, newX);

//---QR Decomposition---

qr(A(2)) := qrHelper(A, id(size(A)), 1);

qrHelper(A(2), Q(2), i) :=
    let
        m := size(A);
        n := size(A[1]);
        
        householder := makeHouseholder(A[i ... m, i]);
        
        H[j,k] := 
                householder[j - i + 1][k - i + 1] when j >= i and k >= i 
            else
                1.0 when j = k else 0.0
            foreach j within 1 ... m,
                    k within 1 ... m;
    in
        [Q,A] when i > (n - 1 when m = n else n)
    else
        qrHelper(mm(H, A), mm(Q, H), i + 1);
    

makeHouseholder(a(1)) := 
    let
        v := [1.0] ++ tail(a / (a[1] + sqrt(sum(a ^ 2)) * sign(a[1])));
        
        H := id(size(a)) - (2.0 / mm([v], transpose([v])))[1,1] * mm(transpose([v]), [v]);
    in
        H;

//---Utilities---

id(n)[i,j] := 1.0 when i = j else 0.0
              foreach i within 1 ... n,
                      j within 1 ... n;
                        
mm(A(2), B(2))[i,j] := sum( A[i] * transpose(B)[j] );