我必须计算大型多元回归模型的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)都是正定的。
I have to calculate Newey-West standard errors for large multiple regression models.
The final step of this calculation is to obtain
nwse = sqrt(diag(N.*inv(X'*X)*Q*inv(X'*X)));This file exchange contribution implements this as
nwse = sqrt(diag(N.*((X'*X)\Q/(X'*X))));This looks reasonable, but in my case (5000x5000 sparse Q and X'*X) it's far too slow for my needs (about 30secs, I have to repeat this for about one million different models). Any ideas how to make this line faster?
Please note that I need only the diagonal, not the entire matrix and that both Q and (X'*X) are positive-definite.
最满意答案
我相信你可以通过明确地进行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, p] = 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,所以我不能保证我把所有的参数都按正确的顺序排列,但我认为这是正确的。
I believe you can save a lot of computation time by explicitly doing an LU factorization, [l, u, p, q] = lu(X'*X); and use those factors when doing the calculations. Also, since X are constant for about 100 models, pre-calculating X'*X will most likely save you some time.
Note that in your case, the most time demanding operation might very well be the sqrt-function.
% 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)))));The first two terms are commonly used:
l the lower triangular matrix
u the upper triangular matrix
Now, p and q are a bit more uncommon.
p is a row permutation matrix used to obtain numerical stability. There is not much performance gain in doing [l, u, p] = lu(M) compared to [l, u] = lu(M) for sparse matrices.
q however offers a significant performance gain. q is a column permutation matrix that is used to reduce the amount of fill when doing the factorization.
Note that the [l, u, p, q] = lu(M) is only a valid syntax for sparse matrices.
As for why using full pivoting as described above should be faster:
Try the following to see the purpose of the column permutation matrix q. It is easier to work with elements that are aligned around the diagonal.
S = sprand(100,100,0.01); [l, u, p] = lu(S); spy(l) figure spy(u)Now, compare it with this:
[ll, uu, pp, qq] = lu(S); spy(ll); figure spy(uu);Unfortunately, I don't have MATLAB here right now, so I can't guarantee that I put all the arguments in the correct order, but I think it's correct.
更多推荐
发布评论