Trace and inverse of big matrices in R and C++

Question

I need to compute two equations related with trace and inverse of a around 30000 x 30000 density matrices. The equations are

-trace( W_i %*% C) 

and

-trace(W_i %*% C %*% W_j C)

I know W_i, W_j and inverse of C. These equations are related with Pearson estimating functions. I am trying to use R and package Matrix, but I couldn't compute the C matrix, using solve() or chol() and chol2inv(). I do not know with is possible using solve() to solve a system of equation and after compute the trace. It is common to use solve function to compute something like C^{-1} W = solve(C, W), but my equation is a little bit different. Any help is welcome. I am thinking about to use RcppArmadillo, but I am not sure that it is able to compute my equations.

Answer 1

You can use RcppArmadillo, but you have to be careful about memory usage. If you save the below code as arma_test.cpp it can be sourced via Rcpp::sourceCpp('w_graph_class.cpp'). Obviously the matrix data is dummy, but it should give you a starting point. Also check out alternative methods to invert C rather than using .i(), such as pinv(), etc.

See this question about large matrices with RcppArmadillo

#include <RcppArmadillo.h>
using namespace Rcpp;

// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]
Rcpp::List arma_calc() {

arma::mat C_inv = arma::mat(30000,30000,arma::fill::randu);
arma::mat W_i = arma::mat(30000,30000,arma::fill::randu);
arma::mat W_j = arma::mat(30000,30000,arma::fill::randu);

double tr_1=-arma::trace(W_i*C_inv.i());

double tr_2=-arma::trace(W_i*C_inv.i()*W_j*C_inv.i());

return Rcpp::List::create(Rcpp::Named("Trace1")=tr_1,Rcpp::Named("Trace2")=tr_2);
}