How to resolve the algorithm QR decomposition step by step in the Axiom programming language
How to resolve the algorithm QR decomposition step by step in the Axiom programming language
Table of Contents
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:
Step by Step solution about How to resolve the algorithm QR decomposition step by step in the Axiom programming language
Source code in the axiom programming language
)abbrev package TESTP TestPackage
TestPackage(R:Join(Field,RadicalCategory)): with
unitVector: NonNegativeInteger -> Vector(R)
"/": (Vector(R),R) -> Vector(R)
"^": (Vector(R),NonNegativeInteger) -> Vector(R)
solveUpperTriangular: (Matrix(R),Vector(R)) -> Vector(R)
signValue: R -> R
householder: Vector(R) -> Matrix(R)
qr: Matrix(R) -> Record(q:Matrix(R),r:Matrix(R))
lsqr: (Matrix(R),Vector(R)) -> Vector(R)
polyfit: (Vector(R),Vector(R),NonNegativeInteger) -> Vector(R)
== add
unitVector(dim) ==
out := new(dim,0@R)$Vector(R)
out(1) := 1@R
out
v:Vector(R) / a:R == map((vi:R):R +-> vi/a, v)$Vector(R)
v:Vector(R) ^ n:NonNegativeInteger == map((vi:R):R +-> vi^n, v)$Vector(R)
solveUpperTriangular(r,b) ==
n := ncols r
x := new(n,0@R)$Vector(R)
for k in n..1 by -1 repeat
index := min(n,k+1)
x(k) := (b(k)-reduce("+",subMatrix(r,k,k,index,n)*x.(index..n)))/r(k,k)
x
signValue(r) ==
R has (sign: R -> Integer) => coerce(sign(r)$R)$R
zero? r => r
if sqrt(r*r) = r then 1 else -1
householder(a) ==
m := #a
u := a + length(a)*signValue(a(1))*unitVector(m)
v := u/u(1)
beta := (1+1)/dot(v,v)
scalarMatrix(m,1) - beta*transpose(outerProduct(v,v))
qr(a) ==
(m,n) := (nrows a, ncols a)
qm := scalarMatrix(m,1)
rm := copy a
for i in 1..(if m=n then n-1 else n) repeat
x := column(subMatrix(rm,i,m,i,i),1)
h := scalarMatrix(m,1)
setsubMatrix!(h,i,i,householder x)
qm := qm*h
rm := h*rm
[qm,rm]
lsqr(a,b) ==
dc := qr a
n := ncols(dc.r)
solveUpperTriangular(subMatrix(dc.r,1,n,1,n),transpose(dc.q)*b)
polyfit(x,y,n) ==
a := new(#x,n+1,0@R)$Matrix(R)
for j in 0..n repeat
setColumn!(a,j+1,x^j)
lsqr(a,y)
m := matrix [[12, -51, 4], [6, 167, -68], [-4, 24, -41]];
qr m
x := vector [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10];
y := vector [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321];
polyfit(x, y, 2)
qr m
+ 6 69 58 +
|- - --- --- |
| 7 175 175 |
| | +- 14 - 21 14 +
| 3 158 6 | | |
[q= |- - - --- - ---|,r= | 0 - 175 70 |]
| 7 175 175| | |
| | + 0 0 - 35+
| 2 6 33 |
| - - -- -- |
+ 7 35 35 +
Type: Record(q: Matrix(AlgebraicNumber),r: Matrix(AlgebraicNumber))
polyfit(x, y, 2)
[1,2,3]
Type: Vector(AlgebraicNumber)
You may also check:How to resolve the algorithm Sum digits of an integer step by step in the REXX programming language
You may also check:How to resolve the algorithm Flipping bits game step by step in the AutoHotkey programming language
You may also check:How to resolve the algorithm Arbitrary-precision integers (included) step by step in the Lua programming language
You may also check:How to resolve the algorithm Miller–Rabin primality test step by step in the Elixir programming language
You may also check:How to resolve the algorithm Harshad or Niven series step by step in the Java programming language