I've difficulites with the NonlinearModelFit function.

In principle Mathematica should be able to deal with complex data.

E.g. if I define the following table

set = Table[{i, (3*i^2*I + 1*I)}, {i, 0, 10}]

and perform NonlinearModelFit, Mathematica has no problem finding the right values

NonlinearModelFit[set, (a + b*x^2*I), {a, b}, x]
(* {a -> 0. + 1. I, b -> 3. + 0. I} *)

But I need to fit some fractional functions, so I tried to fit just the inverse of the above set of data (I don't change any syntax, besides putting 1/ in front)

set = Table[{i, 1/(3*i^2*I + 1*I)}, {i, 0, 10}]
NonlinearModelFit[set, 1/(a + b*x^2*I), {a, b}, x]

Then Mathematica comes up with an error:

FindFit::nrlnum: The function value {1.-1. I,0.5-4.5 I,<<7>>,0.000152393-244.012 I,
<<1>>} is not a list of real numbers with dimensions {11} at {a,b} = {1.,1.}. >>

So now I'm confused- there is no syntax or misspelling problem – is Mathematica just not able to perform such a fractional fit? Or are there any mathematical restrictions, which I have not considered? I'm using Mathematica 7.


3 Answers 3


This is, at least in principle, a duplicate of jVincent's answer and the one I gave here. The general approach has been suggested by various people over the years, although I first encountered it here courtesy of Daniel Lichtblau the first time I needed to fit several datasets simultaneously.

I've been meaning to post this package for a while, ideally after having generalized it further, but given that generalization is somewhat complicated, and yet this remains a common question, I think on balance it's probably worthwhile to post the code as it stands. Despite certain limitations (listed below), it seems good enough for most applications that require fitting complex data.


SetAttributes[TransformedParameter, HoldRest];

Unprotect[TransformedFit]; ClearAll[TransformedFit];

Unprotect[ComplexFit]; ClearAll[ComplexFit];


(* Transform numeric quantities rather than renaming them *)
TransformedParameter[t_, num_?NumericQ] := t[num];

(* Generate unique symbols for each transformed parameter -- this avoids difficulties
caused by overzealous common subexpression elimination when the models are compiled *)
TransformedParameter[t_, p_] := TransformedParameter[t, p] =
   With[{sym = Unique["TransformedParameter$", Temporary]},
    (* Unset memo when cleared to facilitate garbage collection *)
    sym /: clear : (Clear | ClearAll | Remove)[___, sym, ___] :=
     clear /; (TransformedParameter[t, p] =.; True);

    (* Display as the parameter by itself if the transformation is Identity *)
    sym /: MakeBoxes[sym, form_] :=
     With[{boxes = MakeBoxes[p, form]},
       InterpretationBox[boxes, sym]
      ] /; t === Identity;

    sym /: MakeBoxes[sym, form_] :=
     With[{boxes = MakeBoxes[t[p], form]},
      InterpretationBox[boxes, sym]


$FitFunctions = If[$VersionNumber >= 7,
   {FindFit, NonlinearModelFit},

Options[TransformedFit] = {
   "FitFunction" -> First[$FitFunctions],
   "Transformation" -> Identity,
   "ParameterTransformation" -> Identity,
   "Hold" -> False

TransformedFit::cons = 
  "The constraint(s), `1`, should be given in terms of the transformed parameters only.";

    data_, {model_, cons_}, pars_, vars_,
   ] /; Internal`DependsOnQ[cons, pars] :=
  Message[TransformedFit::cons, cons];

(* Deal with data given as ordinate values only *)
   data_?VectorQ, model_, pars_, vars_,
   opts : OptionsPattern[{TransformedFit, Sequence @@ $FitFunctions}]
  ] :=
   Transpose[{Range@Length[data], data}], model, pars, vars,

   data_?MatrixQ, {model_, cons_} | {model_} | model_, pars_, vars_,
   opts : OptionsPattern[{TransformedFit, Sequence @@ $FitFunctions}]
  ] :=
    fitFunction = If[MemberQ[$FitFunctions, #], #, First[$FitFunctions]] & @ OptionValue["FitFunction"],
    transformations = List@OptionValue["Transformation"]~Flatten~1
     transformedParameters, unusedParameterMask, parameterRules,
      abscissae = Take[data, All, {1, -2}],
      ordinates = Take[data, All, {-1}]
     transformedData = {
         ConstantArray[Range@Length[transformations], {Length[abscissae], 1}]~Transpose~{2, 3, 1},
         ConstantArray[abscissae, Length[transformations]],
        }~Flatten~{{2, 3}, {1, 4}};

    transformedParameters = Outer[
       transformations, {pars}~Flatten~1
      ] // Transpose;
      (* Original and transformed parameters without initial guesses *)
      originalParameterNames = Replace[pars, {p_, __?NumericQ} :> p, {1}],
      transformedParameterNames = Replace[transformedParameters, {p_, __?NumericQ} :> p, {2}]
        (* Representations of the original parameters in terms of their transformations *)
        parameterRepresentations = OptionValue["ParameterTransformation"] @@@ transformedParameterNames
       unusedParameterMask = MapThread[
         Composition[Thread, Unevaluated, FreeQ],
         {parameterRepresentations, transformedParameterNames}
       Clear @@  Flatten@Pick[transformedParameterNames, unusedParameterMask];
       parameterRules = Thread[originalParameterNames -> parameterRepresentations];

      reparameterizedModel = model /. parameterRules,
      KroneckerDelta = If[Equal[##], 1, 0] & (* compilable *)
     transformedModel = Inner[
        #1[reparameterizedModel] KroneckerDelta[i, #2] &,
        transformations, Range@Length[transformations]

    (* PERFORM FIT *)
    If[TrueQ@OptionValue["Hold"], Composition[Hold, fitFunction], fitFunction][
     {transformedModel, cons},
     Pick[transformedParameters, unusedParameterMask, False]~Flatten~1,
     {i, vars}~Flatten~1,
     FilterRules[{opts}, Options[fitFunction]]


coordinateSystemRules["Cartesian"] = Sequence[
   "Transformation" -> {Re, Im},
   "ParameterTransformation" -> (#1 + I #2 &)
coordinateSystemRules["Polar"] = Sequence[
   "Transformation" -> {Abs, Arg},
   "ParameterTransformation" -> (#1 Exp[I #2] &)
coordinateSystemRules["Real"] = Sequence[
   "Transformation" -> {Re, Im},
   "ParameterTransformation" -> (#1 &)
coordinateSystemRules["Imaginary"] = Sequence[
   "Transformation" -> {Re, Im},
   "ParameterTransformation" -> (I #2 &)
(* Default to Cartesian coordinates *)
coordinateSystemRules[_] = coordinateSystemRules["Cartesian"];

Options[ComplexFit] = {
   "CoordinateSystem" -> Automatic

   data_, model_, pars_, vars_,
   opts : OptionsPattern[{ComplexFit, TransformedFit, Sequence @@ $FitFunctions}]
  ] :=
   data, model, pars, vars,
   FilterRules[{opts}, Except["CoordinateSystem" | "Transformation" | "ParameterTransformation"]]




Package (.m) and notebook files are also available.

The primary limitations are:

  • the Weights option is not (directly) supported, because it isn't clear to me how one should transform the weights in general when splitting a complex-valued function into a multivalued real mapping
  • the returned FittedModel objects still contain a reference to an index, i, that labels the coordinates (e.g. real line/imaginary line, modulus/argument, or whatever applies to any other coordinate system one may choose), because the structure of these objects is undocumented and I didn't yet figure out how to split them up
  • the transformation is done quite rigidly and is not currently versatile enough to cater for all foreseeable situations

Anyway, let's give it a try:

  Table[{i, I + 3*i^2 I}, {i, 0, 10}],
  a + b*x^2 I, {a, b}, x,
  "FitFunction" -> NonlinearModelFit

Fit to fractional complex data in Cartesian coordinates

Or, in polar coordinates:

  Table[{i, I + 3*i^2 I}, {i, 0, 10}],
  a + b*x^2 I, {a, b}, x,
  "FitFunction" -> NonlinearModelFit,
  "CoordinateSystem" -> "Polar"

Fit to fractional complex data in polar coordinates

And, just for fun, here's an example using FindFit instead of NonlinearModelFit, and where the parameters and the initial guesses of their values are explicitly complex:

 Table[{x, (17.381 + 53.249 I) x^(1.897 + 0.632 I)}, {x, -10, 10}],
 (a x^b), {{a, 20 + 50 I}, {b, 2 + 0.5 I}}, x
(* -> {Re[a] -> 17.381, Im[a] -> 53.249, Re[b] -> 1.897, Im[b] -> 0.632} *)

This is also useful for fitting real-valued data where the model may become erroneously complex-valued for certain values of the parameters. For example, from the other question:

 {{0.0, 100.0}, {0.02, 81.87}, {0.04, 67.03},
  {0.06, 54.88}, {0.08, 44.93}, {0.10, 36.76}},
 a b^t, {a, b}, t,
 "CoordinateSystem" -> "Real"
(* -> {Re[a] -> 100.004, Re[b] -> 0.0000452493} *)
  • $\begingroup$ You could post this here or maybe here for more visibility. $\endgroup$
    – Szabolcs
    Mar 2, 2014 at 19:15
  • $\begingroup$ @Szabolcs thanks. I will post it in István's thread, which somehow I had missed until now. $\endgroup$ Mar 2, 2014 at 19:18
  • $\begingroup$ (+1) One additional limitation: in the most practical cases the conservative assumption is that the number of observations is equal to the number of complex data points, not twice the number. In such cases current implementation of ComplexFit significantly underestimates the confidence intervals and standard errors. $\endgroup$ Jul 15, 2014 at 7:09
  • $\begingroup$ @AlexeyPopkov very good point, and one that I had overlooked; thank you. I would obviously like to address this, but (a) I'm not sure how to without somehow hijacking functions that work on FittedModel objects and (b) the situation is arguably problem-dependent. I do agree with you though that a conservative estimate would be better. $\endgroup$ Jul 15, 2014 at 8:37
  • 1
    $\begingroup$ @AlexeyPopkov yes, that's possible at present (although these two sublists would simply be joined into one) but one cannot guess the weights automatically, so I simply didn't do anything with them. There is nothing stopping someone from using the Weights and VarianceEstimatorFunction options with ComplexFit provided they know how the transformation is done. For this you can specify "Hold" -> True and the transformed FindFit/NonlinearModelFit call is returned held rather than being evaluated. $\endgroup$ Jul 15, 2014 at 9:00

NonlinearModelFit seeks to reduce the sum of the squares of the residuals. In this case, because you have used complex numbers, these residuals are complex values. The residuals need to be real for NonlinearModelFit to work.

I haven´t encountered a situation with complex residuals before. To be quite honest, I am bit confused by the idea. The complex numbers are not an order-able field, but I guess this doesn't really matter...

At any rate, what you probably want to do is to use FindFit instead of NonlinearModelFit. It has an option that allows you to set how the residuals are compared to each other:

FindFit[set, 1/(a + b*x^2*I), {a, b}, x, 
 NormFunction -> (Norm[Abs[#], 2] &)]

If this isn't what you want, you can probably build a fitting function really easily by using NMinimize. You just need to explicitly define how you want to minimize the residuals.

  • 2
    $\begingroup$ I am just confused that NonlinearModelFit works in the first case with (a + b*x^2*I) , but not if I want to fit the inverse- which should be the same in principle. I know that there are different ways to solve this fitting problem - but i would like to use nonlinearmodelfit since it seems to be more effective for more difficult problems (which I have to deal with) and gives you some errors for the parameters $\endgroup$
    – Justus
    Mar 10, 2012 at 18:44

I can't say I don't know why NonlinearModelFit doesn't just work smoothly. As Searke mentions, having complex residuals is quite exotic, however typically you would simply use the norm of the difference. This however does not work still. So what I typically do is to take my complex data {{x[i],yR[i]+I yI[i]}...} and change it into {{1,x[1],yR[1]},{1,x[1],yR[1]}...,{2,x[1],yI[1]},{2,x[1],yI[1]}}. This way the first coordinate indicates whether you have the real or imaginary component, so you can just change your model to depend on one more parameter that switches between these. The first part can be put into a function as:


To select the domain the following function is used, since it avoids the problems with Part printing errors when it's not yet evaluated and also rounding the domain variable:

domainSwitch[p_?NumericQ,domain_] := domain[[Round@p]]

And then your fit is just:

set = Table[{i, 1/(3*i^2*I + 1*I)}, {i, 0, 10}];
res=NonlinearModelFit[domainExtend[set, {Re, Im}], 
     domainSwitch[domain, {Re, Im}]@(1/(a + b*x^2*I)), {a, b}, {domain, x}]

So this now allows you to get all the typical functions of the fittedModel from NonlinearModelFit like "ANOVATable". The downside is that the data is suddenly on a different format, and so is the returned function, so you need to take this into account when plotting the results:

 Plot[{res["Function"][1, x], res["Function"][2, x]}, {x, 0, 10}],
           Cases[res["Data"], {1, _, _}][[1 ;;, 2 ;; 3]],
           Cases[res["Data"], {2, _, _}][[1 ;;, 2 ;; 3]]

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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