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 std.stdio, std.math, std.algorithm, std.traits,
       std.typecons, std.numeric, std.range, std.conv;
 
template elementwiseMat(string op) {
    T[][] elementwiseMat(T)(in T[][] A, in T B) pure nothrow {
        if (A.empty)
            return null;
        auto R = new typeof(return)(A.length, A[0].length);
        foreach (immutable r, const row; A)
            R[r][] = mixin("row[] " ~ op ~ "B");
        return R;
    }
 
    T[][] elementwiseMat(T, U)(in T[][] A, in U[][] B)
    pure nothrow if (is(Unqual!T == Unqual!U)) {
        assert(A.length == B.length);
        if (A.empty)
            return null;
        auto R = new typeof(return)(A.length, A[0].length);
        foreach (immutable r, const row; A) {
            assert(row.length == B[r].length);
            R[r][] = mixin("row[] " ~ op ~ "B[r][]");
        }
        return R;
    }
}
 
alias mSum = elementwiseMat!q{ + },
      mSub = elementwiseMat!q{ - },
      pMul = elementwiseMat!q{ * },
      pDiv = elementwiseMat!q{ / };
 
bool isRectangular(T)(in T[][] mat) pure nothrow {
    return mat.all!(r => r.length == mat[0].length);
}
 
T[][] matMul(T)(in T[][] a, in T[][] b) pure nothrow
in {
    assert(a.isRectangular && b.isRectangular &&
           a[0].length == b.length);
} body {
    auto result = new T[][](a.length, b[0].length);
    auto aux = new T[b.length];
    foreach (immutable j; 0 .. b[0].length) {
        foreach (immutable k; 0 .. b.length)
            aux[k] = b[k][j];
        foreach (immutable i; 0 .. a.length)
            result[i][j] = a[i].dotProduct(aux);
    }
    return result;
}
 
Unqual!T[][] transpose(T)(in T[][] m) pure nothrow {
    auto r = new Unqual!T[][](m[0].length, m.length);
    foreach (immutable nr, row; m)
        foreach (immutable nc, immutable c; row)
            r[nc][nr] = c;
    return r;
}
 
T norm(T)(in T[][] m) pure nothrow {
    return transversal(m, 0).map!q{ a ^^ 2 }.sum.sqrt;
}
 
Unqual!T[][] makeUnitVector(T)(in size_t dim) pure nothrow {
    auto result = new Unqual!T[][](dim, 1);
    foreach (row; result)
        row[] = 0;
    result[0][0] = 1;
    return result;
}
 
/// Return a nxn identity matrix.
Unqual!T[][] matId(T)(in size_t n) pure nothrow {
    auto Id = new Unqual!T[][](n, n);
    foreach (immutable r, row; Id) {
        row[] = 0;
        row[r] = 1;
    }
    return Id;
}
 
T[][] slice2D(T)(in T[][] A,
                 in size_t ma, in size_t mb,
                 in size_t na, in size_t nb) pure nothrow {
    auto B = new T[][](mb - ma + 1, nb - na + 1);
    foreach (immutable i, brow; B)
        brow[] = A[ma + i][na .. na + brow.length];
    return B;
}
 
size_t rows(T)(in T[][] A) pure nothrow { return A.length; }
 
size_t cols(T)(in T[][] A) pure nothrow {
    return A.length ? A[0].length : 0;
}
 
T[][] mcol(T)(in T[][] A, in size_t n) pure nothrow {
    return slice2D(A, 0, A.rows - 1, n, n);
}
 
T[][] matEmbed(T)(in T[][] A, in T[][] B,
                  in size_t row, in size_t col) pure nothrow {
    auto C = new T[][](rows(A), cols(A));
    foreach (immutable i, const arow; A)
        C[i][] = arow[]; // Some wasted copies.
    foreach (immutable i, const brow; B)
        C[row + i][col .. col + brow.length] = brow[];
    return C;
}
 
// Main routines ---------------
 
T[][] makeHouseholder(T)(in T[][] a) {
    immutable m = a.rows;
    immutable T s = a[0][0].sgn;
    immutable e = makeUnitVector!T(m);
    immutable u = mSum(a, pMul(e, a.norm * s));
    immutable v = pDiv(u, u[0][0]);
    immutable beta = 2.0 / v.transpose.matMul(v)[0][0];
    return mSub(matId!T(m), pMul(v.matMul(v.transpose), beta));
}
 
Tuple!(T[][],"Q", T[][],"R") QRdecomposition(T)(T[][] A) {
    immutable m = A.rows;
    immutable n = A.cols;
    auto Q = matId!T(m);
 
    // Work on n columns of A.
    foreach (immutable i; 0 .. (m == n ? n - 1 : n)) {
        // Select the i-th submatrix. For i=0 this means the original
        // matrix A.
        immutable B = slice2D(A, i, m - 1, i, n - 1);
 
        // Take the first column of the current submatrix B.
        immutable x = mcol(B, 0);
 
        // Create the Householder matrix for the column and embed it
        // into an mxm identity.
        immutable H = matEmbed(matId!T(m), x.makeHouseholder, i, i);
 
        // The product of all H matrices from the right hand side is
        // the orthogonal matrix Q.
        Q = Q.matMul(H);
 
        // The product of all H matrices with A from the LHS is the
        // upper triangular matrix R.
        A  = H.matMul(A);
    }
 
    // Return Q and R.
    return typeof(return)(Q, A);
}
 
// Polynomial regression ---------------
 
/// Solve an upper triangular system by back substitution.
T[][] solveUpperTriangular(T)(in T[][] R, in T[][] b) pure nothrow {
    immutable n = R.cols;
    auto x = new T[][](n, 1);
 
    foreach_reverse (immutable k; 0 .. n) {
        T tot = 0;
        foreach (immutable j; k + 1 .. n)
            tot += R[k][j] * x[j][0];
        x[k][0] = (b[k][0] - tot) / R[k][k];
    }
 
    return x;
}
 
/// Solve a linear least squares problem by QR decomposition.
T[][] lsqr(T)(T[][] A, in T[][] b) pure nothrow {
    const qr = A.QRdecomposition;
    immutable n = qr.R.cols;
    return solveUpperTriangular(
        slice2D(qr.R, 0, n - 1, 0, n - 1),
        slice2D(qr.Q.transpose.matMul(b), 0, n - 1, 0, 0));
}
 
T[][] polyFit(T)(in T[][] x, in T[][] y, in size_t n) pure nothrow {
    immutable size_t m = x.cols;
    auto A = new T[][](m, n + 1);
    foreach (immutable i, row; A)
        foreach (immutable j, ref item; row)
            item = x[0][i] ^^ j;
    return lsqr(A, y.transpose);
}
 
void main() {
    // immutable (Q, R) = QRdecomposition([[12.0, -51,   4],
    immutable qr = QRdecomposition([[12.0, -51,   4],
                                    [ 6.0, 167, -68],
                                    [-4.0,  24, -41]]);
    immutable form = "[%([%(%2.3f, %)]%|,\n %)]\n";
    writefln(form, qr.Q);
    writefln(form, qr.R);
 
    immutable x = [[0.0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]];
    immutable y = [[1.0, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321]];
    polyFit(x, y, 2).writeln;
}