It appears that we are getting a little sidetracked here, or maybe getting off on the wrong foot.
The code package I offer at jaeckel.org is a minimalist C/C++ interface to what is (mainly) machine-translated Fortran code. Anyone is at liberty to use it any way they wish and have their own frontend to it.
I, for example, have used it with C# frontends, Python frontends, and even a frontend into an in-house developed scripting language (written in the language-language EBNF and compiler-compiled with ANTLR). This enabled the user to specify problems in C# or Python or that in-house language (where we used it, among other things, for commodity futures [and swap rates etc] time series mean reversion parameter calibration, e.g., maximum likelihood based, non-linear variations of calibration etc. etc.). I have used it for yield curve construction, parametric volatility surface calibration to market option prices, FX/credit risk correlation calibration, and myriad of other applications.
As for an 'evolved' C++ example, kindly have a look at the function fit_polynomial() in fit_polynomial.cpp. Its purpose is to demonstrate in Excel the result of a calibration of a polynomial to some scattered data. The function essentially does just this:
fit_polynomial(const XLOper & order, const XLOper & abscissas, const XLOper & ordinates, const XLOper & use_jacobian) {
PolynomialFit poly(abscissas.to_vector(), ordinates.to_vector(), use_jacobian | false);
std::vector<double> coeffs(order.to_int()+1);
NL2SOL nl2sol;
nl2sol.tryToSolve(poly, coeffs, poly.data_length());
...
The syntax 'use_jacobian | false' in this context means that the input XLOper 'use_jacobian' defaults to 'false' (if omitted in the worksheet function in Excel, the value 'false' is used). Anyone who does not like this syntax has no obligation to use it - do your own thing!
The invocation of NL2SOL is done here via the NL2SOL class member function (aka 'method')
double tryToSolve(Functor& functor, std::vector<double> &x, size_t n_f);
which, as you can see, uses a 'std::vector<double> &' data type, and no raw pointers, and the NL2SOL algorithm safely takes care of the function value output vector memory allocation (and clean-up) of size n_f. Obviously, the bespoke functor class PolynomialFit has to comply with the NL2SOL interface, and that is how I wrote it, but you or anyone else can write whatever code they wish. There is no need to criticize this generic interface here - it is provided '
as is', and it is provided as it is for very good reasons (compatibility, simplicity, speed, etc.). NL2SOL itself has other reasons for the raw pointer interface, e.g., it keeps more than one copy of the x[] vector and may invoke the functor or Jacobian on any of those copies as it sees fit. It keeps all memory in one monolithic block which in turn is a performance benefit, so no need to criticize that either.
Btw, thank you for sharing with us that you hardly ever use or need raw pointers these days, and rest assured none of use would have assumed that you do, nor, I presume, do any of us, unless, of course
there is absolutely no other way (without major compromises), such as in the context of interfacing between run time environments or when dealing with machine-translated code. Please forgive my tongue-in-cheek there... we must retain a sense of humour!
NL2SOL has also been used at hedge funds for identification of globally optimal investment choices (the details are obviously confidential).
With the provided interface, it is straightforward to adapt other language frontends, e.g., I can see immediately how one would expose this to Java (based on my experience with C# which is virtually identical to Java in this aspect). I bet that sooner or later someone will have interest in a frontend to Rust.
All such interfacing is necessarily always based on the common denominator between run time environments, and that is the raw data pointer / raw function pointer (which is literally a code execution entry pointer). None of this is affected by the evolution of any language, C++ or otherwise. Ultimately, it's machine-hardware/OS-native-sized integers, and pointers are of course just that.
I am offering my 25(+)-years-curated version of the code '
as is' as a courtesy to the community. The only value I add is the few algorithm/bug fixes that are explained in detail in my notes right beneath my NL2SOL.7z download link at jaeckel.org - please also have a look at what I wrote about the history, other usages of NL2SOL, etc. As for the meaning of NL2SOL's control parameters and return values, I refer to the original authors' documentation, available via their two mathematical articles and their detailed textual explanations in the comments of their source code (see nl2sol.cpp). Regarding the minimalist C++ frontend, e.g., the Jacobian definition, I am not sure how this description could be made any clearer (you could always look up 'column-major' on the web), but anyone is of course at liberty to do so in their own copy of the code:
// The Jacobian function may be left undefined in which case finite-differencing is used.
virtual bool have_jacobian() const { return false; }
// The return value FALSE indicates that x[] was not inside the admissible domain, i.e., that constraints are violated.
// NOTE:
// The Jacobian function must store the output matrix J in the output vector J[] in column-major order (ForTran).
// For example, if n_f is the length of f[], and n_x is the length of x[], then J(k,l) is to be stored in J[k+l*n_f].
// The element J(k,l) represents the derivative
// ∂ f[k]
// ------
// ∂ x[l] .
virtual bool jacobian(const double *x, size_t n_x, size_t n_f,
double *J /* J[k+l*n_f] = ∂f[k]/∂x[l]. The length of J[] must be at least n_f*n_x */){
return false;
}
I also give my advice from experience that NL2SOL is the best I have ever worked with for these many purposes mentioned above (and I know of other teams who have done exhaustive comparisons with other algorithms and came to the conclusion that NL2SOL was superior by a wide margin to all else they tested).
Besides all of the above advocating of NL2SOL, I must caution and/or concede that NL2SOL is
not a method that outright and only uses a gradually collected sequence of vectors to estimate the Hessian. When NL2SOL in its evaluation course considers the Jacobian to be the next best step to make progress, it will evaluate it, either by the user provided function, or by finite-differencing approximation (which is done as a sequence of gradient vector computations). BFGS is closer to the idea of "Hessian estimation", as mentioned by the originator of this thread, but I found BFGS to be much, much, much slower in all applications where a local gradient approximation is admissible, i.e., can be used as an alternative (not all underlying problems allow that: a model function that is based on a random forest, for example, is internally a piecewise constant function with many many different levels, meaning, the local gradient is always zero! A similar issue arises in certain Monte-Carlo method
training stages - see the chapter on Bermudan swaption pricing in my book for an example where I did use BFGS for the same reason).
So, after all this, I point out with remorse that my posting about NL2SOL is, alas, effectively,
off-the-original-topic.
Nevertheless, I am hopeful that it may be of interest and/or use to some.
I finish with the polite request to kindly honour my conclusion:
as always,
positive feedback is welcome.
pointers as function arguments are kind of out-of-fashion, they are a throwback to C. I always used (const) refferences since 1990.