矩阵相乘得到了两个不同的答案

Question

这里有一些简单的代码，展示了我所看到的内容：

A = randn(1,5e6)+1i*randn(1,5e6);
B = randn(1,5e6)+1i*randn(1,5e6);

sum(A.*conj(B)) - A*B'
sum(A.*conj(B)) - mtimes(A,B')
A*B' - mtimes(A,B')

现在，底部显示的三种方法应该都是做同样的事情，所以答案应该是零，对吗？不对！差异很小，但不足以被认为是可以忽略的。此外，随着A和B的长度增加，误差也会增加。

有人知道这些方法之间的实际区别吗？我知道代码中可能有一些捷径，但如果可能的话，我想量化它们。 Matlab是否在任何地方发布了这些差异？我四处寻找，但没有找到任何信息。

- David K

---

即使使用“format long g”，在Octave中它们都是零。 - am304

@am304，既然你有Octave：能试试 A*B' - fliplr(A)*fliplr(B)' 吗？ - Luis Mendo

答案为-1.31512933876365e-009 - 1.06410880107433e-010i，所以并不是零... - am304

@am304 谢谢。这证实了我的猜测，即 OP 中描述的错误是由 Matlab 内部选择不同的操作顺序引起的。 - Luis Mendo

在使用有限精度时，最终会出现误差。记住可能你找到的两个解决方案都是“错误”的是很好的。如果需要更高的精度，可以尝试使用符号工具箱。 - Dennis Jaheruddin

Answer 1

这可能与操作执行的顺序有关。例如，

sum(A.*conj(B)) - fliplr(A)*fliplr(B)'

产生的结果与预期不同

sum(A.*conj(B)) - A*B'

或者更引人注目的是，

A*B' - fliplr(A)*fliplr(B)'

这里的“nonzero result”的意思是返回的结果不为零，并且与你的测试结果同阶。因此，我的猜测是依赖于所使用的方法（sum或*），Matlab 内部执行操作的顺序可能不同，这很可能会导致你观察到的不同舍入误差。

- Luis Mendo

Answer 2

考虑到每个数字的大小，对这个数字进行简单运算时的舍入误差是10^-14级别的。

你有5 * 10^6个数字，因此如果你非常不幸，舍入误差可能会高达5 * 10^-8。

你观察到的误差大小为10^10，处于预期范围内。

请注意，差异不是由于复杂的转置，而是由于求和与矩阵乘积之积造成的。

A = randn(1,5e6)+1i*randn(1,5e6);
B = randn(1,5e6)+1i*randn(1,5e6);

B1 = conj(B); 
B2 = B';

isequal(B1(:),B2(:)) % This returns true

A*transpose(conj(B)) - A*B' % Hence this returns zero

sum(A.*transpose(B')) - A*B' % But this returns something like 1e-10

类似的效果也会发生在非复杂的A和B上：

N=1e6;
A = 1:N; 
B=1:N;

(N * (N + 1) * (2*N + 1))/6 % This will give exactly the right answer
A*B' 
fliplr(A)*fliplr(B)'

请注意，最低的两个答案只相差几百，但它们与正确答案实际上相差超过2000。如果这是一个问题，请考虑使用符号工具箱。那可以让您进行任意精度计算。

- Dennis Jaheruddin