I use LAPACK now, but long ago I wrote code for singular value
decomposition, which includes some code for eigenvalues/vectors of a
symmetric matrix. I have appended it here in case you would like to use it.
Henry Rich
display =: (i.0 0)"_ ((1!:2) 2:)
dispwithts =: 3 : 'display(": 3 }. 6!:0 ''),'' '',y'
NB. y is a symmetric matrix
NB. Result is the Householder tridiagonalization:
NB. diagonal elements ; subdiagonal elements ; xform =. |: orthonormal xform
NB. (Transpose of the usual form, for computational ease)
NB. (|: xform) mp diagonal mp xform reconstructs the original matrix
NB. Equations are from Numerical Recipes in C, 1st edition, pp. 368 ff.
NB. Henry H. Rich, 2001
trisymm =: 3 : 0
dispwithts 'starting trisymm'
diag =. offdiag =. acchu =. $0
m =. y
while. 2 < #m do.
hh =. }. {. m NB. 1st col of m, dropping first element
subm =. (< ,~ <<0) { m NB. m, dropping 1st row & column
mag =. %: ssq =. +/ *: hh NB. mag, magsq of hh
smag =. mag * (0 > {. hh) { 1 _1 NB. +- mag of hh.
NB. This produces the Householder matrix; we choose the sign
NB. to be the same as the sign of the first element to reduce
NB. roundoff error.
u =. (smag + {. hh) 0} hh NB. I - */~ normalized u will be the
NB. Householder matrix: first column of m,
NB. with the first element +- magnitude of
column
h =. % ssq + smag * {. hh NB. reciprocal of h: magsq of modified
column is the magsq of
NB. the unmodified column +- magnitude * 1st element
offdiag =. offdiag , - smag NB. The resulting offdiagonal value will
NB. be -+ mag of unmodified column
p =. h * subm +/@:*"1 u NB. Eqn. 11.2.10
k =. -: h * +/ p * u NB. Eqn 11.2.11
q =. p - k * u NB. Eqn 11.2.12
diag =. diag , (< ,~ <0) { m NB. top-left element stays on the diagonal
m =. subm - (q */ u) + (u */ q) NB. Eqn. 11.2.13
acchu =. acchu , < h , u NB. Accumulate h and u for rollup at end
end.
diag =. diag , (<0 1) |: m
offdiag =. offdiag , (< 0 ; <<0) { m
NB. x is h,u y is matrix #m = #u
NB. Result is M * (Householder matrix derived from h,u), with identity
column & row appended
rollup =. 13 : '(1 , (#y)$0) , 0 ,. y - (y (+/@:*)"1 (}.x)) */ ((}. *
{.) x)'
diag ; offdiag ; rollup&.>/ acchu , < =/~ i. 2
)
NB. y is diagonal elements ; subdiagonal elements ; |: orthonormal xform
NB. Result is eigenvalues ; |: eigenvectors (i. e. 1-cells are eigenvectors)
NB. This is the QL algorithm with implicit shifts, taken from
NB. Numerical Recipes in C pp. 380-381
eigentrisymm =: 3 : 0
dispwithts 'starting eigentrisymm'
NB. We extend the offdiagonal elements with one 0 to make them
NB. match the diagonal and to provide a stopper for the search below.
d =. 0{::y
e =. (1{::y),0
v =. 2{::y NB. Init eigenvectors to xform to date
NB. Loop to produce each eigenvalue. We consider all elements starting at
NB. 'start' - when we compute one eigenvalue, we know that we have
rotated the offdiagonal
NB. start{e to zero, so we put the eigenvalue over start{d, advance
start, and continue
NB. processing the (now smaller) matrix. This is the 'deflation' of d&e.
start =. 0
while. start < #d do.
NB. Find a small off-diagonal element (i. e. zero to machine
precision when compared
NB. with the diagonal elements); the iteration can start there, since
a rotation
NB. of the row containing that offdiagonal will leave it zero.
NB. If the small off-diagonal element is the first, stop iterating
and out
NB. the eigenvalue
while. start < i =. split =. start + ( e ((] =!.0 +) (+
1&(|.!.0))@:|)&(start&}.) d ) i. 1 do.
NB. Calculate the shift for the implicit-shift method. di1 is used
in the loop as
NB. (i+1){d
g =. ( (di1 =. split{d) - (start{d) ) + (-/ (1 0+start){d) ( ] % (+
((* 0&<:)~ >:&.*:))@(% +:) ) start{e
s =. c =. 1 NB. Init the rotation
p =. 0
while. (i =. <: i) >: start do.
NB. We are now performing Givens rotations on rows i and i+1, by
the book
f =. s * i{e
b =. c * i{e
if. (f >:&| g) do.
e =. (f*r =. %: >: *: c =. g % f) (>:i)} e
c =. c * s =. %r
else.
e =. (g*r =. %: >: *: s =. f % g) (>:i)} e
s =. s * c =. %r
end.
d =. (g + p =. s * r =. ( ( (di1 =. i{d) - g =. di1 - p ) * s ) +
+: c * b) (>:i)} d
g =. (c*r) - b
NB. The rotation was given by c and s. Rotate the eigenvectors
correspondingly
v =. (((c,-s),:s,c) +/ . * (0 1+i) { v) (0 1+i)} v
end.
d =. ( (start{d) - p ) start} d NB. Install the calculated eigenvalue
e =. (g,0) (start,split)} e NB. Install the first offdiagonal, and
clear the last, which was trashed
end.
start =. >: start
end.
d ; v
)
NB. Return eigenvalues and eigenvectors of symmetric matrix
NB. y is symmetric matrix
NB. result is eigenvalues ; (|: eigenvectors)
symmeigeninfo =: eigentrisymm @: trisymm
NB. y is n,m$ matrix
NB. Result is y +/ . * |:y (the input matrix is transposed already,
NB. you see). We do this without an explicit transpose, and we take
advantage of the
NB. symmetricity of the result: compute upper half of result, then
reflect to bottom
NB. half
NB. calcata =: 13 : '> (([ , }...@[ ,. ])&.>)/ (<@((+/@:*)"1 {.))\. y'
afterdisplaying 'starting calcata'
NB. This version doesn't use symmetry, but it seems to use less memory
calcata =: +/@:*"1/~
NB. svd in transpose space
NB. y is |: matrix
NB. returns: (|: left eigenvectors),.(|: right eigenvectors),.singular
values
NB. so if the return is x;y;z, (|:x) +/ . * y * z is the original matrix
NB. rows of x are left eigenvectors
NB. rows of z are right eigenvectors
svd =: 3 : 0
NB. For best results y should have more columns than rows
dispwithts 'starting calcata'
'r s'=. symmeigeninfo calcata y
dispwithts 'Inverting r'
t=. %:r
dispwithts 'Calculating u and stitching results'
NB. y =. y;s;t
y =. (y +/ . * ~ s % t),.s,.t
if. 0: dispwithts 'Returning boxed results' do. end.
)
NB. =================================================
NB. test for random matrices of size y (default 3 3), or given matrix
NB. e.g. testsvd 4 5
testsvd=: 3 : 0
diagmat=. * =...@i.@#
if. 2=#$y do. mat=. y
else. mat=. ? (2{.y,3 3) $ 10 end.
x =. svd |: mat
z =. {:"1 x
y =. }:" 1 (c =. - >: {:$mat) {."1 x
x =. c }."1 x
NB. mat should be:
NB. t=. (|:x) +/ . * (diagmat z) +/ . * (y)
t=. (|:x) +/ . * y * (z)
ok1 =. *./ 0.001 > | , mat-t
NB. Furthermore, x +/ . * |:x and y +/ . * |: y should be identity
ok2 =. *./ 0.001 > | , (=/~ i. # x) - calcata x
ok3 =. *./ 0.001 > | , (=/~ i. # y) - calcata y
ok1 *. ok2 *. ok3
)
Igor Zhuravlov wrote:
> В сообщении от Четверг 10 июня 2010 13:34:57 вы написали:
>> Although it would be a lot of work, have you considered implementing any
>> of the eigenvalue routines, like those in dgeev? Even if it were just for
>> symmetric matrices it would be a big step forward.
>
> You are right, there is a lot of work here. But this task is not impossible.
> xGEEx (xGEEV and almost identical to it xGEES) code consists of approx. 7400
> Fortran lines (including comments and excluding similar routines). Approx.
> 20%
> of xGEEx is already implemented [1]. I think xGEEx will be implemented some
> day, by me or somebody else.
>
> As for me, I'm inclining to import the generalized EV routines xGGEx (xGGEV
> and similar xGGES) first (see TODO file in mt addon). xGGEx code consists of
> approx. 9200 Fortran lines (including comments and excluding duplicated
> code).
> But:
> - approx. 40% of xGGEx is already implemented [1];
> - xGGEV will allow me to advance in yet another addon (still in development).
>
> ---
> [1] Callgraphs for LAPACK eigenvalue-related routines are referenced in
> http://www.jsoftware.com/jwiki/Addons/math/mt/MATRIX
>
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
