0
$\begingroup$

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"}}]

enter image description here

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

$\endgroup$
4
  • $\begingroup$ PseudoInverse[mat].yy1, PseudoInverse[mat].yy2 and -PseudoInverse[mat4] . vec are alternative ways to solve $\endgroup$ Nov 25, 2021 at 7:22
  • $\begingroup$ @UlrichNeumann Thanks but they give the same results as my methods... $\endgroup$
    – Hugh
    Nov 25, 2021 at 9:32
  • $\begingroup$ It would be a bad alternative if not! $\endgroup$ Nov 25, 2021 at 9:56
  • $\begingroup$ @UlrichNeumann Do you have any idea about the internals of these codes? How valid is taking the real part of the solutions? $\endgroup$
    – Hugh
    Nov 25, 2021 at 10:38

1 Answer 1

2
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ This method assumes that the errors are independent (both between and within observations) and with a common variance. Do the residuals back up those assumptions? $\endgroup$
    – JimB
    Jun 4, 2022 at 14:55
  • $\begingroup$ @JimB You make a good point about the statistical assumptions. The concern is with data that are complex. If we assume that the complex data is derived via a Fourier transform then I would suggest the the real and imaginary parts are independent. Also, if the real time domain data has independent errors between observations then I think the complex data generated by Fourier transform would also be independent. Such independence would allow for standard least squares methods. However, do you have any insights into this area or can you suggest tests? $\endgroup$
    – Hugh
    Jun 5, 2022 at 19:43
  • $\begingroup$ The answers here stats.stackexchange.com/questions/66088/… might be helpful. $\endgroup$
    – JimB
    Jun 5, 2022 at 21:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.