python scipy leastsq fit with complex numbers

Question

I have a data set of complex numbers, and I'd like to be able to find parameters that best fit the data. Can you fit data in complex numbers using leastsq as implemented by scipy in python?

For example, my code is something like this:

 import cmath
 from scipy.optimize import leastsq
 def residuals(p,y,x):
      L,Rs,R1,C=p
      denominator=1+(x**2)*(C**2)*(R1**2)
      sim=complex(Rs+R1/denominator,x*L-(R1**2)*x*C/denominator)
      return(y-sim)

 z=<read in data, store as complex number>
 x0=np.array[1, 2, 3, 4]
 res = leastsq(residuals,x0, args=(z,x))

However, residuals doesn't like working with my complex number, I get the error:

File "/tmp/tmp8_rHYR/___code___.py", line 63, in residuals
    sim=complex(Rs+R1/denominator,x*L-(R1**_sage_const_2 )*x*C/denominator)
  File "expression.pyx", line 1071, in sage.symbolic.expression.Expression.__complex__ (sage/symbolic/expression.cpp:7112)
TypeError: unable to simplify to complex approximation

I'm guessing that I need to work only with floats/doubles rather than complex numbers. In that case, how can I evaluate the real and complex parts separately and then lump them back together into a single error metric for residuals to return?

Answer 1

The least squares function in scipy wants a real residual returned because it is difficult to compare complex values (e.g. is 1+2j greater or less than 2+1j?). Remember the residual is essentially a measure of the quality of the set of parameters passed in, it tells leastsq how close to the true fit it is.

What you can do is add the error (y-sim) in quadrature, appending these lines after you calculate 'sim' in your residuals function:

a = y-sim
return a.real**2 + a.imag**2

So long as y and sim are both np.array's of complex's then this will work and is relatively efficient.