QR decomposition: Difference between revisions
Exact and arbitrary precision solution
(D entry updated) |
(Exact and arbitrary precision solution) |
||
Line 1,196:
/* Note: the lapack package is a lisp translation of the fortran lapack library */</lang>
For an exact or arbitrary precision solution:<lang maxima>load("linearalgebra")$
load("eigen")$
unitVector(n) := ematrix(n,1,1,1,1);
signValue(r) := block([s:sign(r)],
if s='pos then 1 else if s='zero then 0 else -1);
householder(a) := block([m : length(a),u,v,beta],
u : a + sqrt(a . a)*signValue(a[1,1])*unitVector(m),
v : u / u[1,1],
beta : 2/(v . v),
diagmatrix(m,1) - beta*transpose(v . transpose(v)));
getSubmatrix(obj,i1,j1,i2,j2) :=
genmatrix(lambda([i,j], obj[i+i1-1,j+j1-1]),i2-i1+1,j2-j1+1);
setSubmatrix(obj,i1,j1,subobj) := block([m,n],
[m,n] : matrix_size(subobj),
for i: 0 thru m-1 do
(for j: 0 thru n-1 do
obj[i1+i,j1+j] : subobj[i+1,j+1]));
qr(obj) := block([m,n,qm,rm,i],
[m,n] : matrix_size(obj),
qm : diagmatrix(m,1),
rm : copymatrix(obj),
for i: 1 thru (if m=n then n-1 else n) do
block([x,h],
x : getSubmatrix(rm,i,i,m,i),
h : diagmatrix(m,1),
setSubmatrix(h,i,i,householder(x)),
qm : qm . h,
rm : h . rm),
[qm,rm]);
solveUpperTriangular(r,b) := block([n,x,index,k],
n : second(matrix_size(r)),
x : genmatrix(lambda([a, b], 0), n, 1),
for k: n thru 1 step -1 do
(index : min(n,k+1),
x[k,1] : (b[k,1] - (getSubmatrix(r,k,index,k,n) . getSubmatrix(x,index,1,n,1)))/r[k,k]),
x);
lsqr(a,b) := block([q,r,n],
[q,r] : qr(a),
n : second(matrix_size(r)),
solveUpperTriangular(getSubmatrix(r,1,1,n,n), transpose(q) . b));
polyfit(x,y,n) := block([a,j],
a : genmatrix(lambda([i,j], if j=1 then 1.0b0 else bfloat(x[i,1]^(j-1))),
length(x),n+1),
lsqr(a,y));</lang>Then we have the examples:<lang maxima>(%i) [q,r] : qr(a);
[ 6 69 58 ]
[ - - --- --- ]
[ 7 175 175 ]
[ ] [ - 14 - 21 14 ]
[ 3 158 6 ] [ ]
(%o) [[ - - - --- - --- ], [ 0 - 175 70 ]]
[ 7 175 175 ] [ ]
[ ] [ 0 0 - 35 ]
[ 2 6 33 ]
[ - - -- -- ]
[ 7 35 35 ]
(%i) mat_norm(q . r - a, 1);
(%o) 0
(%i) x : transpose(matrix([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]))$
(%i) y : transpose(matrix([1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321]))$
(%i) fpprec : 30$
(%i) polyfit(x, y, 2);
[ 9.99999999999999999999999999996b-1 ]
[ ]
(%o) [ 2.00000000000000000000000000002b0 ]
[ ]
[ 3.0b0 ]</lang>
=={{header|PARI/GP}}==
| |||