Re: linalg[leastsqrs] in Maple V R4--success!!!

lcw1964 wrote:

The Svd path is the one I must take!

I happy to report success!

My original post was inspired by a desire to adapt the work at
http://www.library.cornell.edu/nr/bookcpdf/c5-13.pdf for my own use in
Maple. In a nutshell, Press et al. provide a routine that uses a
weighted least squares approach to give fairly good approximations to
minimax rational fits where convergence to the positively "best" equal
ripple minimax fit via the Remes algorithm is not required. Indeed, I
was looking for an alternative, since numapprox[minimax], which uses
numapprox[remez], was not finding a fit in a problem where I knew for
sure one should be found.

Apart from some tidying up in the user interface that is still needed,
I have managed to emulate the algorithm in Maple, taking advantage of
Maple's built-in linear algebra capabilities. Due to wisdom obtained
from others in this thread, I have relied on singular value
decomposition rather than the built in leastsqrs routine, since the
matrices involved can be ill-conditioned. I have been able to replicate
the example given in the book, and have been able to generate
approximations of the error function that do not perfectly mimic the
near best minimax ones of Cody--see
http://www.netlib.org/specfun/erf--but come gratifyingly close, so I
know I am on the right track. I have had to modify the routine to
accomodate whether I am minimizing absolute or relative error, but
fortunately this did not prove difficult.

The routine can be sluggish for larger problems where more than 10
parameters are being estimated, and in situations where a rational
approximation is particularly good to start with I have to set Digits
pretty high to mitigate rounding error in the final result, and this
slows things down too. Fortunately, I have found that the advice of
Press et al. holds true--there seems to be little benefit in doing more
than five iterations, and the best choice usually pops up in the second
or third iteration. The algorithm does not beat down all of the error
extrema into pristine equal ripple submission, but it does come
reasonably close and in some smaller cases the result looks equal
ripple to the naked eye, to my gratified amazement.

The bottom line is thanks to assistance in this NG I have, as a
complete rank amateur, enjoyed a little bit of Maple success, and I
thought an expression of pleasure and gratitude was in order.

Les

.

Relevant Pages

  • Re: Optimization help please
    ... few people that can improve such a routine using assembly. ... Have the compiler output assembler code for that routine and see if any ... be read then the CRC can be computed on each 'block' and the algorithm could ... > checksum is cpu and disk intensive, ...
    (comp.lang.asm.x86)
  • Re: Use of random bits
    ... The application I have in mind is the Algorithm P in Knuth's book ... items the first time through the loop, n-1 the second, etc., you could ... routine above: ... Notice how each time through the while loop, I had to use more decimal ...
    (sci.stat.math)
  • Re: my assembler is better than your assembler
    ... don't look at an existing algorithm and then turn ... My 'incredible slow' loop does maximal 31 iterations (2..3 cycles) ... Now you might argue that "this routine is so ... Even with that stack frame my code was considerably *faster* ...
    (alt.lang.asm)
  • Re: Division by rotation on 8 bit machine
    ... but I don't have a routine for the 8051 family ... Meaning, I assume: ... >I want algorithm for division by rotation which gets executed faster. ...
    (comp.arch.embedded)