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


-- demo/rosettacode/QRdecomposition.exw
with javascript_semantics

function matrix_mul(sequence a, b)
    integer arows = ~a, acols = ~a[1],
            brows = ~b, bcols = ~b[1]
    if acols!=brows then return 0 end if
    sequence c = repeat(repeat(0,bcols),arows)
    for i=1 to arows do
        for j=1 to bcols do
            for k=1 to acols do
                c[i][j] += a[i][k]*b[k][j]
            end for
        end for
    end for
    return c
end function

function vtranspose(sequence v)
-- transpose a vector of length m into an mx1 matrix, 
--                       eg {1,2,3} -> {{1},{2},{3}}
    integer l = length(v)
    sequence res = repeat(0,l)
    for i=1 to l do res[i] = {v[i]} end for
    return res
end function

function mat_col(sequence a, integer col)
    integer la = length(a)
    sequence res = repeat(0,la)
    for i=col to la do
        res[i] = a[i,col]
    end for
    return res
end function

function mat_norm(sequence a)
    atom res = 0
    for i=1 to length(a) do
        res += a[i]*a[i]
    end for
    res = sqrt(res)
    return res
end function

function mat_ident(integer n)
    sequence res = repeat(repeat(0,n),n)
    for i=1 to n do
        res[i,i] = 1
    end for
    return res
end function

function QRHouseholder(sequence a)
    integer cols = length(a[1]),
            rows = length(a),
            m = max(cols,rows),
            n = min(rows,cols)
    sequence q, I = mat_ident(m), Q = I, u, v
--
-- Programming note: The code of this main loop was not as easily
-- written as the first glance might suggest. Explicitly setting 
-- to 0 any a[i,j] [etc] that should be 0 but have inadvertently 
-- gotten set to +/-1e-15 or thereabouts may be advisable. The
-- commented-out code was retrieved from a backup and should be
-- treated as an example and not be trusted (iirc, it made no
-- difference to the test cases used, so I deleted it, and then
-- had second thoughts about it a few days later).
--
    for j=1 to min(m-1,n) do
        u = mat_col(a,j)
        u[j] -= mat_norm(u)
        v = sq_div(u,mat_norm(u))
        q = sq_sub(I,sq_mul(2,matrix_mul(vtranspose(v),{v})))
        a = matrix_mul(q,a)
--      for row=j+1 to length(a) do
--          a[row][j] = 0
--      end for
        Q = matrix_mul(Q,q)
    end for

    -- Get the upper triangular matrix R.
    sequence R = repeat(repeat(0,n),m)
    for i=1 to n do -- (logically 1 to m(>=n), but no need)
        for j=i to n do
            R[i,j] = a[i,j]
        end for
    end for
        
    return {Q,R}
end function

constant a = {{12, -51,   4},
              { 6, 167, -68},
              {-4,  24, -41}}

sequence {q,r} = QRHouseholder(a)

ppOpt({pp_Nest,1,pp_IntFmt,"%4d",pp_FltFmt,"%4g",pp_IntCh,false})
?"A"        pp(a)
?"Q"        pp(q)
?"R"        pp(r)
?"Q * R"    pp(matrix_mul(q,r))

function matrix_transpose(sequence mat)
    integer rows = length(mat),
            cols = length(mat[1])
    sequence res = repeat(repeat(0,rows),cols)
    for r=1 to rows do
        for c=1 to cols do
            res[c][r] = mat[r][c]
        end for
    end for
    return res
end function

--?"Q * Q'" pp(matrix_mul(q,matrix_transpose(q))) -- (~1e-16s)
?"Q * Q`"   pp(sq_round(matrix_mul(q,matrix_transpose(q)),1e15))

procedure least_squares()
    sequence x = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10},
             y = {1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321},
             a = repeat(repeat(0,3),length(x))
    for i=1 to length(x) do
        for j=1 to 3 do
            a[i,j] = power(x[i],j-1)
        end for
    end for
    {q,r} = QRHouseholder(a)
    sequence t = matrix_transpose(q),
             b = matrix_mul(t,vtranspose(y)),
             z = repeat(0,3)
    for k=3 to 1 by -1 do
        atom s = 0
        if k<3 then
            for j = k+1 to 3 do
                s += r[k,j]*z[j]
            end for
        end if
        z[k] = (b[k][1]-s)/r[k,k]
    end for
    printf(1,"Least-squares solution:\n")
--  printf(1," %v\n",{z})                   -- {1.0,2.0.3,0}
--  printf(1," %v\n",{sq_sub(z,{1,2,3})})   -- (+/- ~1e-14s)
    printf(1," %v\n",{sq_round(z,1e13)})    -- {1,2,3}
end procedure
least_squares()
