21
$\begingroup$

I have large dataset and need to fit rather complicated function on it with different values of one of its parameters (this parameter must be fixed in every fit). I use the "LevenbergMarquardt" algorithm of FindMinimum as described here (with the only difference that my function is not a "black-box").

I notice that if I set 

Method -> {"LevenbergMarquardt", 
   "Residual" -> Sqrt[2] residualVector[optimVariables], 
   "Jacobian" -> {"Symbolic", EvaluationMonitor :> ++steps}}

then FindMinimum takes 5 minutes to create an optimized form of symbolic jacobian and after this every evaluation of the jacobian takes only 2 seconds.

I need to fit my function with many different values of the parameter and to spend 5 minutes for making the jacobian that is identical for all the fits (with exception for the value of only one parameter) is a waste of time. I tried to compute symbolic form of the jacobian by myself and got identical results with automatic symbolic jacobian. I used the code:

jacobianMatrix[_List?(VectorQ[#, NumberQ] &)] = 
    D[Sqrt[2] residualVector[optimVariables], {optimVariables}]

and 

Method -> {"LevenbergMarquardt", 
   "Residual" -> Sqrt[2] residualVector[optimVariables], 
   "Jacobian" -> {jacobianMatrix[optimVariables], EvaluationMonitor :> ++steps}}

The problem is that each evaluation of this symbolic jacobian with numerical values of parameters takes 4 minutes! Obviously FindMinimum optimizes its internal representation in some way and it gives huge speedup. But FindMinimum creates the jacobian for every fit again. Is it possible to optimize the internal representation manually?

I should stress that I compute with WorkingPrecision higher than MachinePrecision. So probably compilation is not an option.

|improve this question
$\endgroup$
  • 1
    $\begingroup$ It might work better if you make jacobianMatrix a CompiledFunction. That will help if FindMinimum uses MachinePrecision to get close to the optimal solution, and switch to higher precision for fine tuning. $\endgroup$ – Ted Ersek May 18 '13 at 11:25
  • $\begingroup$ @Ted Good idea, but I would prefer to optimize the jacobianMatrix once because high precision significantly decreases the number of steps needed to get the minimum: large data set results in large rounding-off errors in summation... $\endgroup$ – Alexey Popkov May 18 '13 at 11:32
  • 2
    $\begingroup$ performing arithmetic operations in the right way can reduce round off errors. One important example is use of HornerForm (a built in function) for computing polynomials. Most text books on numerical methods show a few other tricks for specific cases. It's all about avoiding catastrophic cancelation, which is the loss of precision due to subtracting numbers that are nearly equal but not close to zero. $\endgroup$ – Ted Ersek May 18 '13 at 12:10
  • 3
    $\begingroup$ Very good question. A good answer will hopefully discuss the Experimental`NumericalFunction, which is a structure produced by FindMinimum from the function and its Jacobian that is optimized for fast numerical evaluation. I know almost nothing about these objects or how to create/use them, but I would like to find out. By the way, I am not sure if it is just a typo, but it is going to be a problem that you used RuleDelayed for the "Jacobian" option in the second case. As it is currently, the 4 minute delay may be due to re-deriving the Jacobian at each point. $\endgroup$ – Oleksandr R. May 19 '13 at 1:30
10
$\begingroup$

On the base of findings in this thread the procedure is as follows.

Assuming that residualVect is the residual vector and initGuess is the initial guess for FindMinimum, let us create an internally optimized residual vector via Experimental`CreateNumericalFunction:

resV = Experimental`CreateNumericalFunction[initGuess[[;; , 1]], 
   Sqrt[2]*residualVect, {Length[residualVect]}, WorkingPrecision -> 20];

Note that we have multiplied the true residual vector by Sqrt[2] in order to get the result consistent with the ordinary FindMinimum syntax where the first argument is equal to Total[residualVect^2].

Now FindMinimum can be used as follows:

FindMinimum[Null, initGuess, Method -> {"LevenbergMarquardt", 
   "Residual" :> resV[initGuess[[;; , 1]]], 
   "Jacobian" :> resV["Jacobian"@initGuess[[;; , 1]]]},
                             WorkingPrecision -> 20]

This approach gives more than one order of magnitude speedup in my application.

|improve this answer
$\endgroup$
5
$\begingroup$

Thanks to the recent post by ilian who uncovered and explained the undocumented Optimization`FindFit`ObjectiveFunction, I have found Optimization`FindFit`ResidualFunction which seems to be ideally suited to the original goals of the question. With the latter we don't even need to construct the residual vector explicitly, we just specify the model function and the data and get the optimized NumericalFunction automatically! A little example from the linked thread:

data = BlockRandom[SeedRandom[12345];
   Table[{x, Exp[-2.3 x/(11 + .4 x + x^2)] + RandomReal[{-.5, .5}]}, {x, 
     RandomReal[{1, 15}, 20]}]];

residual = Optimization`FindFit`ResidualFunction[data, Exp[a x/(b + c x)], {a, b, c}, x]

Experimental`NumericalFunction[{a,b,c},Block[{x={2.69745,5.61891,11.9585,7.02236,4.1302,7.48274,11.3322,10.9067,12.0728,2.47999,6.32617,5.88977,8.63809,11.2304,11.0019,8.65308,4.26835,9.26918,3.04309,9.8615},Optimization`FindFit`y$588={0.972294,0.344088,0.767092,0.460704,0.947346,1.21116,0.439919,0.675085,1.17798,0.961116,0.720896,1.18664,0.605936,0.811138,0.861414,0.670646,0.770535,0.365288,0.649543,0.551157}},E^((a x)/(b+c x))-Optimization`FindFit`y$588],-NumericalFunctionData-]

screenshot

|improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.