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:





A
=


(



12


−
51


4




6


167


−
68




−
4


24


−
41



)




{\displaystyle A={\begin{pmatrix}12&-51&4\\6&167&-68\\-4&24&-41\end{pmatrix}}}


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





H
=
I
−


2


u

T


u



u

u

T




{\displaystyle H=I-{\frac {2}{u^{T}u}}uu^{T}}


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





u
=
a
−
∥
a

∥

2




e

1




{\displaystyle u=a-\|a\|_{2}\;e_{1}}


then the transformation reflects 





a




{\displaystyle {\mathit {a}}}

 onto the first standard basis vector






e

1


=
[
1

0

0

.
.
.

]

T




{\displaystyle e_{1}=[1\;0\;0\;...]^{T}}


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}}

:





u
=
a
+


sign


(

a

1


)
∥
a

∥

2




e

1




{\displaystyle u=a+{\textrm {sign}}(a_{1})\|a\|_{2}\;e_{1}}


and normalize with respect to the first element:





v
=


u

u

1






{\displaystyle v={\frac {u}{u_{1}}}}


The equation for 



H


{\displaystyle H}

 thus becomes:





H
=
I
−


2


v

T


v



v

v

T




{\displaystyle H=I-{\frac {2}{v^{T}v}}vv^{T}}


or, in another form





H
=
I
−
β
v

v

T




{\displaystyle H=I-\beta vv^{T}}


with





β
=


2


v

T


v





{\displaystyle \beta ={\frac {2}{v^{T}v}}}


Applying 





H




{\displaystyle {\mathit {H}}}

 on 





a




{\displaystyle {\mathit {a}}}

 then gives





H

a
=
−


sign


(

a

1


)

∥
a

∥

2




e

1




{\displaystyle H\;a=-{\textrm {sign}}(a_{1})\;\|a\|_{2}\;e_{1}}


and applying 





H




{\displaystyle {\mathit {H}}}

 on the matrix 





A




{\displaystyle {\mathit {A}}}

 zeroes all subdiagonal elements of the first column:






H

1



A
=


(




r

11





r

12





r

13






0


∗


∗




0


∗


∗



)




{\displaystyle H_{1}\;A={\begin{pmatrix}r_{11}&r_{12}&r_{13}\\0&*&*\\0&*&*\end{pmatrix}}}


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:






H

2


=


(



1


0


0




0


H





0





)




{\displaystyle H_{2}={\begin{pmatrix}1&0&0\\0&H&\\0&&\end{pmatrix}}}


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}}}

.






H

n



.
.
.


H

3



H

2



H

1


A
=
R


{\displaystyle H_{n}\;...\;H_{3}H_{2}H_{1}A=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}}}

.






H

1



H

2



H

3



.
.
.


H

n


=
Q


{\displaystyle H_{1}H_{2}H_{3}\;...\;H_{n}=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





R

x
=

Q

T



b


{\displaystyle R\;x=Q^{T}\;b}


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.





R
=


(




R

1






0



)




{\displaystyle R={\begin{pmatrix}R_{1}\\0\end{pmatrix}}}


and the same for the RHS:






Q

T



b
=


(




q

1







q

2





)




{\displaystyle Q^{T}\;b={\begin{pmatrix}q_{1}\\q_{2}\end{pmatrix}}}


Finally, solve the square upper triangular system by back substitution:






R

1



x
=

q

1




{\displaystyle R_{1}\;x=q_{1}}


import <Utilities/Math.sl>;
import <Utilities/Sequence.sl>;
import <Utilities/Conversion.sl>;
 
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] );