Search code examples
matlabfloating-pointprecisionjitscoping

(Matlab) Strange precision loss, while assigning Complex-matrix to local variable

I'm getting this unexpected (and seems to be non-complaint) floating-point precision loss, specifically of Complex values, when storing an expression into a local variable, the same statement runs correctly in a 'script' scope

(currently testing on Matlab-2018a Linux, seems like a bug, haven't tested on other releases yet)

The calculation is storing the expression transpose(a)*conj(a) where a is a complex matrix. Storing the equivalent conj(a'*a) working without a problem, I want to understand the issue though, I've extracted this from a colleague's larger code-base that was getting unexpected results.

In the example, I'm using ishermitian() to see whether the unexpected behavior occurred, since these expressions by definition give Hermitian matrices, and compliant floating-point rounding semantics will keep it hermitian even when losing precision. The same thing doesn't happen if the matrix is "real" btw.

I have the following Minimal-Reproducible-Example: tst_bad() illustrates the ill-behaved version

len = 16;         % generate complex input:
t = 2*pi/len*(0:len-1)';
tt = t + (0:0.1:0.3);
a = hilbert(sin(tt));
% a = rand(len, 4)+1i*rand(len, 4); % alternative input, almost as good

m1 = transpose(a)*conj(a);
res = -ones(1, 4);
res(1) = ishermitian(m1);
[res(2) m2] = tst_bad(a);
[res(3) m3] = tst_good(a);
m4 = single(m1);
res(4) = ishermitian(m4)
display(res)
% 1 0 1 1

function [f, m] = tst_bad(a)
    m = transpose(a)*conj(a); 
    f = ishermitian(m);
    m_iseq = isequal(m,  transpose(a)*conj(a)) % 'true' but :
  % running "isequal(m,  transpose(a)*conj(a))" in the REPL will return 
  % 'false' if running in matlab-debugger
end

function [f, m] = tst_good(a)
    m = conj(a'*a);
    f = ishermitian(m);
end

Has anyone seen any similar behavior? Note that even the diagonal members of the ill m2 matrix - are not real (as expected)

Aftermath:

  1. Important note from @Cris Luengo's reply - 'function' code is mostly JIT'd in current Matlab engines, while 'script' code - is generally not, hence the difference.

  2. Thus, this isn't any variable scoping or property (I've started trying allocating a variable in the 'script-scope' and passing it to the function, etc)

  3. Pretty much makes me want to do numerical-code like that in a language where results of calculations manifest such assumptions into the type-system (such as stating - this is an Hermitian matrix, etc depending on the algebraic structure of scalars), and use this information further to choose algorithms etc.

Solution

  • You can simplify your test to:

    len = 16;
t = 2*pi/len*(0:len-1)';
tt = t + (0:0.1:0.3);
a = hilbert(sin(tt));

m1 = transpose(a)*conj(a);
disp(ishermitian(m1))      % true
m2 = tst_bad(a);
disp(ishermitian(m2))      % false
disp(isequal(m1,m2))       % false

function m = tst_bad(a)
    m = transpose(a)*conj(a);
end

    I've just run this on MATLAB R2018b (online version) and confirm the issue. Computing transpose(a)*conj(a) within a function or outside a function leads to different results. This looks like an issue with the JIT.

    I would suggest you submit this as a bug to the MathWorks. (The "report a bug" link is the rightmost one just underneath the blue bar, you need a valid license to report a bug this way.)