I know this has been brought up several times, but I have to bring it up
yet again...
I am trying to find an algorithm or an implementation for calculating the
eigenvalues and eigenvectors of a symmetric matrix. If anyone has any
pointers, I would appreciate it. The next best alternative for me is to use
MATLAB, but I doubt constructing an interface for MATLAB in lisp is going
to be easy.
Thanks!
Sunil
Sunil Mishra wrote in message ...
>I know this has been brought up several times, but I have to bring it up
>yet again...
>
>I am trying to find an algorithm or an implementation for calculating the
>eigenvalues and eigenvectors of a symmetric matrix. If anyone has any
>pointers, I would appreciate it. The next best alternative for me is to use
>MATLAB, but I doubt constructing an interface for MATLAB in lisp is going
>to be easy.
>
>Thanks!
>
>Sunil
This is plagerism, I am sure. We generalized this fortran function from
_Numerical Recipes_ to in muLisp. The conventions of our generalization
are:
1. every vector space V will be finite-dimensional and over the complex
numbers.
2. matrix means a sqaure matrix with complex elements
3. every linear transfromation T will be from V into itself
4. If T is a linear transformation, then a representation of T in matix A
will be a representation of the type RxxT = A where x is a basis for V
(here xx is subscripts)
[the muLisp code is property of the company I work for thus I can not share
it.] I must be missing something in your request for there are several
routines for latent roots and linear root transformations out there perhaps
not in lisp but many in fortran.
If I have broken some copywrite laws, rules, etiquette I will repent.
btw there was a posting from Leo Saraua a few weeks back. He has a rather
good start on a lisp "solver", "algerba" and "matrix" math package. I down
loaded then from one of the groups web sites. It would be a very good
starting place. I don't know enough common lisp to conjure how difficult
inserting this function would be. I have used his excellent examples to
clean up some of our mulisp math functions. (hmmm this must explain why my
lisp looks like fortram)
SUBROUTINE jacobi(a,n,np,d,v,nrot)
INTEGER n,np,nrot,NMAX
REAL a(np,np),d(np),v(np,np)
PARAMETER (NMAX=500)
INTEGER i,ip,iq,j
REAL c,g,h,s,sm,t,tau,theta,tresh,b(NMAX),z(NMAX)
do 12 ip=1,n
do 11 iq=1,n
v(ip,iq)=0.
11 continue
v(ip,ip)=1.
12 continue
do 13 ip=1,n
b(ip)=a(ip,ip)
d(ip)=b(ip)
z(ip)=0.
13 continue
nrot=0
do 24 i=1,50
sm=0.
do 15 ip=1,n-1
do 14 iq=ip+1,n
sm=sm+abs(a(ip,iq))
14 continue
15 continue
if(sm.eq.0.)return
if(i.lt.4)then
tresh=0.2*sm/n**2
else
tresh=0.
endif
do 22 ip=1,n-1
do 21 iq=ip+1,n
g=100.*abs(a(ip,iq))
if((i.gt.4).and.(abs(d(ip))+
*g.eq.abs(d(ip))).and.(abs(d(iq))+g.eq.abs(d(iq))))then
a(ip,iq)=0.
else if(abs(a(ip,iq)).gt.tresh)then
h=d(iq)-d(ip)
if(abs(h)+g.eq.abs(h))then
t=a(ip,iq)/h
else
theta=0.5*h/a(ip,iq)
t=1./(abs(theta)+sqrt(1.+theta**2))
if(theta.lt.0.)t=-t
endif
c=1./sqrt(1+t**2)
s=t*c
tau=s/(1.+c)
h=t*a(ip,iq)
z(ip)=z(ip)-h
z(iq)=z(iq)+h
d(ip)=d(ip)-h
d(iq)=d(iq)+h
a(ip,iq)=0.
do 16 j=1,ip-1
g=a(j,ip)
h=a(j,iq)
a(j,ip)=g-s*(h+g*tau)
a(j,iq)=h+s*(g-h*tau)
16 continue
do 17 j=ip+1,iq-1
g=a(ip,j)
h=a(j,iq)
a(ip,j)=g-s*(h+g*tau)
a(j,iq)=h+s*(g-h*tau)
17 continue
do 18 j=iq+1,n
g=a(ip,j)
h=a(iq,j)
a(ip,j)=g-s*(h+g*tau)
a(iq,j)=h+s*(g-h*tau)
18 continue
do 19 j=1,n
g=v(j,ip)
h=v(j,iq)
v(j,ip)=g-s*(h+g*tau)
v(j,iq)=h+s*(g-h*tau)
19 continue
nrot=nrot+1
endif
21 continue
22 continue
do 23 ip=1,n
b(ip)=b(ip)+z(ip)
d(ip)=b(ip)
z(ip)=0.
23 continue
24 continue
pause 'too many iterations in jacobi'
return
END
rusty craine wrote in message <············@excalibur.flash.net>...
>
>If I have broken some copywrite laws, rules, etiquette I will repent.
>
Oh well since I am most likely in trouble with the copywrite police anyway
her is the _numercial recipes_ demo driver program for the latent root
function.
rusty craine wrote:
> >If I have broken some copywrite laws, rules, etiquette I will repent.
> Oh well since I am most likely in trouble with the copywrite police anyway
At least theoretically you are (unless the function you posted
is one of the exceptions), since the NumRep copyright is rather
ugly (the preface of the book disallows even the use of
derivative work for other purposes than one's own private
education).
I have no idea if anybody cares -- I just remember somebody long
ago replacing a FFT modelled after what was in NumRep after
learning about the copyright.
Ralf
rusty craine wrote in message <············@excalibur.flash.net>...
>
>
sorry I am going to have to learn to use a mouse someday :)
now lets try this one more time....sorry for posting fortran in a lisp news
group
(btw in math circles _numerical recipes_ has been critised for algorithms
and codeing. None the less it works for us)
PROGRAM xjacobi
C driver for routine jacobi
INTEGER NP,NMAT
PARAMETER(NP=10,NMAT=3)
INTEGER i,ii,j,jj,k,kk,l,ll,nrot,num(3)
REAL ratio,d(NP),v(NP,NP),r(NP)
REAL a(3,3),b(5,5),c(10,10),e(NP,NP)
DATA num/3,5,10/
DATA a/1.0,2.0,3.0,2.0,2.0,3.0,3.0,3.0,3.0/
DATA b/-2.0,-1.0,0.0,1.0,2.0,-1.0,-1.0,0.0,1.0,2.0,
* 0.0,0.0,0.0,1.0,2.0,1.0,1.0,1.0,1.0,2.0,
* 2.0,2.0,2.0,2.0,2.0/
DATA c/5.0,4.3,3.0,2.0,1.0,0.0,-1.0,-2.0,-3.0,-4.0,
* 4.3,5.1,4.0,3.0,2.0,1.0,0.0,-1.0,-2.0,-3.0,
* 3.0,4.0,5.0,4.0,3.0,2.0,1.0,0.0,-1.0,-2.0,
* 2.0,3.0,4.0,5.0,4.0,3.0,2.0,1.0,0.0,-1.0,
* 1.0,2.0,3.0,4.0,5.0,4.0,3.0,2.0,1.0,0.0,
* 0.0,1.0,2.0,3.0,4.0,5.0,4.0,3.0,2.0,1.0,
* -1.0,0.0,1.0,2.0,3.0,4.0,5.0,4.0,3.0,2.0,
* -2.0,-1.0,0.0,1.0,2.0,3.0,4.0,5.0,4.0,3.0,
* -3.0,-2.0,-1.0,0.0,1.0,2.0,3.0,4.0,5.0,4.0,
* -4.0,-3.0,-2.0,-1.0,0.0,1.0,2.0,3.0,4.0,5.0/
do 24 i=1,NMAT
if (i.eq.1) then
do 12 ii=1,3
do 11 jj=1,3
e(ii,jj)=a(ii,jj)
11 continue
12 continue
call jacobi(e,3,NP,d,v,nrot)
else if (i.eq.2) then
do 14 ii=1,5
do 13 jj=1,5
e(ii,jj)=b(ii,jj)
13 continue
14 continue
call jacobi(e,5,NP,d,v,nrot)
else if (i.eq.3) then
do 16 ii=1,10
do 15 jj=1,10
e(ii,jj)=c(ii,jj)
15 continue
16 continue
call jacobi(e,10,NP,d,v,nrot)
endif
write(*,'(/1x,a,i2)') 'Matrix Number ',i
write(*,'(1x,a,i3)') 'Number of JACOBI rotations: ',nrot
write(*,'(/1x,a)') 'Eigenvalues:'
do 17 j=1,num(i)
write(*,'(1x,5f12.6)') d(j)
17 continue
write(*,'(/1x,a)') 'Eigenvectors:'
do 18 j=1,num(i)
write(*,'(1x,t5,a,i3)') 'Number',j
write(*,'(1x,5f12.6)') (v(k,j),k=1,num(i))
18 continue
C eigenvector test
write(*,'(/1x,a)') 'Eigenvector Test'
do 23 j=1,num(i)
do 21 l=1,num(i)
r(l)=0.0
do 19 k=1,num(i)
if (k.gt.l) then
kk=l
ll=k
else
kk=k
ll=l
endif
if (i.eq.1) then
r(l)=r(l)+a(ll,kk)*v(k,j)
else if (i.eq.2) then
r(l)=r(l)+b(ll,kk)*v(k,j)
else if (i.eq.3) then
r(l)=r(l)+c(ll,kk)*v(k,j)
endif
19 continue
21 continue
write(*,'(/1x,a,i3)') 'Vector Number',j
write(*,'(/1x,t7,a,t18,a,t31,a)')
* 'Vector','Mtrx*Vec.','Ratio'
do 22 l=1,num(i)
ratio=r(l)/v(l,j)
write(*,'(1x,3f12.6)') v(l,j),r(l),ratio
22 continue
23 continue
write(*,*) 'press RETURN to continue...'
read(*,*)
24 continue
END
Sunil Mishra <·······@hustle.cc.gatech.edu> writes:
> I know this has been brought up several times, but I have to bring it up
> yet again...
>
> I am trying to find an algorithm or an implementation for calculating the
> eigenvalues and eigenvectors of a symmetric matrix. If anyone has any
> pointers, I would appreciate it. The next best alternative for me is to use
> MATLAB, but I doubt constructing an interface for MATLAB in lisp is going
> to be easy.
>
> Thanks!
>
> Sunil
Thanks a lot for your pointers, all! I have decided on a solution. I have
started working on FFI's to LAPACK. Now if only I can figure out which
libraries I need to include when loading LAPACK object code into lispworks
on IRIX 6.3...
Sunil
From: Christopher Stacy
Subject: Re: Matrix computations in lisp?
Date:
Message-ID: <uvhkt70ba.fsf@pilgrim.com>
MACSYMA includes a MATLAB translator and implements the standard libraries
(as well, of course, as many other things that MATLAB can't do...)
Sunil Mishra wrote:
> I know this has been brought up several times, but I have to bring it up
> yet again...
>
> I am trying to find an algorithm or an implementation for calculating the
> eigenvalues and eigenvectors of a symmetric matrix. If anyone has any
> pointers, I would appreciate it. The next best alternative for me is to use
> MATLAB, but I doubt constructing an interface for MATLAB in lisp is going
> to be easy.
XLispStat does what you want, and is free. You can get it from
http://www.stat.ucla.edu/develop/lisp/xlisp/xlisp-stat/
-----------------------
XLISP-PLUS version 3.0
Portions Copyright (c) 1988, by David Betz.
Modified by Thomas Almy and others.
XLISP-STAT Release 3.50 (Beta).
Copyright (c) 1989-1994, by Luke Tierney.
Initialization may take a moment.
> (help 'eigenvalues)
loading in help file information - this will take a minute ...done
EIGENVALUES [function-doc]
Args: (a)
Returns list of eigenvalues of square, symmetric matrix A
NIL
> (help 'eigenvectors)
EIGENVECTORS [function-doc]
Args: (a)
Returns list of eigenvectors of square, symmetric matrix A
NIL
>
>
>
> Thanks!
>
> Sunil
--
Le Hibou (mo bheachd fh�in: my own opinion)
"it's just that in C++ and the like, you don't trust _anybody_,
and in CLOS you basically trust everybody. the practical result
is that thieves and bums use C++ and nice people use CLOS."
-- Erik Naggum
> Sunil Mishra wrote:
>
> > I know this has been brought up several times, but I have to bring it up
> > yet again...
> >
> > I am trying to find an algorithm or an implementation for calculating the
> > eigenvalues and eigenvectors of a symmetric matrix. If anyone has any
> > pointers, I would appreciate it. The next best alternative for me is to use
> > MATLAB, but I doubt constructing an interface for MATLAB in lisp is going
> > to be easy.
I have translated a number of routines from Numerical Recipes in C into R4RS
Scheme: JACOBI, SIMPLEX, BRENT, MNBRAK, LINMIN, and POWELL. They are part of
my QobiScheme package, available free from my home page. Because you have
specifically asked for eigenvalue/eigenvector computations, I've enclosed
those below.
Jeff (http://www.neci.nj.nec.com/homepages/qobi)
-------------------------------------------------------------------------------
(define (jacobi a)
(unless (and (= (matrix-rows a) (matrix-columns a))
(every-n (lambda (i)
(every-n (lambda (j)
(= (matrix-ref a i j) (matrix-ref a j i)))
(matrix-rows a)))
(matrix-rows a)))
(panic "Can only compute eigenvalues/eigenvectors of a symmetric matrix"))
(let* ((a (map-vector (lambda (row) (map-vector identity row)) a))
(n (matrix-rows a))
(d (make-vector n))
(v (make-matrix n n 0.0))
(b (make-vector n))
(z (make-vector n 0.0)))
(for-each-n (lambda (ip)
(matrix-set! v ip ip 1.0)
(vector-set! b ip (matrix-ref a ip ip))
(vector-set! d ip (matrix-ref a ip ip)))
n)
(let loop ((i 0))
(when (> i 50) (panic "Too many iterations in JACOBI"))
(let ((sm (sum (lambda (ip)
(sum (lambda (ir)
(let ((iq (+ ip ir 1)))
(abs (matrix-ref a ip iq))))
(- n ip 1)))
(- n 1))))
(unless (zero? sm)
(let ((tresh (if (< i 3) (/ (* 0.2 sm) (* n n)) 0.0)))
(for-each-n
(lambda (ip)
(for-each-n
(lambda (ir)
(let* ((iq (+ ip ir 1))
(g (* 100.0 (abs (matrix-ref a ip iq)))))
(cond
((and (> i 3)
(= (+ (abs (vector-ref d ip)) g) (abs (vector-ref d ip)))
(= (+ (abs (vector-ref d iq)) g) (abs (vector-ref d iq))))
(matrix-set! a ip iq 0.0))
((> (abs (matrix-ref a ip iq)) tresh)
(let* ((h (- (vector-ref d iq) (vector-ref d ip)))
(t (if (= (+ (abs h) g) (abs h))
(/ (matrix-ref a ip iq) h)
(let ((theta (/ (* 0.5 h) (matrix-ref a ip iq))))
(if (negative? theta)
(- (/ (+ (abs theta)
(sqrt (+ (* theta theta) 1.0)))))
(/ (+ (abs theta)
(sqrt (+ (* theta theta) 1.0))))))))
(c (/ (sqrt (+ (* t t) 1.0))))
(s (* t c))
(tau (/ s (+ c 1.0)))
(h (* t (matrix-ref a ip iq))))
(define (rotate a i j k l)
(let ((g (matrix-ref a i j))
(h (matrix-ref a k l)))
(matrix-set! a i j (- g (* s (+ h (* g tau)))))
(matrix-set! a k l (+ h (* s (- g (* h tau)))))))
(vector-set! z ip (- (vector-ref z ip) h))
(vector-set! z iq (+ (vector-ref z iq) h))
(vector-set! d ip (- (vector-ref d ip) h))
(vector-set! d iq (+ (vector-ref d iq) h))
(matrix-set! a ip iq 0.0)
(for-each-n (lambda (j)
(cond ((< j ip) (rotate a j ip j iq))
((< ip j iq) (rotate a ip j j iq))
((< iq j) (rotate a ip j iq j)))
(rotate v j ip j iq))
n))))))
(- n ip 1)))
(- n 1)))
(for-each-n (lambda (ip)
(vector-set! b ip (+ (vector-ref b ip) (vector-ref z ip)))
(vector-set! d ip (vector-ref b ip))
(vector-set! z ip 0.0))
n)
(loop (+ i 1)))))
(for-each-n
(lambda (i)
(let ((k i)
(p (vector-ref d i)))
(for-each-n
(lambda (l)
(let* ((j (+ i l 1)))
(when (>= (vector-ref d j) p) (set! k j) (set! p (vector-ref d j)))))
(- n i 1))
(unless (= k i)
(vector-set! d k (vector-ref d i))
(vector-set! d i p)
(for-each-n (lambda (j)
(let ((p (matrix-ref v j i)))
(matrix-set! v j i (matrix-ref v j k))
(matrix-set! v j k p)))
n))))
(- n 1))
(list d (transpose v))))
(define (eigenvalues a) (first (jacobi a)))
(define (eigenvectors a) (second (jacobi a)))
(define (vector->diagonal-matrix v)
(let ((m (make-matrix (vector-length v) (vector-length v) 0.0)))
(for-each-n (lambda (i) (matrix-set! m i i (vector-ref v i)))
(vector-length v))
m))
(define (identity-matrix n) (vector->diagonal-matrix (make-vector n 1.0)))
(define (clip-eigenvalues a v)
(let* ((j (jacobi a))
(e (second j)))
(m* (transpose e)
(m* (vector->diagonal-matrix (map-vector max v (first j))) e))))
;;; The following two routines are limited to 2-by-2 matricies.
(define (eigenvector-angle1 m)
(if (and (< (abs (matrix-ref m 1 0)) *epsilon*)
(< (abs (matrix-ref m 0 1)) *epsilon*))
(if (> (matrix-ref m 1 1) (matrix-ref m 0 0)) half-pi 0.0)
(atan (matrix-ref m 1 0)
(- (vector-ref (eigenvalues m) 0) (matrix-ref m 1 1)))))
(define (eigenvector-angle2 m)
(if (and (< (abs (matrix-ref m 1 0)) *epsilon*)
(< (abs (matrix-ref m 0 1)) *epsilon*))
(if (<= (matrix-ref m 1 1) (matrix-ref m 0 0)) half-pi 0.0)
(atan (matrix-ref m 1 0)
(- (vector-ref (eigenvalues m) 1) (matrix-ref m 1 1)))))
>>>>> "Sunil" == Sunil Mishra <·······@hustle.cc.gatech.edu> writes:
Sunil> I know this has been brought up several times, but I have to bring it up
Sunil> yet again...
Sunil> I am trying to find an algorithm or an implementation for calculating the
Sunil> eigenvalues and eigenvectors of a symmetric matrix. If anyone has any
Sunil> pointers, I would appreciate it. The next best alternative for me is to use
Sunil> MATLAB, but I doubt constructing an interface for MATLAB in lisp is going
Sunil> to be easy.
You can, of course, use a foreign function interface to LAPACK. I've
done this with CMUCL and LAPACK. Works nicely.
Ray