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:
- What is the best way to tackle this problem?
- Is the
LeastSquares
solution just as good if I take the real part?
Note that my actual problem will have very large matrices.
Thanks
PseudoInverse[mat].yy1
,PseudoInverse[mat].yy2
and-PseudoInverse[mat4] . vec
are alternative ways to solve $\endgroup$