CODEWITHSUNDEEP
c++matlabmatrixeigen
Adrian
Adrian

Reputation: 1188

How to calculate diff of sparse matrix in C++ using Eigen

I'm attempting to implement the below MATLAB code (from Eilers, 2003) in C++ using the Eigen library.

m = length(y);
E = speye(m);
D = diff(E, d);
C = chol(E + lambda * D' * D);
z = C\(C'\y);

I can create a sparse matrix easily enough:

Eigen::SparseMatrix<int> E(m,m);
E.setIdentity();

But what's the equivalent of diff in Eigen for a sparse matrix?

I could alternatively assume an order of differences value of 3, giving the coefficients of -1, 3, -3, 1. For an initial matrix of size 10 x 10, this would look like: 

-1 3 -3 1 0 0 0 0 0 0
0 -1 3 -3 1 0 0 0 0 0 
0 0 -1 3 -3 1 0 0 0 0
0 0 0 -1 3 -3 1 0 0 0
0 0 0 0 -1 3 -3 1 0 0
0 0 0 0 0 -1 3 -3 1 0
0 0 0 0 0 0 -1 3 -3 1

If I can't use diff then how can I efficiently generate the above sparse matrix?

Upvotes: 3

Views: 1675

Answers (1)

Tim B
Tim B

Reputation: 3143

MATLAB's diff simply finds the difference between adjacent elements in the matrix (defaulting to row-wise I think), so an easy way to achieve this is Eigen is through simple matrix subtraction, you just need to make sure you generate the right matrices.

You need a function that, given an m x n matrix E, will create two matrices E1 and E2. E1 is simply the first m-1 rows of E, E2 the last m-1 rows (E without the first row). Subtracting E1 from E2 should then give you what you want.

Something like this:

Eigen::SparseMatrix<double> diff(Eigen::SparseMatrix<double> E) {
    Eigen::SparseMatrix<double> E1 = E.block(0, 0, E.rows()-1, E.cols());
    Eigen::SparseMatrix<double> E2 = E.block(1, 0, E.rows()-1, E.cols());
    return E2 - E1;
}

Although I've not tested this, the idea should be correct. You can apply it recursively to get the nth difference approximation.

This operates row-wise, if you want the column wise diff it should be trivial to adapt this method.

Upvotes: 3

Related Questions