Matlab中将复矩阵赋值给局部变量时出现异常精度丢失问题

I ran into this exact edge case back in the MATLAB 2018a Linux days—super frustrating when algebraic guarantees break due to compiler quirks. Let's break down what's happening, why it occurs, and how to work around it:

Problem Overview #

When computing transpose(a)*conj(a) (where a is a complex matrix) inside a function, the resulting matrix loses its hermitian property: diagonal elements end up with non-zero imaginary parts, and ishermitian() returns false. But run the exact same line in a script or the REPL, and it behaves as expected (returning a proper hermitian matrix). The equivalent expression conj(a'*a) works correctly in both contexts, and real matrices never exhibit this issue.

Minimal Reproducible Example #

Here's the code to confirm the bug (as you provided):

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) % Output: 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' inside function, but returns 'false' if checked in debugger REPL
end

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

Root Cause #

The key difference is JIT compilation: MATLAB automatically JIT-compiles functions to speed up execution, but scripts are typically interpreted (or compiled with looser optimizations). The JIT compiler in 2018a Linux seems to reorder or optimize floating-point operations in transpose(a)*conj(a) in a way that breaks the algebraic symmetry required for hermitianness.

Theoretically, transpose(a)*conj(a) must be hermitian: taking the conjugate transpose gives (A^T \overline{A})^H = \overline{(A^T \overline{A})^T} = \overline{\overline{A}^T A} = A^T \overline{A}. But the JIT's optimizations introduce asymmetric floating-point errors, leading to tiny non-zero imaginary parts on the diagonal (which ishermitian() picks up as a failure).

Workarounds & Fixes #

• Use the algebraically equivalent safe expression: Replace transpose(a)*conj(a) with conj(a'*a). Since a' is the conjugate transpose, a'*a is already hermitian, and taking its conjugate leaves it unchanged. This expression is handled correctly by the JIT compiler.
• Cast to single precision (if acceptable): As shown in your example, converting the result to single() restores the hermitian property—likely because the single-precision floating-point model has different error behavior that cancels out the asymmetry.
• Disable JIT for the problematic function (not recommended for performance): You can add feature('JIT', 'off'); at the start of the function, but this will slow down execution significantly.

Language Recommendation for Algebraic Safety #

If you're frustrated by issues where numerical code breaks algebraic guarantees, languages with first-class support for matrix algebraic structures are a great fit:

• Julia: Has built-in Hermitian and Symmetric types that the compiler uses to enforce structural properties during computations. This eliminates most cases where optimizations break algebraic symmetry.
• Python with SciPy: Offers scipy.linalg.Hermitian and similar classes that preserve structure through operations, avoiding these kinds of low-level compiler bugs.

内容的提问来源于stack exchange，提问作者oyd11