Solve Linear Equations NaN return

[shoneill]

I've implemented a modified Linear Levenberg Marquardt function for gauss fitting of spectral data, but have problems from time to time where I simply can't fit a specific data set.

I've traced the problem to the function "Solve Linear Equations" outputting NaN for one or more solution vectors (Not all) which then plays havok with my function.

I've located a data set which is initially OK, but which after solving the linear equations has NaN as the second element of the solution vector array. Changing the matrix type to something other than "generic" gets rid of the NaN for this particular data set, but gives problems with others which work with the "generic" matrix type.

Since my mathematics lectures in college are a few years away now, I'm asking someone out there if they can help me with my problem. I can supply a "corrupting" data set if this would help, or maybe someone knows this problem and has a work-around?

I've thought of leaving the matrix type to "generic" and checking the output to see if there's a NaN present, and if so, then re-calculate with another matrix type until there's no NaN. This is most likely not very sound mathematically, but it might reduce the number of rogue data sets.

Thanks in advance

Shane.

[altenbach]

Shane,

I would be interested in looking into your data and VIs. How complex is your spectral data? How many gaussians?

I've never used "Levenberg Marquardt.vi" much (don't like it ;)). I always use vis based on "Nonlinear Lev-Mar fit.vi" and have written a few gaussian fitting routines (even 2D) that seem very stable. Do you still have my private e-mail address?

[shoneill]

Hi Altenbach,

The spectra can have up to 20 peaks.

I don't really like using the Levenberg Marquardt VI delivered with LabVIEW either, it's painfully slow. The version I've re-written is really a significant re-write throwing out all formula parsing and replacing it with coded Gauss (and other) fits.

What's the difference between linear and non-linear Lev-Mar functions (apart from the obvious). Is it possible to substitute a non-linear equation solve instead of the linear equation solve, or is the difference more complex than that. Do I need derivative equations for the non-linear version?

Of course I could look through the non-linear Lev-mar, but that'll probably take a while.

My code is painfully under-documented at the moment, and it's still a work in progress (trying to optimise a spectrometer calibration based on NIST spectral lines). As such, I don't really favour the option of sending my code. The basic structure is similar to the LV-own Levenberg Marquardt (just optimised).

Do you think there's a reason why the non-linear should not suffer from the same problem as the linear?

Thanks

Shane

[altenbach]

(Somebody correct me if I'm wrong!) I think both routines are basically the same, but written for a different input model (Formula vs. subVI). AFAIK, there is no "linear" Levenberg Marquardt algorithm.

Anyway, if you have multiple independent gaussian peaks, your system is ambiguous by nature, for example if you swap two peaks in the fitting function, you'll get exactly the same fit quality. You need good starting values for it to gravitate to a good solution. It seems some peaks might merge or wander off into the woods in your case. Put some graph probes deep inside your VI to see how the fitting function develops over time.

Forget about the derivative formula calculation for your model, it would be way too complicated. Numerical calculation works fine.

I was actually more interested in your data, not you code. Are all peaks well separated and of similar width or is there significant variation and overlap (e.g. a few narrow gaussians on top of a very wide gaussian).

[shoneill]

Hi Altenbach,

The data can have a significant amount of overlap. Typically three to four groups of peaks fully resolved, but each group can have up to 4-5 peaks in them, partially visible only as shoulders.

The start value "problem" is more or less solved (but as usual, not yet finished).

I've got an approximate first-guess peak find running which does a pretty good job of finding the peaks I need, and I then fit the height and FWHM first before doing a full optimisation (The Sub-VI approach allows much more flexible fitting models).

The fitting works well on almost all data sets. There are simply some spectra which (although visibly hardly different to others which work perfectly) do not fit, i.e. generate the NaN response from the "solve linear equations" function. It's a numeric problem (bug?), I'm sure, and not strictly a peak location or resolution problem. then again, I may be wrong.

I refer to the standard "LevMar" VI as being linear, because it assumes a linear relationship between the variable variation and the end mse used for optimisation. This is where the "Solve linear Equations" comes in. Since the relationship is almost certainly not linear (foe example when peaks overlap), I thought maybe the non-linear coefficient guess may yield better results. I've had a quick look through the non-linear LEV-MAR function, but don't understand it yet ot the extent I understand the linear one. It does indeed seem to take a slightly different approach (once you look past the whole parsing code of the "linear" function.

I'll need some time to get some understanding of the non-linear code.

Attached are some example spectra (One which works, and one which doesn't).

Thanks again

Shane

[DSPGuy]

Hi Shane, Altenbach,

Thought I might jump in and comment on the Nonlinear Lev-Mar.vi and Levenberg Marquardt.vi, and their differences. As Altenbach indicates, the "Levenberg Marquardt.vi" is just an alternate implementation of the same algorithm, but allowing the model to be expressed using formula strings. The formula string is parsed and stored as a parse-tree of operators, and constants. Its performance is due to the fact that every time the model is evaluated the parse tree must be evaluated. There is also a slight difference in the way the coefficients are updated, but this should not generally make any difference in the final results.

Both are implementations of the same algorithm, and this algorithm is a weighted-least-squares fit to a nonlinear model. After the algorithm has finished identifying the coefficients, there is a call to "Inverse Matrix.vi". This is to compute the covariance matrix. Please refer to Numerical Recipes in C for a good explanation and derivation of this detail, and the algorithm in general. The Lev-Mar algorithm is definitely sensitive to the initial guess coefficients. If these are too far off the algorithm can converge to a local minimum that is quite far from the global minimum.

I can't comment on your specific VI (modified Levenberg Marquardt) as I don't have a copy of it but in general you might get a NaN value from a linear system solution if the condition number of the matrix is high, perhaps indicating a degenerate rank condition. In other words, some of the rows or columns may be or nearly be linear combinations of each other. In this case, you might try solving the linear system using the pseudo-inverse. This might not be as accurate, but is usually much more robust.

You might also try and limit the data that you supply to the curve fit routine. For example, you mentioned that you can have 20 peaks in your data. Assuming they don't overlap too much, you might be able to remove the data that is far from a given peak, and let the fitting routine operate only on data that is near the peak.

[shoneill]

Thanks DSPGuy,

I've looked at the algorithm again, and the inverse matrix is called ony AFTER the final coefficients have been found. My VI has a problem DURING the coefficient search (loop) that it suddenly throws in a few NaN into the optimisation process (Which is generally unwelcome). As such, I think replacing the inverse matrix with the pseudo inverse matrix will not help.

I need to get out the old mathematics books and have a look at what's going on in the "solve Linear Equations" VI. Maybe I can code something else myself.....

Thanks.

Shane.

[shoneill]

OK,

I've just tried something Q&D. I've slapped a new function in after the linear equation VI which replaces any NaN with 0 (then to be added to the original value), basically removing any optimisation for that variable in this loop iteration. Then, in the next loop, the other non-NaN variables have changed slightly, and the variable which last time gave NaN is now correctly calculated. Go figure.

Quick and Dirty I know, but hey, it seems to work.

I'll test many more spectra next week and let you all know how it works.

Who needs maths when you've got dirty tricks up your sleeves 😉

Thanks guys

Shane.

[DSPGuy]

Shane,

Your approach of "ignoring" a result with NaNs and allowing the optimization to proceed is a well established and recognized approach to dealing with numerical difficulties of all kinds. If a cost function is discontinuous in places, sometimes the best way to deal with it is to simply ignore the strange results within a particular step and try to get the algorithm to step past the discontinuity.

When I mentioned using matrix pseudo-inverse it was in the context of solving the linear system, not the final inversion used to estimate the covariance matrix. A linear system solve can take the form Ax=b, where A and b are known, and we are solving for x. One way to solve this is to explicitly invert A (inv(A)), and premultiply each side by the inverse. The equation becomes Inv(A)Ax=Inv(A)b. Since Inv(A)A = I (identity matrix), the equation simplifies to x=Inv(A)b. So to solve a linear system, invert the matrix, and perform a matrix vector multiply with the inverted matrix and the known vector.

If the system is ill conditioned, the matrix inverse may yield NaNs or other exceptional values (Inf, -Inf). To avoid difficulties due to rank deficiency, the Inv(A) operation can be replaced by the pseudoInv(A). This yields a least squares solution to the deficient system.

I have coded a replacement for the linear system solve that specifically checks for the rank deficient error code from the Linear System Solve.vi, and also checks for the presence of NaN values in the solution vector. If either condition is True, then the VI resolves the system using the pseudoinverse. This might be better conditioned than simply keeping the previous solution. Then again, it might not. Anyway, hope this helps.

[shoneill]

DSPGuy,

wow, that sounds great. I'd love to try out the VI you posted, but I'm using LV 6.1.

Any chance of posting it in 6.1 form? Please?

I really need to get my old maths books out to understand what you just wrote. 😞

Thanks

Shane.