Обратная сумма двух матриц

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