Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I have a black-box multiargument multiparametric function of the type SRD[dataPoint_List,params_List] which accepts experimental data along with the parameters of the model and returns a vector with the first element being the squared relative deviation (SRD) between the computed model value and experimental value for the given data point and the rest of the elements being the gradient of the SRD in the given point with respect to params. I wish to minimize the sum Sum[First[SRD[dataPoint, params]], {dataPoint, dataPoints}] with respect to params using the Levenberg-Marquardt algorithm and the computed gradient of SRD with respect to params (equal to Sum[Rest[SRD[dataPoint, params]], {dataPoint, dataPoints}]). What is the best way to do this in Mathematica?

share|improve this question
2  
I'll only note that Method -> "LevenbergMarquardt" can take "Jacobian" and "Residual" options, which you should look into. –  J. M. Feb 18 '12 at 8:33
add comment

2 Answers 2

I had a similar situation once, here's a prototype similar to your problem that I wrote once on the Mathematica forum of LinkedIn, you can adapt it to your needs, the gradient is computed using the "Residual" option of the Method option.

targets={1,2};

f[{x_?NumericQ,y_?NumericQ}]:={(x-3)^5,2 (y-5)^3};
objectiveFunction[{params___?NumericQ}]:=Norm[f[{params}]-targets,2];
R[{params___?NumericQ}]:=f[{params}]-targets;

OptimVariables={Symbol["xx"],Symbol["yy"]};
initialGuess={1,-1};

FindMinimum[
   objectiveFunction[OptimVariables]
   ,
   {OptimVariables,initialGuess}\[Transpose]
   ,
   Method->{"LevenbergMarquardt","Residual"->R[OptimVariables]}
] 
share|improve this answer
add comment
up vote 2 down vote accepted

Here is the solution based on Faysal's code:

steps=0; 
FindMinimum[Null,
 {optimVariables, initialGuess}\[Transpose],
 Method -> {"LevenbergMarquardt", 
   "Residual" -> Sqrt[2] residualVector[optimVariables], 
   "Jacobian" -> {Sqrt[2] jacobianMatrix[optimVariables], EvaluationMonitor :> ++steps}, 
   "StepControl" -> {"TrustRegion", "StartingScaledStepSize" -> 1/1000, "MaxScaledStepSize" -> 1/10, "AcceptableStepRatio" -> 1/3}}]

Note that it is recommended to use exact numbers for the parameters of the "TrustRegion" method because these parameters are used inside of the algorithm without any check for consistency with WorkingPrecision. I should also note that the actual residual vector and jacobian must be multiplied by Sqrt[2] for having FindMinimum returning the minimum equal to

residualVector[optimVariables].residualVector[optimVariables]

and not to

residualVector[optimVariables].residualVector[optimVariables]/2

as it is by default.

The jacobian may be calculated automatically by the following code:

jacobianMatrix[_] = D[residualVector[optimVariables], {optimVariables}]

One can restrict the jacobian to be evaluated for numerical values only by defining it as:

jacobianMatrix[_List?(VectorQ[#, NumberQ] &)] = 
    D[residualVector[optimVariables], {optimVariables}]
share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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