为什么x\y比(x'*x)\(x'*y)慢得多？

Question

对于NxP矩阵x和yN> 1的Nx1向量,这两个表达式

x \ y                                       -- (1)

和

(x' * x) \ (x' * y)                         -- (2)

两者都计算b矩阵方程的解

x * b = y

在最小二乘意义上,即使数量

norm(y - x * b)

最小化.表达式(2)使用经典算法来解决普通最小二乘回归,其中\运算符的左手参数是正方形.它相当于写作

inv(x' * x) * (x' * y)                      -- (3)

但它使用的算法在数值上更稳定.事实证明,(3)比(2)适度快,即使(2)不必产生作为副产物的逆矩阵,但我可以接受给定额外的数值稳定性.

然而,一些简单的时间(N = 100,000和P = 30)表明表达式(2)比表达式(1)快5倍以上,即使(1)具有更大的灵活性来选择所使用的算法!例如,对(1)的任何调用都可以调度X的大小,在N> P的情况下,它可以减少到(2),这会增加很少的开销,但肯定不会花费5倍更长的时间.

表达式(1)中发生了什么导致它花了这么长时间？

---

编辑:这是我的时间

x = randn(1e5, 30);
y = randn(1e5,1);
tic, for i = 1:100; x\y; end; t1=toc;
tic, for i = 1:100; (x'*x)\(x'*y); end; t2=toc;
assert( abs(norm(x\y) - norm((x'*x)\(x'*y))) < 1e-10 );
fprintf('Speedup: %.2f\n', t1/t2)

加速:5.23

Answer 1

你知道在你的测试中

size(x) == [1e5  30]    but   size(x'*x) == [30  30]
size(y) == [1e5   1]    but   size(x'*y) == [30   1]

这意味着进入mldivide函数的矩阵大小不同4个数量级!这将导致确定使用哪种算法相当大和重要的任何开销(并且可能还在两个不同的问题上运行相同的算法).

换句话说,你有一个有偏见的测试.为了进行公平的测试,请使用类似的东西

x = randn(1e3);
y = randn(1e3,1);

我发现(最差的5次运行):

Speedup: 1.06  %// R2010a
Speedup: 1.16  %// R2010b
Speedup: 0.97  %// R2013a

......差异几乎消失了.

但是,这确实表明,如果你确实有一个与观察数量相比维度较低的回归问题,那么首先进行乘法真的是值得的:)

mldivide是一个万能的,真的很棒.但通常,有知识有关的问题可能做出更具体的解决方案,如预乘,预调节,lu,qr,linsolve,幅度等下单要快.

Answer 2

尽管（1）有更大的灵活性来选择所使用的算法！例如，对 (1) 的任何调用都可以只按 X 的大小进行调度，在 N>P 的情况下，它可以减少到 (2)，这会增加少量的开销，但肯定不会花费 5 倍更长。

不是这种情况。选择使用哪种算法可能会花费大量开销，特别是与诸如此类的相对较小输入的计算相比。在这种情况下，因为 MATLAB 可以看到您有x'*x，所以它知道其中一个参数必须既是方形的又是对称的（是的 - 即使在解析器级别，线性代数的知识也内置在 MATLAB 中），并且可以直接调用中适当的代码路径之一\。

我不能说这是否完全解释了您所看到的时间差异。我想进一步调查，至少通过：

1. 确保将代码放入函数中，并预热该函数以确保 JIT 参与 - 然后尝试相同的操作feature('accel', 'off')来消除 JIT 的影响
2. 在更大范围的输入大小上进行尝试，以检查“算法选择开销”与计算时间相比有何贡献。