Constraining the least squares fitting in python

Question

I want to solve the following in a least-squares sense:

H = dot(A, B) + dot(A.conj(), C)

where the complex matrices H, B and C are known. The remaining complex matrix A (with its complex conjugate) is being searched.

I tried computing it with the (Python) numpy function:

x, res, r, singval = np.linalg.lstsq(np.vstack((B, C)), H)

However the results are not in a shape that I want( --> array((A, A.conj())).

How can I solve this?

Answer 1

The easiest way I found is separating the value of A into real and imaginary part:

A = U + 1j*V

therefore:

H = dot(U, B+C) + 1j*dot(V, B-C)

and with the use of scipy.optimize.lstsqr, the model is defined as:

def model(B, C, UV):
    U = UV[:len(UV)//2]
    V = UV[len(UV)//2:]
    H = np.dot(U, B+C) + 1j*np.dot(V, B+C)
    return H.view(float)

and the residuals as:

def residuals(params, B, C, H):
    UV = params
    diff = model(B, C, UV) - H.view(float)
    return diff.flatten()

and the result is obtained as follows:

params, cov = optimize.leastsq(residuals, x0 = np.ones(len(UV), dtype = float), 
                               args=(B, C, H))