Subject: Segfault 5.9.1 on large matrix operations
From: Daniel Lakeland
Date: Fri, 25 Nov 2005 10:47:37 -0800
I have a surveying traverse with variable precision in measurements.
I was trying to do a constrainted least squares optimization of the
correction factors to make the closed loop actually closed using
lagrange multipliers.
Essentially what I have is a weighted least squares objective
function, and 2 constraints (total delta x = 0 and total delta y = 0).
--------------------------- code below ----------
vars: append(makelist(concat(da,i),i,1,length(dxes)),
makelist(concat(dl,i),i,1,length(dxes)));
traveqn1: sum((dists[i]+concat(dl,i))*cos(orientsrad[i]+concat(da,i)),i,1,length(dxes));
traveqn2: sum((dists[i]+concat(dl,i))*sin(orientsrad[i]+concat(da,i)),i,1,length(dyes));
objective: sum(1/wa[i]^2*concat(da,i)^2,i,1,length(dxes))
+ sum(1/dsigs[i]^2*concat(dl,i)^2,i,1,length(dxes));
/* here are my lagrange equations, l1 and l2 are the multipliers*/
objgrad1:append(makelist(diff(objective-l1*traveqn1,concat(dl,i)),i,1,nvars),
makelist(diff(objective-l1*traveqn1,concat(da,i)),i,1,nvars));
objgrad2:append(makelist(diff(objective-l2*traveqn2,concat(dl,i)),i,1,nvars),
makelist(diff(objective-l2*traveqn2,concat(da,i)),i,1,nvars));
/*map(assume,makelist(concat(dl,i) > 0,i,1,nvars));*/
lagrsystem: append(makelist(objgrad1[i],i,1,nvars),
makelist(objgrad2[i],i,1,nvars),[traveqn1,traveqn2]);
load("mymnewton");
/*optimize(mnewton);*/
NEWTONEPSILON:1.0e-6;
NEWTONMAXITER:5;
sparse:true;
mnewton(float(lagrsystem),
append(flatten(makelist([concat(da,i),concat(dl,i)],i,1,nvars)),
[l1,l2]),
append(makelist(.001,i,1,nvars),makelist(.3,i,1,nvars),[10,10]));
-------------------------
It turns out there's a contributed package lsquares which does
basically the same thing (that is, calls mnewton).
mymnewton is just the standard mnewton package with a few print
statements to give me feedback. It starts the first iteration but
never finishes it.
I ran it last night and when I woke up this morning it had segfaulted
(5.9.1 on debian without the maxima-share package... I upgraded this
morning to 5.9.2 with the share package)
It used a ridiculous amount of RAM (the method is not efficient) but
it didn't have a problem allocating (on my laptop it complained "could
not allocate" so i moved it to the desktop with 512MB RAM and 1Gig
swap where it segfaulted instead).
If anyone wants the complete code I can post a link.
In the mean time, does anyone have a suggestion of a more efficient
method within maxima?
I'm not that familiar with multidim newtons method, but it looks like
the method is
deltas: solve_the_matrix_equation_for_deltas(jacobian . deltas = guesses);
guesses: guesses + deltas;
my guess is that solve_the_matrix_equation_for_deltas can be written
much more efficiently than the method in mnewton (which calculates the
symbolic inverse, pre-multiplies by the symbolic equations,
substitutes the guesses and tries to evaluate the whole thing to a
floating point vector.
Any suggestions?
--
Daniel Lakeland
dlakelan@street-artists.org
http://www.street-artists.org/~dlakelan