在Matlab中加速计算inv(X'*X)*Q*inv(X'*X)？

Question

我需要为大型多元回归模型计算Newey-West标准误差。

计算的最后一步是获得：

nwse = sqrt(diag(N.*inv(X'*X)*Q*inv(X'*X)));

这个文件交换贡献将其实现为：

nwse = sqrt(diag(N.*((X'*X)\Q/(X'*X))));  

这看起来很合理，但对于我来说（5000x5000的稀疏矩阵Q和X'* X），它对我的需求来说太慢了（约30秒，我需要为大约一百万个不同的模型重复此操作）。有什么办法可以使这行代码更快吗？

请注意，我只需要对角线，而不是整个矩阵，并且Q和（X'* X）都是正定的。

- Trisoloriansunscreen

---

只是随机尝试一下...您能否评估以下表达式并验证两件事：1.它给出了相同的结果。2.它显着提高了速度。((X'X)*((sparseinv(Q)*X')*X)*X。从这里获取sparseinv。 此外，您可以尝试在math.stackexchange.com上发布此问题。确保您以数学方式表述。可能有一些技术可以利用稀疏性和正定属性。 - Autonomous

@Robert P.：大约每100个模型有相同的X，所以我可以缓存inv(X'*X)（还有更优雅的方法吗？） - Trisoloriansunscreen

使用 inv 从来不是正确的方法... 我稍后会回复你 =) X 是一个向量，对吗？ - Stewie Griffin

@Robert P.：X是一个矩阵。X'*X与Q具有相同的维度。 - Trisoloriansunscreen

另一个选择是Cholesky分解。 - Autonomous

Answer 1

我认为你可以通过显式进行LU分解来节省大量计算时间，[l, u, p, q] = lu(X'*X);并在进行计算时使用这些因子。此外，由于X在约100个模型中都是常量，预先计算X'*X很可能会节省你一些时间。

请注意，在你的情况下，最耗时的操作很可能是sqrt函数。

% Constant for every 100 models or so:
M = X'*X;
[l, u, p, q] = lu(M);

% Now, I guess this should be quite a bit faster (I might have messed up the order):
nwse = sqrt(diag(N.*(q * ( u \ (l \ (p * Q))) * q * (u \ (l \ p)))));

前两个术语是常用的： l：下三角矩阵 u：上三角矩阵

现在，p和q则不太常见。 p是行置换矩阵，用于获得数字稳定性。对于稀疏矩阵，使用[l、u、p] = lu(M)与使用[l、u] = lu(M)相比，性能提升不大。

然而，q可以显著提高性能。q是列置换矩阵，用于在进行因式分解时减少填充量。

请注意，[l、u、p、q] = lu(M)仅适用于稀疏矩阵的有效语法。

---

关于为什么如上所述使用完全主元消元法应该更快：尝试以下操作以查看列置换矩阵 q 的目的。围绕对角线对齐的元素更易处理。

S = sprand(100,100,0.01);
[l, u, p] = lu(S);
spy(l)
figure
spy(u)

现在，将其与此进行比较：

[ll, uu, pp, qq] = lu(S);
spy(ll);
figure
spy(uu);

很不幸，我现在没有MATLAB，所以无法保证我将所有参数放置在正确的顺序中，但我认为是正确的。

- Stewie Griffin

Answer 2

在参考了Robert_P的有益答案和Parag的评论后，我发现以下方法对于我的大规模稀疏数据是最快的：

L=chol(X'*X,'lower'); L=full(L);     
invXtX = L'\(L\ speye(size(X,2))); 
nwse = sqrt(N.*sum(invXtX.*(Q*invXtX)));  

最后一行有效地计算了对角线，灵感来自这里。

- Trisoloriansunscreen

---

为什么要执行 L = full(L)？如果跳过它不是更快吗？ - Stewie Griffin

很奇怪，但稀疏矩阵稍微慢一些。 - Trisoloriansunscreen

嗯，这并不一定奇怪，这取决于矩阵的稀疏程度。将完整的矩阵表示为稀疏矩阵将非常非常慢，因此如果矩阵不是非常稀疏，我认为这样做是正确的。我猜平衡点应该在20％的非零值左右，但这只是一个猜测。 - Stewie Griffin

你尝试过以下方法吗？使用稀疏矩阵计算invXtX，然后在执行invXtX = full(invXtX)。也就是说，在计算平方根时只使用完整矩阵。 - Stewie Griffin

我尝试过了，但效果并不好。我猜在我的特定情况下，L 矩阵并不是那么稀疏。 - Trisoloriansunscreen