我需要计算一个帽子矩阵(来自线性回归)。标准 R 代码为:
H <- tcrossprod(tcrossprod(X, solve(crossprod(X))), X)
与 X
作为一个相对较大的矩阵(即 1e5*100),这条线必须运行数千次。我知道最受限制的部分是逆计算,但叉积也可能很耗时。有没有更快的替代方法来执行这些矩阵运算?我尝试了 Rcpp 并查看了几篇文章,但我测试的任何替代方案都比较慢。也许我没有正确编写我的 C++ 代码,因为我不是高级 C++ 程序员。
谢谢!
最佳答案
逐行追逐此代码有点困难,因为 R 代码的设置有点复杂。但请继续阅读下面的指示。
重要的是这个话题已经被讨论了很多次:发生的事情是 R 将它发送到 BLAS (Basic Linear Algebra Subprogram)。和 LAPACK (Linear Algebra PACKage)图书馆。其中包含人类已知的最高效的代码。通常,您不能通过重写来获得 yield 。
可以通过将一种 BLAS/LAPACK 实现切换为另一种来获得性能差异——网上也有很多很多帖子。 R 本身带有已知正确但最慢的所谓“引用 BLAS”。您可以切换到 Atlas、OpenBLAS、MKL 等,具体取决于您的操作系统;关于如何操作的说明在一些 R manuals 中。安装时附带的。
为了完整性,每个文件 src/main/names.c
命令 %*%
、crossprod
和 tcrossprod
均引用do_matprod
。这在文件 src/main/array.c 中并根据参数类型进行大量参数检查、排列和分支,但例如一个路径然后调用
F77_CALL(dsyrk)(uplo, trans, &nc, &nr, &one, x, &nr, &zero, z, &nc
FCONE FCONE);
这是this LAPACK function .对于所有其他人来说,它本质上是相同的,因此这不太可能成为您优化的场所。
关于r - R中计算大矩阵逆的最快方法,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/64501106/