从给定的稀疏矩阵中提取对角矩阵

Question

使用RcppEigen我只想提取稀疏矩阵的对角线作为稀疏矩阵。看起来很简单 - 下面是我的尝试，none 提供了我想要的结果。请注意，尝试 5 无法编译且无法正常工作。这是我使用的一些资源； Rcpp Gallery, KDE Forum and in the same post KDE Forum (2), Eigen Sparse Tutorial and 。感觉自己很接近了。。。可能不是。。。我让专家来决定。

// [[Rcpp::depends(RcppEigen)]]
#include <RcppEigen.h>
#include <Eigen/SparseCore>

// [[Rcpp::export]]
Eigen::SparseMatrix<double> diag_mat1(Eigen::Map<Eigen::SparseMatrix<double> > &X){
  // cannot access diagonal of mapped sparse matrix
  const int n(X.rows());
  Eigen::VectorXd dii(n);
  for (int i = 0; i < n; ++i) {
     dii[i] = X.coeff(i,i);
  }
  Eigen::SparseMatrix<double> ans(dii.asDiagonal());
  return ans;
}

// [[Rcpp::export]]
Eigen::SparseMatrix<double> diag_mat2(Eigen::SparseMatrix<double> &X){
  Eigen::SparseVector<double> dii(X.diagonal().sparseView());
  Eigen::SparseMatrix<double> ans(dii);
  return ans;
}

// [[Rcpp::export]]
Eigen::SparseMatrix<double> diag_mat3(Eigen::SparseMatrix<double> &X){
  Eigen::VectorXd dii(X.diagonal());
  Eigen::SparseMatrix<double> ans(dii.asDiagonal());
  ans.pruned(); //hoping this helps
  return ans;
}

// [[Rcpp::export]]
Eigen::SparseMatrix<double> diag_mat4(Eigen::SparseMatrix<double> &X){
  Eigen::SparseMatrix<double> ans(X.diagonal().asDiagonal());
  return ans;
}

// [[Rcpp::export]]
Eigen::SparseMatrix<double> diag_mat5(Eigen::SparseMatrix<double> &X){
  struct keep_diag {
    inline bool operator() (const int& row, const int& col, const double&) const
    { return row==col; }
  };
  Eigen::SparseMatrix<double> ans(X.prune(keep_diag()));
  return ans;
}


/***R
library(Matrix)
set.seed(42)
nc <- nr <- 5
m  <- rsparsematrix(nr, nc, nnz = 10)
diag_mat1(m)
diag_mat2(m)
diag_mat3(m)
diag_mat4(m)

*/

编辑：添加了每次尝试给出的结果；

> diag_mat1(m)
5 x 5 sparse Matrix of class "dgCMatrix"

[1,] 0  .     . . .  
[2,] . -0.095 . . .  
[3,] .  .     0 . .  
[4,] .  .     . 2 .  
[5,] .  .     . . 1.5
> diag_mat2(m)
5 x 1 sparse Matrix of class "dgCMatrix"

[1,]  .    
[2,] -0.095
[3,]  .    
[4,]  2.000
[5,]  1.500
> diag_mat3(m)
5 x 5 sparse Matrix of class "dgCMatrix"

[1,] 0  .     . . .  
[2,] . -0.095 . . .  
[3,] .  .     0 . .  
[4,] .  .     . 2 .  
[5,] .  .     . . 1.5
> diag_mat4(m)
5 x 5 sparse Matrix of class "dgCMatrix"

[1,] 0  .     . . .  
[2,] . -0.095 . . .  
[3,] .  .     0 . .  
[4,] .  .     . 2 .  
[5,] .  .     . . 1.5

EDIT2：添加了所需的输出；

5 x 5 sparse Matrix of class "dgCMatrix"               
[1,] .  .     . . .  
[2,] . -0.095 . . .  
[3,] .  .     . . .  
[4,] .  .     . 2 .  
[5,] .  .     . . 1.5

回答灵感来自 Aleh；

Eigen::SparseMatrix<double> diag_mat6(Eigen::Map<Eigen::SparseMatrix<double> > &X){
  const int n(X.rows());
  Eigen::SparseMatrix<double> dii(n, n);
  for (int i = 0; i < n; ++i) {
    if (X.coeff(i,i) != 0.0 ) dii.insert(i, i) = X.coeff(i,i);
  }
  dii.makeCompressed();
  return dii;
}

Answer 1

我更喜欢 RcppArmadillo 因为它通常表现得比 RcppEigen 更像 R。

对于您的问题，使用 RcppArmadillo，您可以：

// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>

// [[Rcpp::export]]
arma::sp_mat extractDiag(const arma::sp_mat& x) {

  int n = x.n_rows;
  arma::sp_mat res(n, n);

  for (int i = 0; i < n; i++)
    res(i, i) = x(i, i);

  return res;
}

根据@mtall 的建议，您可以简单地使用：

// [[Rcpp::export]]
arma::sp_mat extractDiag3(const arma::sp_mat& x) {
  return arma::diagmat(x);
}

如果你真的想在 Eigen 中这样做，来自 the documentation，我想出了：

// [[Rcpp::export]]
Eigen::SparseMatrix<double> extractDiag2(Eigen::Map<Eigen::SparseMatrix<double> > &X){

  int n = X.rows();
  Eigen::SparseMatrix<double> res(n, n);
  double d;

  typedef Eigen::Triplet<double> T;
  std::vector<T> tripletList;
  tripletList.reserve(n);
  for (int i = 0; i < n; i++) {
    d = X.coeff(i, i);
    if (d != 0) tripletList.push_back(T(i, i, d));
  }
  res.setFromTriplets(tripletList.begin(), tripletList.end());

  return res;
}

Answer 2

我认为你只需要跳过对角线上的零个元素：

 for (int i = 0; i < n; ++i) {
     if (X.coeff(i,i) != 0.0)
        dii[i] = X.coeff(i,i);
     }
  }

从给定的稀疏矩阵中提取对角矩阵

Extract diagonal matrix from a given sparse matrix

r

rcpp

eigen