How to resolve the algorithm QR decomposition step by step in the Rascal programming language
How to resolve the algorithm QR decomposition step by step in the Rascal 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 Rascal programming language
Source code in the rascal programming language
import util::Math;
import Prelude;
import vis::Figure;
import vis::Render;
public rel[real,real,real] QRdecomposition(rel[real x, real y, real v] matrix){
//orthogonalcolumns
oc = domainR(matrix, {0.0});
for (x <- sort(toList(domain(matrix)-{0.0}))){
c = domainR(matrix, {x});
o = domainR(oc, {x-1});
for (n <- [1.0 .. x]){
o = domainR(oc, {n-1});
c = matrixSubtract(c, matrixMultiplybyN(o, matrixDotproduct(o, c)/matrixDotproduct(o, o)));
}
oc += c;
}
Q = {};
//from orthogonal to orthonormal columns
for (el <- oc){
c = domainR(oc, {el[0]});
Q += matrixNormalize({el}, c);
}
//from Q to R
R= matrixMultiplication(matrixTranspose(Q), matrix);
R= { | <- R};
println("Q:");
iprintlnExp(Q);
println();
println("R:");
return R;
}
//a function that takes the transpose of a matrix, see also Rosetta Code problem "Matrix transposition"
public rel[real, real, real] matrixTranspose(rel[real x, real y, real v] matrix){
return { | <- matrix};
}
//a function to normalize an element of a matrix by the normalization of a column
public rel[real,real,real] matrixNormalize(rel[real x, real y, real v] element, rel[real x, real y, real v] column){
normalized = 1.0/nroot((0.0 | it + v*v | <- column), 2);
return matrixMultiplybyN(element, normalized);
}
//a function that takes the dot product, see also Rosetta Code problem "Dot product"
public real matrixDotproduct(rel[real x, real y, real v] column1, rel[real x, real y, real v] column2){
return (0.0 | it + v1*v2 | <- column1, <- column2, y1==y2);
}
//a function to subtract two columns
public rel[real,real,real] matrixSubtract(rel[real x, real y, real v] column1, rel[real x, real y, real v] column2){
return { | <- column1, <- column2, y1==y2};
}
//a function to multiply a column by a number
public rel[real,real,real] matrixMultiplybyN(rel[real x, real y, real v] column, real n){
return { | <- column};
}
//a function to perform matrix multiplication, see also Rosetta Code problem "Matrix multiplication".
public rel[real, real, real] matrixMultiplication(rel[real x, real y, real v] matrix1, rel[real x, real y, real v] matrix2){
if (max(matrix1.x) == max(matrix2.y)){
p = { | <- matrix1, <- matrix2};
result = {};
for (y <- matrix1.y){
for (x <- matrix2.x){
v = (0.0 | it + v | <- p, x==x2 && y==y1, x1==y2 && y2==x1);
result += ;
}
}
return result;
}
else throw "Matrix sizes do not match.";
}
// a function to visualize the result
public void displayMatrix(rel[real x, real y, real v] matrix){
points = [box(text(""), align(0.3333*(x+1),0.3333*(y+1)),shrink(0.25)) | <- matrix];
render(overlay([*points], aspectRatio(1.0)));
}
//a matrix, given by a relation of .
public rel[real x, real y, real v] matrixA = {
<0.0,0.0,12.0>, <0.0,1.0, 6.0>, <0.0,2.0,-4.0>,
<1.0,0.0,-51.0>, <1.0,1.0,167.0>, <1.0,2.0,24.0>,
<2.0,0.0,4.0>, <2.0,1.0,-68.0>, <2.0,2.0,-41.0>
};
You may also check:How to resolve the algorithm Keyboard input/Obtain a Y or N response step by step in the QB64 programming language
You may also check:How to resolve the algorithm Read entire file step by step in the Zig programming language
You may also check:How to resolve the algorithm Comma quibbling step by step in the F# programming language
You may also check:How to resolve the algorithm Sierpinski carpet step by step in the Mathematica/Wolfram Language programming language
You may also check:How to resolve the algorithm FizzBuzz step by step in the hexiscript programming language