Scipy's leastsq with complex numbers

Question

I'm trying to use scipy.optimize.leastsq with complex numbers. I know there are some questions about this already but I still can't get my simple example working, which is complaining about casting from complex to real numbers.

If I did it right the solution to the below should be x=[1+1j,2j]:

import numpy as np
from scipy.optimize import leastsq

def cost_cpl(x,A,b):
    return np.abs(np.dot(A,x)-b)

A=np.array([[1,1],[2j,0]],dtype=np.complex128)
b=np.array([1+3j,-2+2j],dtype=np.complex128)
x,r=leastsq(cost_cpl,np.array([0+0j,0+0j]),args=(A,b))
print x
print r

But I'm getting

TypeError: Cannot cast array data from dtype('complex128') to dtype('float64') according to the rule 'safe'

EDIT: If I change the first guess from np.array([0+0j,0+0j]) to np.array([0,0]) the function runs but I get the wrong answer (a real one).

Answer 1

Since leastsq() can only accept real numbers, you need to use .view() method to convert between real array and complex array.

import numpy as np
from scipy.optimize import leastsq

def cost_cpl(x, A, b):
    return (np.dot(A, x.view(np.complex128)) - b).view(np.double)

A = np.array([[1,1],[2j,0]],dtype=np.complex128)
b = np.array([1+3j,-2+2j],dtype=np.complex128)
init = np.array([0.0, 0.0, 0.0, 0.0])
x, r = leastsq(cost_cpl, init, args=(A, b))
print(x.view(np.complex128))

output:

array([  1.00000000e+00+1.j,   4.96506831e-16+2.j])