1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
|
comment
Single Value Decomposition Tools
endcomment
function svdkernel (A)
## kernel(a) computes the kernel of the quadratic matrix a.
## This is using the singular value decomposition and works
## for real matrices only.
## The vectors spanning the kernel are orthonormal.
{B,w,V}=svd(A);
i=nonzeros(w~=0);
if cols(i)==0; return zeros(cols(V),1);
else return V[:,i];
endif;
endfunction
function svdimage (A)
## image(a) computes the image of the quadratic matrix a.
## This is using the singular value decomposition and works
## for real matrices only.
## The vectors spanning the image are orthonormal.
{B,w,V}=svd(A);
i=nonzeros(!(w~=0));
if cols(i)==0; return zeros(1,cols(A));
else return B[:,i];
endif;
endfunction
function svdcondition (A)
## condition(A) returns the condition number based on
## a singular value decompostion of A.
## 0 means singularity.
## A must be real.
{B,w,V}=svd(A);
if any(w~=0); return 0; endif;
return max(abs(w))/min(abs(w));
endfunction
function svddet (A)
## det(A) returns the determinant based on a
## singular value decomposition of A.
## A must be real.
{B,w,V}=svd(A);
return prod(w);
endfunction
function svdeigenvalues (A)
## For a symmetrical A, this returns the eigenvalues of A.
## For a non-symmetrical A, the singular values.
## A must be real.
{B,w,V}=svd(A);
return w;
endfunction
function svdeigenspace (a,l)
## svdeigenspace(A,l) returns the eigenspace of A to the eigenvaue l.
k=svdkernel(a-l*id(cols(a)));
if k==0; error("No eigenvalue found!"); endif;
si=size(k);
loop 1 to si[2];
x=k[:,index()];
k[:,index()]=x/sqrt(x'.x);
end;
return k;
endfunction
function svdsolve (A,b)
## Solve A.x=b with smallest norm for x.
## A must be real.
## This function can be used instead of the fit function.
{B,w,V}=svd(A);
i=nonzeros(!(w~=0));
if (cols(i)>0); w[i]=1/w[i]; endif;
return V.diag(size(V),0,w).(B'.b);
endfunction
|
