两个矩阵之和的反函数

Question

我想实现一个代码来计算两个矩阵之和的逆。我的算法是递归的，我需要使用循环for()我试图在R，但我的代码是非常缓慢的。然后，我正在尝试使用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×15000

最佳

Answer 1

您是否尝试过解决（）？

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 

Answer 2

我更舒服征和可以得到一些速度，而不会改变算法：

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