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


signature RADCATFIELD = sig
type real
val zero : real
val one : real
val + : real * real -> real
val - : real * real -> real
val * : real * real -> real
val / : real * real -> real
val sign : real -> real
val sqrt : real -> real
end
 
functor QR(F: RADCATFIELD) = struct
structure A = struct
local
    open Array
in
fun unitVector n = tabulate (n, fn i => if i=0 then F.one else F.zero)
fun map f x = tabulate(length x, fn i => f(sub(x,i)))
fun map2 f (x, y) = tabulate(length x, fn i => f(sub(x,i),sub(y,i)))
val op + = map2 F.+
val op - = map2 F.-
val op * = map2 F.*
fun multc(c,x) = array(length x,c)*x
fun dot (x,y) = foldl F.+ F.zero (x*y)
fun outer f (x,y) =
    Array2.tabulate Array2.RowMajor (length x, length y,
				     fn (i,j) => f(sub(x,i),sub(y,j)))
fun copy x = map (fn x => x) x
fun fromVector v = tabulate(Vector.length v, fn i => Vector.sub(v,i))
fun slice(x,i,sz) =
    let	open ArraySlice
	val s = slice(x,i,sz)
    in Array.tabulate(length s, fn i => sub(s,i)) end
end
end
structure M = struct
local
    open Array2
in
fun map f x = tabulate RowMajor (nRows x, nCols x, fn (i,j) => f(sub(x,i,j)))
fun map2 f (x, y) =
    tabulate RowMajor (nRows x, nCols x, fn (i,j) => f(sub(x,i,j),sub(y,i,j)))
fun scalarMatrix(m, x) = tabulate RowMajor (m,m,fn (i,j) => if i=j then x else F.zero)
fun multc(c, x) = map (fn xij => F.*(c,xij)) x
val op + = map2 F.+
val op - = map2 F.-
fun column(x,i) = A.fromVector(Array2.column(x,i))
fun row(x,i) = A.fromVector(Array2.row(x,i))
fun x*y = tabulate RowMajor (nRows x, nCols y,
			     fn (i,j) => A.dot(row(x,i), column(y,j)))
fun multa(x,a) = Array.tabulate (nRows x, fn i => A.dot(row(x,i), a))
fun copy x = map (fn x => x) x
fun subMatrix(h, i1, i2, j1, j2) =
    tabulate RowMajor (Int.+(Int.-(i2,i1),1),
		       Int.+(Int.-(j2,j1),1),
		       fn (a,b) => sub(h,Int.+(i1,a),Int.+(j1,b)))
fun transpose m = tabulate RowMajor (nCols m,
				     nRows m,
				     fn (i,j) => sub(m,j,i))
fun updateSubMatrix(h,i,j,s) =
    tabulate RowMajor (nRows s, nCols s, fn (a,b) => update(h,Int.+(i,a),Int.+(j,b),sub(s,a,b)))
end
end
fun toList a =
    List.tabulate(Array2.nRows a, fn i => List.tabulate(Array2.nCols a, fn j => Array2.sub(a,i,j)))
fun householder a =
    let open Array
	val m = length a
	val len = F.sqrt(A.dot(a,a))
	val u = A.+(a, A.multc(F.*(len,F.sign(sub(a,0))), A.unitVector m))
	val v = A.multc(F./(F.one,sub(u,0)), u)
	val beta = F./(F.+(F.one,F.one),A.dot(v,v))
    in
	M.-(M.scalarMatrix(m,F.one), M.multc(beta,A.outer F.* (v,v)))
    end
fun qr mat =
    let open Array2
	val (m,n) = dimensions mat
	val upperIndex = if m=n then Int.-(n,1) else n
	fun loop(i,qm,rm) = if i=upperIndex then {q=qm,r=rm} else 
			    let val x = A.slice(A.fromVector(column(rm,i)),i,NONE)
				val h = M.scalarMatrix(m,F.one)
				val _ = M.updateSubMatrix(h,i,i,householder x)
			    in
				loop(Int.+(i,1), M.*(qm,h), M.*(h,rm))
			    end
    in
	loop(0, M.scalarMatrix(m,F.one), mat)
    end
fun solveUpperTriangular(r,b) =
    let open Array
	val n = Array2.nCols r
	val x = array(n, F.zero)
	fun loop k =
	    let val index = Int.min(Int.-(n,1),Int.+(k,1))
		val _ = update(x,k,
			       F./(F.-(sub(b,k),
				       A.dot(A.slice(x,index,NONE),
					     A.slice(M.row(r,k),index,NONE))),
				   Array2.sub(r,k,k)))
	    in
		if k=0 then x else loop(Int.-(k,1))
	    end
    in
	loop (Int.-(n,1))
    end
fun lsqr(a,b) =
    let val {q,r} = qr a
	val n = Array2.nCols r
    in
	solveUpperTriangular(M.subMatrix(r, 0, Int.-(n,1), 0, Int.-(n,1)),
			     M.multa(M.transpose(q), b))
    end
fun pow(x,1) = x
  | pow(x,n) = F.*(x,pow(x,Int.-(n,1)))  
fun polyfit(x,y,n) =
    let open Array2
	val a = tabulate RowMajor (Array.length x,
				   Int.+(n,1),
				   fn (i,j) => if j=0 then F.one else
					       pow(Array.sub(x,i),j))
    in
	lsqr(a,y)
    end
end