Least Squares Minimization Complex Numbers

Question

I've been using my Matlab, but it's my vision to eventually switch over to doing all of my analysis in Python since it is an actual programming language and a few other reasons.

The recent problem I've been trying to tackle is to do least squares minimization of complex data. I'm an engineer and we deal with complex impedance pretty often, and I'm trying to use curve fitting to fit a simple circuit model to measured data.

The impedance equation is as follows:

Z(w) = 1/( 1/R + j*w*C ) + j*w*L

I'm then trying to find the values of R, C, and L such that the least squares curve is found.

I've tried using the optimization package such as optimize.curve_fit or optimize.leastsq, but they don't work with complex numbers.

I then tried making my residual function return the magnitude of the complex data, but this did not work either.

Answer 1

Referring to unutbu answer's, there is no need to reduce the available information by taking the magnitude squared in function residuals because leastsq does not care whether the numbers are real or complex, but only that they are are expressed as a 1D array, preserving the integrity of the functional relationship.

Here is a replacement residuals function:

def residuals(params, w, Z):
    R, C, L = params
    diff = model(w, R, C, L) - Z
    z1d = np.zeros(Z.size*2, dtype = np.float64)
    z1d[0:z1d.size:2] = diff.real
    z1d[1:z1d.size:2] = diff.imag
    return z1d

This is the only change that need be made. The answers for the seed 2013 are: [2.96564781, 1.99929516, 4.00106534].

The errors relative to to unutbu answer's are reduced by significantly more than an order of magnitude.

Answer 2

Ultimately, the goal is to reduce the absolute value of sum of the squared differences between the model and the observed Z:

abs(((model(w, *params) - Z)**2).sum())

My original answer suggested applying leastsq to a residuals function which returns a scalar representing the sum of the squares of the real and imaginary differences:

def residuals(params, w, Z):
    R, C, L = params
    diff = model(w, R, C, L) - Z
    return diff.real**2 + diff.imag**2

Mike Sulzer suggested using a residual function which returns a vector of floats.

Here is a comparison of the result using these residual functions:

from __future__ import print_function
import random
import numpy as np
import scipy.optimize as optimize
j = 1j

def model1(w, R, C, L):
    Z = 1.0/(1.0/R + j*w*C) + j*w*L
    return Z

def model2(w, R, C, L):
    Z = 1.0/(1.0/R + j*w*C) + j*w*L
    # make Z non-contiguous and of a different complex dtype
    Z = np.repeat(Z, 2)
    Z = Z[::2]
    Z = Z.astype(np.complex64)
    return Z

def make_data(R, C, L):
    N = 10000
    w = np.linspace(0.1, 2, N)
    Z = model(w, R, C, L) + 0.1*(np.random.random(N) + j*np.random.random(N))
    return w, Z

def residuals(params, w, Z):
    R, C, L = params
    diff = model(w, R, C, L) - Z
    return diff.real**2 + diff.imag**2

def MS_residuals(params, w, Z):
    """
    https://stackoverflow.com/a/20104454/190597 (Mike Sulzer)
    """
    R, C, L = params
    diff = model(w, R, C, L) - Z
    z1d = np.zeros(Z.size*2, dtype=np.float64)
    z1d[0:z1d.size:2] = diff.real
    z1d[1:z1d.size:2] = diff.imag
    return z1d

def alt_residuals(params, w, Z):
    R, C, L = params
    diff = model(w, R, C, L) - Z
    return diff.astype(np.complex128).view(np.float64)

def compare(*funcs):
    fmt = '{:15} | {:37} | {:17} | {:6}'
    header = fmt.format('name', 'params', 'sum(residuals**2)', 'ncalls')
    print('{}\n{}'.format(header, '-'*len(header)))
    fmt = '{:15} | {:37} | {:17.2f} | {:6}'
    for resfunc in funcs:
        # params, cov = optimize.leastsq(resfunc, p_guess, args=(w, Z))
        params, cov, infodict, mesg, ier = optimize.leastsq(
            resfunc, p_guess, args=(w, Z),
            full_output=True)
        ssr = abs(((model(w, *params) - Z)**2).sum())
        print(fmt.format(resfunc.__name__, params, ssr, infodict['nfev']))
    print(end='\n')

R, C, L = 3, 2, 4
p_guess = 1, 1, 1
seed = 2013

model = model1
np.random.seed(seed)
w, Z = make_data(R, C, L)
assert np.allclose(model1(w, R, C, L), model2(w, R, C, L))

print('Using model1')
compare(residuals, MS_residuals, alt_residuals)

model = model2
print('Using model2')
compare(residuals, MS_residuals, alt_residuals)

yields

Using model1
name            | params                                | sum(residuals**2) | ncalls
------------------------------------------------------------------------------------
residuals       | [ 2.86950167  1.94245378  4.04362841] |              9.41 |     89
MS_residuals    | [ 2.85311972  1.94525477  4.04363883] |              9.26 |     29
alt_residuals   | [ 2.85311972  1.94525477  4.04363883] |              9.26 |     29

Using model2
name            | params                                | sum(residuals**2) | ncalls
------------------------------------------------------------------------------------
residuals       | [ 2.86590332  1.9326829   4.0450271 ] |              7.81 |    483
MS_residuals    | [ 2.85422448  1.94853383  4.04333851] |              9.78 |    754
alt_residuals   | [ 2.85422448  1.94853383  4.04333851] |              9.78 |    754

So it appears which residual function to use may depend on the model function. I'm at a loss to explain the difference in result given the similarity of model1 and model2.