Anonymous user
QR decomposition: Difference between revisions
no edit summary
(Updated D entry) |
No edit summary |
||
Line 1,072:
m <- lm(y ~ x + xx)
coef(m)</lang>
=={{header|Rascal}}==
[[File:qrresult.jpeg||200px|thumb|right]]
This function applies the Gram Schmidt algorithm. Q is printed in the console, R can be printed or visualized.
<lang Rascal>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= {<x,y,toReal(round(v))> | <x,y,v> <- 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 {<y, x, v> | <x, y, v> <- 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 | <x,y,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 | <x1,y1,v1> <- column1, <x2,y2,v2> <- 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 {<x1,y1,v1-v2> | <x1,y1,v1> <- column1, <x2,y2,v2> <- 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 {<x,y,v*n> | <x,y,v> <- 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 = {<x1,y1,x2,y2, v1*v2> | <x1,y1,v1> <- matrix1, <x2,y2,v2> <- matrix2};
result = {};
for (y <- matrix1.y){
for (x <- matrix2.x){
v = (0.0 | it + v | <x1, y1, x2, y2, v> <- p, x==x2 && y==y1, x1==y2 && y2==x1);
result += <x,y,v>;
}
}
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("<v>"), align(0.3333*(x+1),0.3333*(y+1)),shrink(0.25)) | <x,y,v> <- matrix];
render(overlay([*points], aspectRatio(1.0)));
}
//a matrix, given by a relation of <x-coordinate, y-coordinate, value>.
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>
};</lang>
Example using visualization
<pre>rascal>displayMatrix(QRdecomposition(matrixA))
Q:
{
<1.0,0.0,-0.394285714285714285714285714285714285714285714285714285714285714285713300>,
<0.0,0.0,0.857142857142857142857142857142857142857142857142857142857142857142840>,
<0.0,1.0,0.428571428571428571428571428571428571428571428571428571428571428571420>,
<0.0,2.0,-0.285714285714285714285714285714285714285714285714285714285714285714280>,
<2.0,0.0,-0.33142857142857142857142857142857142857142857142857142857142857142858800>,
<1.0,2.0,0.171428571428571428571428571428571428571428571428571428571428571428571000>,
<2.0,2.0,-0.94285714285714285714285714285714285714285714285714285714285714285719000>,
<1.0,1.0,0.902857142857142857142857142857142857142857142857142857142857142857140600>,
<2.0,1.0,0.03428571428571428571428571428571428571428571428571428571428571428571600>
}
See R in picture</pre>
=={{header|Tcl}}==
| |||