R中计算大矩阵逆的最快方法

Question

我需要计算一个帽子矩阵（如线性回归）。标准 R 代码将是：

H <- tcrossprod(tcrossprod(X, solve(crossprod(X))), X)

用X是相对大的矩阵（即1E5 * 100），并且该线具有运行数千次。我知道最大的限制部分是逆计算，但叉积也可能很耗时。有没有更快的替代方法来执行这些矩阵运算？我尝试了 Rcpp 并查看了几篇文章，但我测试的任何替代方法都较慢。也许我没有正确编写 C++ 代码，因为我不是高级 C++ 程序员。

谢谢！

Answer 1

逐行跟踪代码有点困难，因为 R 代码的设置有点复杂。但请继续阅读，下面的指针。

重要的部分是这个话题已经讨论过很多次了：发生的事情是 R 将它分派到BLAS（基本线性代数子程序）和LAPACK（线性代数 PACKage）库。其中包含人类已知的最有效的代码。一般来说，你不能通过重写来获得它。

通过将一个 BLAS/LAPACK 实现切换为另一个实现，可以获得性能差异——网上也有很多很多帖子。R 本身带有所谓的“参考 BLAS”，已知是正确的，但最慢。您可以切换到 Atlas、OpenBLAS、MKL 等，具体取决于您的操作系统；安装随附的一些R 手册中提供了有关如何执行此操作的说明。

为了完整起见，每个文件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); 

这是这个 LAPACK 函数。对于所有其他人来说，它本质上是相同的，因此这不太可能成为您进行优化的场所。