44 lines
1.5 KiB
Text
44 lines
1.5 KiB
Text
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));
|