levenberg-marquardt non-linear fit implemented in Lua + FFI

[Francesco Abbate <francesco.bbt@...>]

Hi all,

I'm glad to announce that I've successfully implemented the
levenberg-marquardt algorithm for non-linear least squares fit in
LuaJIT2 using the FFI calls. Everything is already in the gsl shell
git repository in the luajit2 branch, the implementation file is
'num/lmfit.lua.in'. The algorithm is already tested and debugged using
the NIST datasets ENSO, Thurber and Hahn1 (
http://www.itl.nist.gov/div898/strd/nls/nls_main.shtml ). The test can
be launched with dofile('benchmarks/lmfit/nist_test.lua') and you will
get also some nice plots as a bonus ;-)


The algorithm is just a translation of the GSL algorthm in Lua using
the FFI calls to the same GSL functions. The FFI was quite important
to have access to all the functions of linear algebra provided in GSL.
The key linear algebra functions used by the algorithm are the
function related to QR decomposition used to perform a "smart"
inversion of the matrix (J' J + \lambda D). The GSL algorithm itself
appears to be itself an adaptation of the fortran algorithm provided
in the minpack.

The accuracy of the Lua+FFI is perfect in all the tests that I've
made, so the results are really satisfying. I didn't do any benchmark
for the execution speed because even for the most complex tests the
execution time was negligible. I may possibly perform some tests in
the future but this is not a priority.

Otherwise, I've also implemented in GSL shell the ODE integrator
Runge-Kutta-Fehlberg of order 4(5) and the Runge-Kutta Prince-Dormand
of order 8(9). Both implementation outperform C in the case of system
of small dimensions. For the moment I don't have a working
implementation in array form but I'm planning to do it very soon.

I've also implemented the QAG quadrature algorithm for function
integration. This latter in an adaptive algorithm based on the
Gauss-Kronrod quadrature rules. I've also implemented the non-adaptive
Gauss-Kronrod integratr, the qng algorithm. The performance are on the
par with C but in this case C wins with a small marge.

For the curious I've implemented also a code to calculate on the fly
the gauss-kronrod coefficients, the code in the kronrod branch of the
git repository.

For the moment I'm preparing a new release of GSL shell where I will
introduce all the JIT related innovations. I would greatly appreciate
any help here, there are a lot of things to do for the release. I fear
that as usual I will be alone to do all the hard work, but that's it,
I'm getting used to that... :-)

I was also planning to write an article on that but I would like to
get the advice of the most experienced people here...

Best regards,
--
Francesco