Can anyone tell me the most effective way to extract the diagonal of a sparse dynamically allocated matrix and store it in a vector? I have tried using the diagonal() command, but I can't seem to make it work. Thanks for the help.

Cheers

Which Eigen version? VectorXd d = A.diagonal() with A a SparseMatrix<double> is working for me.

Is it possible to extract the diagonal of a SparseMatrix as a SparseVector rather than a full vector? Something like:

Code: Select all

SparseVector<double> d = A.diagonal()


Though currently that doesn't seem to work.

I can do

Code: Select all

SparseVector<double> d = A.diagonal().sparseView();


but I guess this is first constructing a dense vector and then pruning it, which will be unnecessarily costly if A is enormous and sparse along the diagonal.

Code: Select all

A.diagonal().sparseView();


is actually a good option as 1) A.diagonal() returns a lightweight expression and 2) I cannot think about a faster implementation than what this version is doing internally.

Thanks for the replies. I have a further question regarding this topic. What I really need to do, is as follows. I have a square, sparse matrix A as well as sparse matrices B and C. I want to form the following product:

B * D * C

where D is the sparse diagonal matrix with the diagonal of A. I want to do this as efficiently as possible. Any suggestions on how to form this product would be greatly appreciated. One approach I tried was

B * VectorXd(A.diagonal()).asDiagonal() * C,

however, (please correct me if I am wrong) I believe VectorXd(A.diagonal()).asDiagonal() forms a dense diagonal matrix. Thanks again for the help!

No need for the VectorXd temporary. (A.diagonal()).asDiagonal() is a Diagonal matrix for which products are optimized, regardless wether A is sparse or dense.

Thanks again for the help. However, when I try to do something like

SparseMatrix<double, RowMajor> S = B * ( A.diagonal( ) ).asDiagonal( ) * C;

I get the following compiler error from visual studio 2010:

c:\program files\eigen_3_20\eigen\src\sparsecore\sparsecwisebinaryop.h(60): error C2338: THE_STORAGE_ORDER_OF_BOTH_SIDES_MUST_MATCH
1> c:\program files\eigen_3_20\eigen\src\sparsecore\sparsecwisebinaryop.h(56) : while compiling class template member function 'Eigen::CwiseBinaryOpImpl<BinaryOp,Lhs,Rhs,StorageKind>::CwiseBinaryOpImpl(void)'
1> with
1> [
1> BinaryOp=Eigen::internal::scalar_product_op<double>,
1> Lhs=const Eigen::Block<const Eigen::SparseMatrix<double,1>,1,-1,true>,
1> Rhs=const Eigen::Transpose<const Eigen::Diagonal<const Eigen::SparseMatrix<double,1>>>,
1> StorageKind=Eigen::Sparse
1> ]

Works for me with this complete example (g++):

Code: Select all

#include <Eigen/Sparse>
using namespace Eigen;
int main() {
  SparseMatrix<double> A, B, C;

  SparseMatrix<double, RowMajor> S = B * ( A.diagonal( ) ).asDiagonal( ) * C;
}


What about other operations with the diagonal as a sparse diagonal matrix? For example what if I want to subtract the diagonal from a matrix?

I'd like to compute:

Code: Select all

A = (A - A.diagonal().asDiagonal()).eval();

Instead I do:

Code: Select all

SparseMatrix<double> I(A.rows(),A.cols());
A = (A - I*A.diagonal().asDiagonal()).eval();


But is there a better way to mix `asDiagonal()` with addition and subtraction? (or even just a better way to remove the diagonal from a sparse matrix?).

I did not check, but A.diagonal() should be writable for sparse matrices too. Therefore, it is simpler to do A.diagonal() = ...; For your particular example it is much faster to use the prune method:

Code: Select all

struct keep_diag {
   inline bool operator() (const Index& row, const Index& col, const Scalar&) const
    { return row!=col;  }
};
A.prune(keep_diag());


unless you want to keep the explicit zeros on the diagonal.

Should it be { return row==col; } instead.

Also compiler complains about type Scalar/Index, so I replace them with double/int.

Yes, it's rather return row==col, and Index/Scalar were place-holders for the actual types.