How to solve complex least squares with real unknowns

Question

I am looking at the case where my least squares problem is of the form

$ y=A x $

where y is a vector of measured complex values, A is the model or design matrix, also complex, and x are the model coefficients which are real. How do I formulate the problem in Mathematica to get real values for my model coefficients?

Here is a minimum working example of what I have looked at.

First I make a noise free set of artificial data where my unknows are the values 1 to 7 and use the ordinary LeastSquares program to find my unknowns.

mat = RandomComplex[{-10 (1 + I), 10 (1 + I)}, {20,7}]; (* Design matrix *)
uk = Range[7]; (* Known unknowns *)
yy1 = mat . uk; (* Measured values *)
sol1 = LeastSquares[mat, yy1] (* Solution *)


(* {1. + 1.93795*10^-15 I, 2. - 2.22045*10^-16 I, 3. - 5.3385*10^-16 I, 
 4. - 2.02397*10^-15 I, 5. + 4.70461*10^-16 I, 6. - 5.55476*10^-16 I, 
 7. + 5.52099*10^-17 I} *)

As expected the real part of the solution is perfect and the imaginary part is numerical noise.

Now I put some noise into the measured values and use LeastSquares again.

yy2 = yy1 +  RandomComplex[{-1 - I, 1 + I}, 20];(* Measured values with noise *)
sol2 = LeastSquares[mat, yy2] (* Solution *)

(* {1.02302 + 0.0244589 I, 2.00445 - 0.051833 I, 2.97412 + 0.0227336 I, 
 4.0041 - 0.00438793 I, 5.01066 - 0.0132584 I, 5.96674 + 0.0293913 I, 
 7.00924 + 0.00691532 I} *)

As expected there are now complex solutions although the imaginary parts are small.

I now attempt to do the matrix algebra the long way to get a solution with real values. There is probably a better way of doing the matrix algebra (please make suggestions). I finish by comparing the solutions from the three methods.

  xx = Array[x, 7]; (* Unknown unknowns *)
    ee = Sum[
ComplexExpand[(mat[[n]] . xx - yy2[[n]]) Conjugate[mat[[n]] . xx - yy2[[n]]]],
         {n, 20}]; (* Sum of errors *)
    mat3 = D[ee, #] & /@ xx; (* Take derivaties *)
    {vec, mat4} = CoefficientArrays[mat3, xx];(* Find LHS and design matrix *)
    sol3 = -Inverse[mat4] . vec ;(* Solve *)
    TableForm[Transpose[{sol1, sol2, sol3}], 
     TableHeadings -> {None, {"Exact", "Complex Least Squares","Real Least Squares"}}]

So I have got my solution to have real values. However, I note that the errors in the real part of the LeastSquares solution are very similar to the errors in my real solution.

Two questions:

1. What is the best way to tackle this problem?
2. Is the LeastSquares solution just as good if I take the real part?

Note that my actual problem will have very large matrices.

Thanks

Answer 1

I am answering my own question. This turns out to be easier than I thought. Consider the problem where we have

$ y=A x $

with y a column of complex measurements, A a matrix of complex values and x a vector of real unknowns. The objective is to find the real values of x via a least squares process.

Split the complex values into real and imaginary parts

$ y_r+i y_i=x \left(A_r+i A_i\right) $

where the subscripts r, i mean real and imaginary parts. This may be written as two equations which may be put into matrix form as

$ \left( \begin{array}{c} y_r \\ y_i \\ \end{array} \right)=\left( \begin{array}{c} A_r \\ A_i \\ \end{array} \right) x $

With this configuration the LeastSquares function may be applied directly.

Here is an example.

nn = 100; (* Number of measurements*)
nk = 7; (* Number of unknowns *)
A = RandomComplex[{-10 (1 + I), 10 (1 + I)}, {nn, 
   nk}]; (* Design matrix *)
uk = Range[nk]; (* Known unknowns *)
y = A . uk + 
  RandomComplex[{-1 - I, 1 + I}, 
   nn]; (* Measured values with complex noise *)
Ar = Re[A];
Ai = Im[A];
yr = Re[y];
yi = Im[y];
yy = Join[yr, yi];
AA = Join[Ar, Ai];
LeastSquares[AA, yy]

(* {0.997568, 1.99928, 3.00159, 3.99664, 5.00797, 5.98345, 7.00619} *)

So this works well.