c++ - R 和 C++ 中大矩阵的跟踪和逆

Question

我需要计算两个与大约 30000 x 30000 密度矩阵的迹和逆相关的方程。方程是

-trace( W_i %*% C) 

和

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

我知道 W_i、W_j 和 C 的逆。这些方程与 Pearson 估计函数有关。我正在尝试使用 R 和包 Matrix，但我无法计算 C 矩阵，使用 solve() 或 chol() 和chol2inv()。我不知道是否可以使用 solve() 求解方程组并计算迹线。通常使用 solve 函数来计算类似 C^{-1} W = solve(C, W) 的东西，但我的等式有点不同。欢迎任何帮助。我正在考虑使用 RcppArmadillo，但我不确定它是否能够计算我的方程。

Answer 1

你可以使用 RcppArmadillo，但你必须小心内存使用。如果将以下代码保存为 arma_test.cpp，则可以通过 Rcpp::sourceCpp('w_graph_class.cpp') 获取它。显然矩阵数据是虚拟的，但它应该给你一个起点。另请检查替代方法来反转 C 而不是使用 .i()，例如 pinv() 等。

查看此 question 关于使用 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);
}