迭代 MappedSparseMatrix 特征中的非零元素
Iterating non-zero elements in MappedSparseMatrix eigen
我正在将一个小的稀疏矩阵(用于测试)传递给 R 中的 C++ 函数。该矩阵属于 class dgCMatrix,如下所示:
5 x 5 sparse Matrix of class "dgCMatrix"
[1,] . . . . .
[2,] 1 1 . . .
[3,] . . . . .
[4,] . . 1 . .
[5,] . 1 . . .
我正在按照文档 here 中的说明迭代此矩阵。
我的函数打印出迭代器的值和行索引、列索引。
C++函数定义如下:
#include <RcppEigen.h>
// [[Rcpp::depends(RcppEigen)]]
using Eigen::MappedSparseMatrix;
using Eigen::SparseMatrix;
using Eigen::VectorXi;
using Eigen::Map;
using namespace Rcpp;
using namespace std;
// [[Rcpp::export]]
void createRec(RObject sparse_mat, IntegerVector sparse_vec) {
const MappedSparseMatrix<int> spmat(as<MappedSparseMatrix<int> >(sparse_mat));
long int nrow = spmat.rows();
long int ncol = spmat.cols();
NumericVector sim(nrow);
for(int k=0;k<spmat.outerSize();k++){
for(SparseMatrix<int,Eigen::ColMajor>::InnerIterator it(spmat,k);it;++it){
cout<<"k="<<k<<endl;
cout<<"value="<<it.value()<<endl;
cout<<"it.row="<<it.row()<<endl;
cout<<"it.col="<<it.col()<<endl;
cout<<"index="<<it.index()<<endl;
}
}
}
对于上面给出的矩阵,打印了以下结果:
k=0
value=156148016
it.row=66211520
it.col=0
index=66211520
k=1
value=0
it.row=0
it.col=1
index=0
k=1
value=1
it.row=4
it.col=1
index=4
k=2
value=1
it.row=3
it.col=2
index=3
1.) k=0对应的值有什么解释吗?这些可能是由于以错误的方式传递矩阵造成的吗?
2.) k 在 outerSize 上迭代,它等于 5,为什么它不迭代 k=3,4?考虑到它是一个稀疏矩阵,迭代器会出现这种行为。
每当您看到像 156148016
或 66211520
这样的非常大的数字时,很可能是您遇到了未定义的行为 (UB) 或某个值未正确初始化。在这种情况下,它是后者。具体来说,dgCMatrix
class' 的基础类型是 double
而不是 int
.
The dgCMatrix
class is a class of sparse numeric matrices in the compressed, sparse, column-oriented format. In this implementation the non-zero elements in the columns are sorted into increasing row order. dgCMatrix
is the "standard" class for sparse numeric matrices in the Matrix
package.
因此,当您尝试创建到底层 RObject
内存位置的映射时,需要一个额外的步骤来在请求的 不同[=88] 中重新创建对象=] 类型。添加 const
项后,我敢打赌这些条目会按预期进行,因为编译器可能会在内存中保留中间对象。
所以,变化如下:
MappedSparseMatrix<int> spmat(as<MappedSparseMatrix<int> >(sparse_mat));
至:
MappedSparseMatrix<double> spmat(as<MappedSparseMatrix<double> >(sparse_mat));
应该足够了。
linked example使用了SparseMatrix
矩阵,这里你使用的是MappedSparseMatrix
但没有为第二个循环设置合适的MappedSparseMatrix::InnerIterator
。
因此,我们有:
for(SparseMatrix<int,Eigen::ColMajor>::InnerIterator it(spmat,k);it;++it){
要去:
for(MappedSparseMatrix<double>::InnerIterator it(spmat,k);it;++it){
另外,请注意 SparseMatrix<int, Eigen::ColMajor>::InnerIterator
中不需要使用 Eigen::ColMajor
,因为 that is the default initialization。所以,我删除了这条语句。
关于你的第二个问题,关于k
的迭代。
k
对两个 k=3,4
进行 迭代 ,但这些列中没有元素。因此,输出 k
的内部循环会调用 not。
如果我们在外循环和内循环中放置两个 k
声明性输出语句,这很容易看出。
例如
for(int k = 0; k < spmat.outerSize(); ++k) {
Rcpp::Rcout << "Overall k = " << k << std::endl << std::endl;
for(MappedSparseMatrix<double>::InnerIterator it(spmat,k); it; ++it) {
Rcpp::Rcout << "Inner k = " << k << std::endl;
}
}
避免using namespace std;
添加命名空间有时会产生意想不到的后果,尤其是像 std
.
这样大的后果
根据上面的观点并稍微简化您的示例,我们有以下简单的工作示例:
#include <RcppEigen.h>
// [[Rcpp::depends(RcppEigen)]]
using Eigen::MappedSparseMatrix;
using Eigen::SparseMatrix;
using Eigen::VectorXi;
using Eigen::Map;
// [[Rcpp::export]]
void createRec(Rcpp::RObject sparse_mat) {
MappedSparseMatrix<double> spmat(Rcpp::as<MappedSparseMatrix<double> >(sparse_mat));
long int nrow = spmat.rows();
Rcpp::NumericVector sim(nrow);
for(int k = 0; k < spmat.outerSize(); ++k) {
Rcpp::Rcout << "Overall k = " << k << std::endl << std::endl;
for(MappedSparseMatrix<double>::InnerIterator it(spmat,k); it; ++it) {
Rcpp::Rcout << "Inner k = " << k << std::endl
<< "value = " << it.value() << std::endl
<< "it.row = " << it.row() << std::endl
<< "it.col = " << it.col() << std::endl
<< "index = " << it.index() << std::endl;
}
}
}
/***R
# Setup values
id_row = c(2, 2, 4, 5)
id_col = c(1, 2, 3, 2)
vals = rep(1,4)
# Make the matrix
x = sparseMatrix(id_row, id_col, x = vals, dims = c(5, 5))
# Test the function
createRec(x)
*/
输出:
Overall k = 0
Inner k = 0
value = 1
it.row = 1
it.col = 0
index = 1
Overall k = 1
Inner k = 1
value = 1
it.row = 1
it.col = 1
index = 1
Inner k = 1
value = 1
it.row = 4
it.col = 1
index = 4
Overall k = 2
Inner k = 2
value = 1
it.row = 3
it.col = 2
index = 3
Overall k = 3
Overall k = 4
有关 Eigen 和 Rcpp 中稀疏矩阵的更多详细信息,您不妨阅读 Soren Hojsgaard 和 Doug Bates 的 Rcpp Gallery: Using iterators for sparse vectors and matrices。
我正在将一个小的稀疏矩阵(用于测试)传递给 R 中的 C++ 函数。该矩阵属于 class dgCMatrix,如下所示:
5 x 5 sparse Matrix of class "dgCMatrix"
[1,] . . . . .
[2,] 1 1 . . .
[3,] . . . . .
[4,] . . 1 . .
[5,] . 1 . . .
我正在按照文档 here 中的说明迭代此矩阵。 我的函数打印出迭代器的值和行索引、列索引。
C++函数定义如下:
#include <RcppEigen.h>
// [[Rcpp::depends(RcppEigen)]]
using Eigen::MappedSparseMatrix;
using Eigen::SparseMatrix;
using Eigen::VectorXi;
using Eigen::Map;
using namespace Rcpp;
using namespace std;
// [[Rcpp::export]]
void createRec(RObject sparse_mat, IntegerVector sparse_vec) {
const MappedSparseMatrix<int> spmat(as<MappedSparseMatrix<int> >(sparse_mat));
long int nrow = spmat.rows();
long int ncol = spmat.cols();
NumericVector sim(nrow);
for(int k=0;k<spmat.outerSize();k++){
for(SparseMatrix<int,Eigen::ColMajor>::InnerIterator it(spmat,k);it;++it){
cout<<"k="<<k<<endl;
cout<<"value="<<it.value()<<endl;
cout<<"it.row="<<it.row()<<endl;
cout<<"it.col="<<it.col()<<endl;
cout<<"index="<<it.index()<<endl;
}
}
}
对于上面给出的矩阵,打印了以下结果:
k=0
value=156148016
it.row=66211520
it.col=0
index=66211520
k=1
value=0
it.row=0
it.col=1
index=0
k=1
value=1
it.row=4
it.col=1
index=4
k=2
value=1
it.row=3
it.col=2
index=3
1.) k=0对应的值有什么解释吗?这些可能是由于以错误的方式传递矩阵造成的吗?
2.) k 在 outerSize 上迭代,它等于 5,为什么它不迭代 k=3,4?考虑到它是一个稀疏矩阵,迭代器会出现这种行为。
每当您看到像 156148016
或 66211520
这样的非常大的数字时,很可能是您遇到了未定义的行为 (UB) 或某个值未正确初始化。在这种情况下,它是后者。具体来说,dgCMatrix
class' 的基础类型是 double
而不是 int
.
The
dgCMatrix
class is a class of sparse numeric matrices in the compressed, sparse, column-oriented format. In this implementation the non-zero elements in the columns are sorted into increasing row order.dgCMatrix
is the "standard" class for sparse numeric matrices in theMatrix
package.
因此,当您尝试创建到底层 RObject
内存位置的映射时,需要一个额外的步骤来在请求的 不同[=88] 中重新创建对象=] 类型。添加 const
项后,我敢打赌这些条目会按预期进行,因为编译器可能会在内存中保留中间对象。
所以,变化如下:
MappedSparseMatrix<int> spmat(as<MappedSparseMatrix<int> >(sparse_mat));
至:
MappedSparseMatrix<double> spmat(as<MappedSparseMatrix<double> >(sparse_mat));
应该足够了。
linked example使用了SparseMatrix
矩阵,这里你使用的是MappedSparseMatrix
但没有为第二个循环设置合适的MappedSparseMatrix::InnerIterator
。
因此,我们有:
for(SparseMatrix<int,Eigen::ColMajor>::InnerIterator it(spmat,k);it;++it){
要去:
for(MappedSparseMatrix<double>::InnerIterator it(spmat,k);it;++it){
另外,请注意 SparseMatrix<int, Eigen::ColMajor>::InnerIterator
中不需要使用 Eigen::ColMajor
,因为 that is the default initialization。所以,我删除了这条语句。
关于你的第二个问题,关于k
的迭代。
k
对两个 k=3,4
进行 迭代 ,但这些列中没有元素。因此,输出 k
的内部循环会调用 not。
如果我们在外循环和内循环中放置两个 k
声明性输出语句,这很容易看出。
例如
for(int k = 0; k < spmat.outerSize(); ++k) {
Rcpp::Rcout << "Overall k = " << k << std::endl << std::endl;
for(MappedSparseMatrix<double>::InnerIterator it(spmat,k); it; ++it) {
Rcpp::Rcout << "Inner k = " << k << std::endl;
}
}
避免using namespace std;
添加命名空间有时会产生意想不到的后果,尤其是像 std
.
根据上面的观点并稍微简化您的示例,我们有以下简单的工作示例:
#include <RcppEigen.h>
// [[Rcpp::depends(RcppEigen)]]
using Eigen::MappedSparseMatrix;
using Eigen::SparseMatrix;
using Eigen::VectorXi;
using Eigen::Map;
// [[Rcpp::export]]
void createRec(Rcpp::RObject sparse_mat) {
MappedSparseMatrix<double> spmat(Rcpp::as<MappedSparseMatrix<double> >(sparse_mat));
long int nrow = spmat.rows();
Rcpp::NumericVector sim(nrow);
for(int k = 0; k < spmat.outerSize(); ++k) {
Rcpp::Rcout << "Overall k = " << k << std::endl << std::endl;
for(MappedSparseMatrix<double>::InnerIterator it(spmat,k); it; ++it) {
Rcpp::Rcout << "Inner k = " << k << std::endl
<< "value = " << it.value() << std::endl
<< "it.row = " << it.row() << std::endl
<< "it.col = " << it.col() << std::endl
<< "index = " << it.index() << std::endl;
}
}
}
/***R
# Setup values
id_row = c(2, 2, 4, 5)
id_col = c(1, 2, 3, 2)
vals = rep(1,4)
# Make the matrix
x = sparseMatrix(id_row, id_col, x = vals, dims = c(5, 5))
# Test the function
createRec(x)
*/
输出:
Overall k = 0
Inner k = 0
value = 1
it.row = 1
it.col = 0
index = 1
Overall k = 1
Inner k = 1
value = 1
it.row = 1
it.col = 1
index = 1
Inner k = 1
value = 1
it.row = 4
it.col = 1
index = 4
Overall k = 2
Inner k = 2
value = 1
it.row = 3
it.col = 2
index = 3
Overall k = 3
Overall k = 4
有关 Eigen 和 Rcpp 中稀疏矩阵的更多详细信息,您不妨阅读 Soren Hojsgaard 和 Doug Bates 的 Rcpp Gallery: Using iterators for sparse vectors and matrices。