r语言 - 有没有更快的方法可以在 Rcpp/c++ 中进行这种 Cholesky 分解

Question

这是我的马尔可夫链蒙特卡罗（MCMC）算法中出现的问题。我觉得这个问题很常见，尤其是在分层高斯模型中。因此，如果有更有效的解决方案，那就太好了。所以问题是这样的：

我xi有许多正整数向量，对于从 1 到 n 的i，一个 p.s.d. 矩阵A和一个 p.s.d. 矩阵B。对于每个xi，我都想计算以下 Cholesky 分解：

chol（ kron（ diagmat（ xi）， A ） + B ）

所以多元高斯的协方差矩阵kron( diagmat( xi ), A ) + B，我想从这个高斯中采样，因此需要它的乔列斯基分解。A和B的维度不小，而且我有一个很大的n因此计算上述所有xi的Cholesky分解真的很耗时。以下是我使用 RcppArmadillo 编写的Rcpp函数：

#include <cmath>
#include <Rmath.h>
#include "RcppArmadillo.h"
// [[Rcpp::depends(RcppArmadillo)]]
using namespace arma;
using namespace Rcpp;
// [[Rcpp::export]]
mat test_C(mat A, mat B, mat X){
  int n_x = X.n_cols;
  int d_B = B.n_rows;
  mat sample(d_B, n_x);
  mat mS_chol_inv;
  for (int i = 0; i < n_x; i++){
    mS_chol_inv = inv(trimatu( chol(kron(diagmat( X.col(i) ),A) + B) ));
    sample.col(i) = mS_chol_inv*randn(d_B);
  }
  return(sample);
}

我还使用以下代码将其与 R 对应项进行比较来测试计算效率：

test_R <- function(A,B,X){
  n_x <- ncol(X)
  d_B <- ncol(B)
  res <- sapply(1:n_x, function(x){
    mS_chol <- chol( kronecker( diag(X[,x]),A ) + B )
    return( mS_chol%*%as.matrix( rnorm(d_B) ) )
  })
  return(res)
}
# Simulate Data
R1 <- matrix(rnorm(24*2),24,2)
A <- R1%*%t(R1) + 0.1*diag(24)
R2 <- matrix(rnorm(264*2),264,2)
B <- R2%*%t(R2) + 0.1*diag(264)
X <- matrix(rpois(11*2178, 5),11,2178)
res <- benchmark(res_R <- test_R(A, B, X),
             res_C <- test_C(A, B, X),
             columns=c("test", "replications", "elapsed", "relative"),
             order="relative",
             replications = 2)

结果如下

> print(res) 
                      test replications elapsed relative
1 res_R <- test_R(A, B, X)            2  18.920    1.000
2 res_C <- test_C(A, B, X)            2  20.724    1.095

可以看出，单次运行大约需要 10 秒，这在 MCMC 算法中根本不可行。此外，由于chol()主导了计算复杂性，因此使用Rcpp而不是纯R的改进是微不足道的。但也许我没有写出最有效的代码？那么有什么建议吗？

由于chol()内部的矩阵非常结构化，唯一变化的是xi，也许有一些我不知道的代数技巧可以有效地解决这个问题？我已将其作为数学下的线性代数问题发布，这里是链接。不幸的是，到目前为止，我还没有收到任何解决方案，人们确实指出这是令人尴尬的平行。

任何关于代码/代数的建议都会有所帮助！提前感谢您的时间。

Answer 1

你可以

尝试从这里使用Microsoft R，艺术家以前称为RRO。它与多线程英特尔 MKL 库（查找和安装它的位置相同）集成，并且在英特尔硬件矩阵上操作非常快。

免责声明：YMMV

Answer 2

（Rcpp） Eigen or Armadillo for Cholesky？

我不能将此作为对@zdebruine答案的评论，但我想弄清楚使用 RcppArmadillo/Eigen 的 Cholesky 最佳实践是什么。

注意：这个SO问题是"最快的Cholesky使用犰狳或本征"的顶级Google搜索结果。我认为人们评估两者的优缺点很有用

这是基准测试的第一次尝试（这里的CPU是AMD Ryzen 9 5950X）。我将根据建议进行编辑和更新。

sessionInfo()回报：

R version 4.1.2 (2021-11-01)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.3 LTS
Matrix products: default
BLAS/LAPACK: /opt/intel/compilers_and_libraries_2019.5.281/linux/mkl/lib/intel64_lin/libmkl_rt.so

源：

arma_source <- '
//[[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
//[[Rcpp::export]]
arma::mat arma_chol(const arma::mat& X){
  return arma::chol(X, "lower");
}
'
eigen_source <- '
//[[Rcpp::depends(RcppEigen)]]
#include <RcppEigen.h>
//[[Rcpp::export]]
Eigen::MatrixXd eigen_chol(const Eigen::MatrixXd& X){
  Eigen::LLT<Eigen::MatrixXd> Xllt(X); 
  return Xllt.matrixL();
}
'
Rcpp::sourceCpp(code=arma_source)
Rcpp::sourceCpp(code=eigen_source)

nr <- 2000
Gen <- matrix(rnorm(nr^2), ncol=nr)
X <- crossprod(Gen)
RhpcBLASctl::blas_set_num_threads(1)
RhpcBLASctl::omp_set_num_threads(1)
rbenchmark::benchmark(
  L_arma = arma_chol(X),
  L_eigen = eigen_chol(X),
  U_R = chol(X),
  replications=100)

结果：

     test replications elapsed relative user.self sys.self user.child sys.child
1  L_arma          100   4.770    1.056     4.635    0.136          0         0
2 L_eigen          100  12.908    2.856    12.768    0.141          0         0
3     U_R          100   4.519    1.000     4.403    0.116          0         0

我不熟悉使特征代码更有效的方法。我觉得很奇怪，这完全是有区别的。如果我允许多线程，差异就会扩大，因为 Eigen 在我当前的配置中只在 1 个线程上运行。R的chol()与arma相同这一事实并不奇怪，因为两者都使用MKL的BLAS的BLA。

如果我们添加#define EIGEN_USE_BLAS我们将得到以下内容

     test replications elapsed relative user.self sys.self user.child sys.child
1  L_arma          100   5.556    1.039     4.958    0.597          0         0
2 L_eigen          100   6.257    1.170     5.621    0.631          0         0
3     U_R          100   5.347    1.000     4.555    0.790          0         0

这几乎使所有内容都使用相同的 BLAS 函数

Answer 3

改用 RcppEigen 将使您受益匪浅。

这是因为：

• 本征有一个更快的乔列斯基求解器
• Eigen 允许通过引用进行子视图，并具有从这些子视图求解 Cholesky 的内置优化
• Eigen 有大量专门针对 LLT 分解的底层优化，您可以在不知情的情况下
从中受益。

就您可以期望的性能提升而言，我只能谈谈我的用例，我使用 LLT 分解分解高达 100x100 的矩阵，但犰狳几乎慢了 50%！

看看本征和犰狳的源代码！犰狳通过通用模板运行每个cholesky到LAPACK/BLAS"辅助"库。Eigen 在内部运行决策树以选择优化的实现，然后尽可能通过参考内存中的内部结构直接求解。

dr：对于分解，你不会轻易击败艾根。