在 Rcpp/c++ 中有更快的方法来进行 Cholesky 分解吗？
Is there a faster way to do this Cholesky factorization in Rcpp/c++?
这是我的马尔可夫链 Monte Carlo (MCMC) 算法中出现的问题。我觉得这个问题很常见，尤其是在分层高斯模型中。因此，如果有更有效的解决方案，那就太好了。所以问题是这样的：
我有很多正整数向量xi，对于i从1到n，一个p.s.d。矩阵 A 和一个 p.s.d。矩阵 B。对于每个 xi，我想计算以下 Cholesky 分解：
chol( kron( diagmat( xi ), A ) + B )
所以 kron( diagmat( xi ), A ) + B 是多元高斯的协方差矩阵，我想从这个高斯中采样，因此需要它的 Cholesky 分解。 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，也许有一些我不知道的代数技巧可以有效地解决这个问题？我已将其作为线性代数问题发布在 数学 下，这里是 link。不幸的是，到目前为止我还没有收到任何解决方案，人们确实指出这是令人尴尬的并行。
code/algebra 上的任何建议都会有所帮助！提前感谢您的宝贵时间。
您可以尝试使用来自 here 的 Microsoft R，前身为 RRO 的艺术家。它集成了多线程 Intel MKL 库（在同一个地方可以找到并安装它）并且在 Intel 硬件上矩阵操作非常快。
免责声明：YMMV
你会从使用 RcppEigen 中受益匪浅。
这是因为：
- Eigen 具有更快的 Cholesky 求解器
- Eigen 允许通过引用子视图，并内置优化从这些子视图求解 Cholesky
- Eigen 有大量专门针对 LLT 分解的底层优化，您可以在不知情的情况下从中受益。
就您可以预期的性能提升而言，我只能谈谈我的用例，我使用 LLT 分解将矩阵分解到 100x100，但犰狳几乎 50% 慢!
查看 Eigen 和 Armadillo 的源代码！ Armadillo 通过通用模板将每个 cholesky 运行到 LAPACK/BLAS“辅助”库。 Eigen在内部运行一个决策树来选择一个优化的实现，然后尽可能直接参考内存中的内部结构来求解。
tl;dr：对于分解，你不会轻易击败 Eigen。
(Rcpp) Cholesky 的 Eigen 或犰狳？
我不能post这是对@zdebruine 的回答的评论，但我想弄清楚 Cholesky 使用 RcppArmadillo/Eigen 的最佳做法是什么。
注意：这个 SO 问题是“使用 Armadillo 或 Eigen 的最快 Cholesky”搜索结果中排名靠前的 Google。我认为这对人们评估两者的优缺点很有用
这是基准测试的第一次尝试（CPU 这里是 AMD Ryzen 9 5950X）。我会根据建议进行编辑和更新。
sessionInfo() returns:
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 代码效率的方法。我觉得奇怪，完全不同。如果我允许多线程，差异会增加，因为 Eigen 仅在我当前配置中的 1 个线程上运行。
R 的 chol() 与 arma 的相同这一事实不足为奇，因为两者都使用 MKL 的 BLAS
如果我们添加 #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 函数
这是我的马尔可夫链 Monte Carlo (MCMC) 算法中出现的问题。我觉得这个问题很常见，尤其是在分层高斯模型中。因此，如果有更有效的解决方案，那就太好了。所以问题是这样的：
我有很多正整数向量xi，对于i从1到n，一个p.s.d。矩阵 A 和一个 p.s.d。矩阵 B。对于每个 xi，我想计算以下 Cholesky 分解：
chol( kron( diagmat( xi ), A ) + B )
所以 kron( diagmat( xi ), A ) + B 是多元高斯的协方差矩阵，我想从这个高斯中采样，因此需要它的 Cholesky 分解。 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，也许有一些我不知道的代数技巧可以有效地解决这个问题？我已将其作为线性代数问题发布在 数学 下，这里是 link。不幸的是，到目前为止我还没有收到任何解决方案，人们确实指出这是令人尴尬的并行。
code/algebra 上的任何建议都会有所帮助！提前感谢您的宝贵时间。
您可以尝试使用来自 here 的 Microsoft R，前身为 RRO 的艺术家。它集成了多线程 Intel MKL 库（在同一个地方可以找到并安装它）并且在 Intel 硬件上矩阵操作非常快。
免责声明：YMMV
你会从使用 RcppEigen 中受益匪浅。
这是因为：
- Eigen 具有更快的 Cholesky 求解器
- Eigen 允许通过引用子视图，并内置优化从这些子视图求解 Cholesky
- Eigen 有大量专门针对 LLT 分解的底层优化，您可以在不知情的情况下从中受益。
就您可以预期的性能提升而言，我只能谈谈我的用例，我使用 LLT 分解将矩阵分解到 100x100，但犰狳几乎 50% 慢!
查看 Eigen 和 Armadillo 的源代码！ Armadillo 通过通用模板将每个 cholesky 运行到 LAPACK/BLAS“辅助”库。 Eigen在内部运行一个决策树来选择一个优化的实现，然后尽可能直接参考内存中的内部结构来求解。
tl;dr：对于分解，你不会轻易击败 Eigen。
(Rcpp) Cholesky 的 Eigen 或犰狳？
我不能post这是对@zdebruine 的回答的评论，但我想弄清楚 Cholesky 使用 RcppArmadillo/Eigen 的最佳做法是什么。
注意：这个 SO 问题是“使用 Armadillo 或 Eigen 的最快 Cholesky”搜索结果中排名靠前的 Google。我认为这对人们评估两者的优缺点很有用
这是基准测试的第一次尝试（CPU 这里是 AMD Ryzen 9 5950X）。我会根据建议进行编辑和更新。
sessionInfo() returns:
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 代码效率的方法。我觉得奇怪，完全不同。如果我允许多线程，差异会增加，因为 Eigen 仅在我当前配置中的 1 个线程上运行。
R 的 chol() 与 arma 的相同这一事实不足为奇，因为两者都使用 MKL 的 BLAS
如果我们添加 #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 函数