# How to create internally optimized expression for computing with high WorkingPrecision?

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.

-
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. –  Ted Ersek May 18 '13 at 11:25
@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... –  Alexey Popkov May 18 '13 at 11:32
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. –  Ted Ersek May 18 '13 at 12:10
Very good question. A good answer will hopefully discuss the ExperimentalNumericalFunction, 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. –  Oleksandr R. May 19 '13 at 1:30
@Oleksandr RuleDelayed was a typo here, I use Rule both for "Residual" and "Jacobian". –  Alexey Popkov May 19 '13 at 11:20

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

Assuming that residualVect[vars] is the residual vector where vars is a list of variables of complete (long) model, and we have a short model where one of the vars is fixed to numerical value. This fixed variable along with its value we will store as a rule in the variable set. Let us define the residual vector and Jacobian as NumericalFunctions:

ClearAll[resV, resVF, jac];
resV = ExperimentalCreateNumericalFunction[vars, Sqrt[2] residualVect[vars],
{Length[residualVect[vars]]}, WorkingPrecision->20];
resVF /; VectorQ[vars /. set, NumberQ] := resV[vars /. set];
jac /; VectorQ[vars /. set, NumberQ] := Drop[resV["Jacobian"[vars /. set]], None,
First@Position[vars, set[[1]]]];


Note that we must multiply the true residual vector by Sqrt[2] in order to get result consistent with ordinary output of FindMinimum where the first argument would be residualVect[vars].residualVect[vars].

Now FindMinimum can be used as follows:

FindMinimum[Null, initGuess, Method -> {"LevenbergMarquardt",
"Residual" -> resVF,
"Jacobian" -> {jac, EvaluationMonitor :> ++steps}},
WorkingPrecision -> 20]


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

-