两个矩阵之和的逆
Inverse of sum of two matrices
我正在尝试实现一个代码来计算两个矩阵之和的逆。我的算法是递归的,我需要使用我在R中尝试过的循环for()
,但我的代码非常慢。然后,我尝试使用RcppArmadillo,但我的代码非常非常慢。我觉得我做错了什么。让我展示一下我的R代码。
mySolveR <- function(A,B){
ncol = dim(B)[1]
ZERO.B <- Matrix(0,ncol = ncol, nrow = ncol)
invCi <- A
for(i in 1:ncol){
ZERO.B[,i] <- B[,i]
gi <- 1/(1 + sum(diag(ZERO.B%*%invCi)))
invCi <- invCi - gi*(invCi%*%ZERO.B%*%invCi)
ZERO.B[,i] <- 0
}
return(invCi)}
现在我的C++代码使用RcppArmadillo。
src <- '
Rcpp::NumericMatrix Ac(A); // creates Rcpp matrix from SEXP
Rcpp::NumericMatrix Bc(B);
int n = Ac.nrow(), k = Ac.ncol();
arma::mat A(Ac.begin(), n, k, false); // reuses memory and avoids extra copy
arma::mat B(Bc.begin(), n, k, false);
arma::mat Z(n,k);
Z.zeros();
arma::mat invCi = A;
for( int i = 0 ; i < n ; i++){
Z.col(i) = B.col(i);
double gi = 1/(1 + trace(Z*invCi));
invCi = invCi - gi*(invCi*Z*invCi);
Z.zeros() ;
}
return wrap(invCi);'
我正在使用内联包来编译我的函数。
mySolveCpp <- cxxfunction(signature(A = "numeric", B = "numeric"),
src, plugin="RcppArmadillo")
现在考虑以下简单的例子,
A <- diag(5)
B <- matrix(c(1,-1,0,0,0, -1, 2, -1,0,0, 0,-1,2,-1,0,
0,0,-1,2,-1, 0,0,0,-1,1),5,5)
用我的函数计算A+B 的倒数
mySolveCpp(A,B)
mySolveR(A,B)
在这个小示例中,您可以看到我的函数工作得很好。但我想把这个算法应用于一个大约15000 x 15000的矩阵。在这种情况下,我的R代码不起作用,而我的C++代码非常慢,需要花费数小时来计算倒数。我想知道是否有可能改进我的C++代码来处理大矩阵,如15000 x 15000。
最佳
您尝试过solve()吗?
A <- diag(5)
B <- matrix(c(1,-1,0,0,0, -1, 2, -1,0,0, 0,-1,2,-1,0,0,0,-1,2,-1, 0,0,0,-1,1),5,5)
solve(A+B)
对于稀疏矩阵对象:
As=Matrix(A)
Bs=Matrix(B)
solve(As+Bs)
5 x 5 Matrix of class "dsyMatrix"
[,1] [,2] [,3] [,4] [,5]
[1,] 0.61818182 0.23636364 0.09090909 0.03636364 0.01818182
[2,] 0.23636364 0.47272727 0.18181818 0.07272727 0.03636364
[3,] 0.09090909 0.18181818 0.45454545 0.18181818 0.09090909
[4,] 0.03636364 0.07272727 0.18181818 0.47272727 0.23636364
[5,] 0.01818182 0.03636364 0.09090909 0.23636364 0.61818182
我对Eigen更满意,并且可以在不改变算法的情况下获得一些加速:
src2 <- '
using Eigen::Map;
using Eigen::MatrixXd;
using Rcpp::as;
const Map<MatrixXd> A(as<Map<MatrixXd> >(AA));
const Map<MatrixXd> B(as<Map<MatrixXd> >(BB));
const int n = A.rows(), k = A.cols();
MatrixXd Z(n,k), C(n,k);
const MatrixXd Z0 = Z.setZero();
MatrixXd invCi = A;
double gi;
for( int i = 0 ; i < n ; i++){
Z.col(i) = B.col(i);
C = Z*invCi;
gi = 1/(1 + C.trace());
invCi -= gi*(invCi*C);
Z=Z0;
}
return wrap(invCi);'
mySolveCpp2 <- cxxfunction(signature(AA = "matrix", BB = "matrix"),
src2, plugin="RcppEigen")
set.seed(42)
A <- matrix(rnorm(1e4), 1e2)
B <- matrix(rnorm(1e4), 1e2)
all.equal(
mySolveCpp(A,B),
mySolveCpp2(A,B))
#[1] TRUE
library(microbenchmark)
microbenchmark(mySolveCpp(A,B),
mySolveCpp2(A,B), times=10)
#Unit: milliseconds
# expr min lq median uq max neval
# mySolveCpp(A, B) 129.51222 129.62216 132.68336 136.67307 137.43591 10
# mySolveCpp2(A, B) 46.76913 47.26311 47.96435 50.12505 61.82288 10
相关文章:
- 通过插槽和信号在不同线程中的两个qt对象之间进行通信
- 为什么具有静态存储持续时间的同一内联变量在包含在 VS2017 编译的两个翻译单元中时会构造和销毁两次
- 大小为 4*2 和 2*4 的两个矩阵的矩阵乘法
- 在C++中将返回unique_ptr和shared_ptr的两个工厂方法组合为一个
- 一个 ui 成员的 Qt 连接和同一连接中的两个信号
- 查找和比较两个字符串中的字符
- 具有单个函数的两个点的距离和坡度
- 使用聚合创建和关联两个不同的对象成员
- 给定仅包含布尔类型成员的结构的两个对象 s1 和 s2,只要 s1 的成员为 true,请检查 s2 的每个成员是否为真
- 比较C或C++中浮点值的两个和
- 将或(||)运算符应用于只有0和1的两个向量
- 在 GCC 4.6 和 4.7 上模板模板扣除的两个不同结果
- 调试:跟踪(和diffing)同一程序的两个版本的函数调用树
- 同一运算符和类的两个重载函数
- 在整数数据类型和连续两个字符类型之后.第 2 个字符的数据类型跳过..为什么
- 具有相同名称和不同模板参数的两个结构如何工作
- 为什么使用带有 gcc 和模板函数的两个或多个基于 typedef 的静态断言时会出现"conflicting declaration"错误?
- 使用Qt的两个按钮和标签之间的基本交互
- 使用类和标头的两个错误
- 如何在C++中为矢量数据结构中的二维数组获取用户输入,而不使用行和列上的两个for循环