numerics - Is it possible to use the LevenbergMarquardt algorithm for fitting a black-box residual function?

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?

Answer

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}]