首页
学习
活动
专区
圈层
工具
发布
社区首页 >问答首页 >(Matlab)当将复矩阵赋值给局部变量时,会导致精度损失。

(Matlab)当将复矩阵赋值给局部变量时,会导致精度损失。
EN

Stack Overflow用户
提问于 2019-02-19 15:26:43
回答 1查看 116关注 0票数 3

在将表达式存储到局部变量中时,我得到了这个意外的(似乎是非抱怨的)浮点精度损失,特别是复杂值的浮点精度损失,相同的语句在“脚本”作用域中正确运行。

(目前在Matlab-2018aLinux上进行测试,似乎是个bug,还没有在其他版本上进行测试)

计算是存储表达式transpose(a)*conj(a),其中a是一个复矩阵。存储等效的conj(a'*a)没有问题,但我想了解这个问题,我从同事的更大的代码库中提取了这个问题,得到了意想不到的结果。

在这个示例中,我使用ishermitian()查看是否发生了意外行为,因为这些表达式根据定义给出了Hermitian矩阵,而兼容的浮点舍入语义即使在失去精度时也将保持它为hermitian。如果矩阵是“真实”的,同样的事情不会发生。

我有以下最小可重现性示例:tst_bad()演示了不正确的版本。

代码语言:javascript
复制
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

有人看到过类似的行为吗?请注意,即使是病态m2矩阵的对角线成员-也不是真实的(如预期的)。

善后:

  1. @Cris Luengo的回复--“函数”代码在当前的Matlab引擎中主要是JIT'd,而'script‘代码--的重要提示通常不是,因此区别在于此。
  2. 因此,这不是任何变量作用域或属性(我已经开始尝试在“脚本-作用域”中分配一个变量并将其传递给函数,等等)
  3. 这让我很想用这样的语言进行数值编码,在这种语言中,计算结果将这样的假设显示到类型系统中(例如,这是Hermitian矩阵,等等,这取决于标量的代数结构),并进一步使用这些信息来选择算法等等。
EN

回答 1

Stack Overflow用户

回答已采纳

发布于 2019-02-19 17:10:23

您可以将测试简化为:

代码语言:javascript
复制
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

我刚刚在MATLAB R2018b (在线版本)上运行了这个程序,并确认了这个问题。在函数内或函数外部计算transpose(a)*conj(a)会导致不同的结果。这似乎是JIT的一个问题。

我建议你将其作为bug提交给MathWorks。( "report“链接是最右边的一个,就在蓝色栏下面,您需要一个有效的许可证才能以这种方式报告bug。)

票数 3
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接:

https://stackoverflow.com/questions/54769719

复制
相关文章

相似问题

领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档