MATLAB 2016a and 2017a give different results when multiplying by the complex conjugate

I am running some code and found that multiplying a complex vector x by x' gives me slightly different results between MATLAB2016a (gives me answer I am looking for) and 2017a (has a complex component on the diagonal when I am not expecting it). If I troubleshoot with randn(3)+randn(3)*i I have no problems (with rng seeded). Is it due to precision? An example x vector that will give me trouble is below.

x=

4.07165372227704e-09 - 1.57176600601260e-08i

-1.39876117489891e-07 - 1.17688977471939e-07i

2.19114480796579e-07 + 1.33125665568044e-08i

3.68460108310824e-08 + 9.08472574818831e-08i

This gives me further trouble even if I run the result of Rxx=x*x' from 2016a on 2017a when I use the eig() function, eig(Rxx), so I am hoping it's a parameter I can change that will fix both results in MATLAB.

Thanks, CJ

Answer 1

[Christine Tobler]

This is wrong behavior: the diagonal of matrix Rxx should have zero imaginary part. We will work on fixing this for a future release. Thank you very much for letting us know about this.

Note the issue specifically happens when computing the outer product of a vector with itself. For matrices, the diagonal is still real:

 >> X = [x, zeros(4, 1)];
 >> isreal(diag(x*x'))
 ans =
   logical
    0
 >> isreal(diag(X*X'))
 ans =
   logical
    1

As a workaround, remove the imaginary part of the diagonal elements explicitly:

 Rxx = Rxx - 1i*diag(imag(diag(Rxx)));

This should also resolve the issue in EIG: The problem there would have been that EIG uses a different algorithm when the input matrix is hermitian, and ishermitian(A) requires the diagonal of A to be real.

[Christine Tobler]

A colleague is looking into this, and suggested a simpler workaround for R2017a:

 if iscolumn(x)
     Rxx = x.*x';  % Uses implicit expansion behavior
 else
     Rxx = x*x';
 end

[CJ Werner]

I used the if else method (didn't try the first one). That seemed to fix everything! Thanks.