Segfault 5.9.1 on large matrix operations



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