QR decomposition: Difference between revisions
Content added Content deleted
(Updated D entry) |
No edit summary |
||
Line 1,072: | Line 1,072: | ||
m <- lm(y ~ x + xx) |
m <- lm(y ~ x + xx) |
||
coef(m)</lang> |
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}}== |
=={{header|Tcl}}== |