So ... from the example and tests it looks like you nailed this on the head? Read to merge now?

I think it is peculiar that Armadillo allows for the column pointers to be decreasing in this case. You may want to check with the author if that is intentional.

It is quite common to have code that checks the last element of the column pointers (position n + 1 in R and index n in C++) to determine the number of non-zeros in the sparse matrix.

There is some slight difference in column pointers between Armadillo and Matrix pkg in R.

In a matrix below:

3x5sparseMatrixofclass"dgCMatrix"              [1,]1....[2,].1...[3,]..1..

Thep slot is0 1 2 3 3 3 and thecol_ptrs in Armadillo is0 1 2 3 0 0.

I think this can be intentional since we only have three nonzero elements.

I will double check with the author.

Any news on checking with Conrad or other experts on sparse matrices?

The sparse matrix in armadillo is by Ryan (our MLPACK hero) and he is on travel.

@thirdwing Where are we on this and the off-list reply from Ryan? What's next to move this forward?

I am afraid we didn't get any reply from Ryan.

So what is next? Is it conceivable that Armadillo and Matrix have different ideas about corner cases of the DSC (or other) formats? Shall we document that / create test cases?

Every description of the compressed sparse row/column format that I have read states that the non-zeros for row/column j are in positionsp[j] top[j+1]-1 of the pointer array. SeeTim Davis's book orWikipedia

Another way of saying this is that, in 0-based indices the pointer array is recursively defined as

p[0] = 0p[j] = p[j-1] + number of non-zeros in column j-1

It is common to determine the number of non-zeros in each column as the successive differences in the column pointers, which is whyp is usually defined to be non-decreasing. As I mentioned, this leads to the last element of the column pointer vector being the number of non-zeros in the matrix.

To me the solution is to convince Conrad to modify the Armadillo representation so that it is in keeping with other sparse matrix representations.

Thanks for chiming in@dmbates. I am trying to find a private email thread where Ryan described the scheme to (to KK, who forwarded to me). I think it said the same.

@rcurtin Do you mind if I paste your definition in here?

Sure, no problem from my end.

From an earlier email by@rcurtin :

An answer to your two questions:

For a matrix generated by "arma::speye(3, 5)",

(1) Is the length of the col_ptrs always n_cols + 1? Even there are
only three nonzero elements?

Yes, this is the 'compressed sparse column' format. Technically in this
implementation, the length is always n_cols + 2.

The column pointers hold the index of the first element in each column.

So, col_ptrs[0] is always 0 (note below), col_ptrs[1] is the index of
the first nonzero element in column 1, and so forth. col_ptrs[n_cols]
is always n_nonzero.

col_ptrs[n_cols + 1] is a "sentinel" value used to terminate loops, and
should always be numeric_limits::max().
> Note: since col_ptrs[0] is always 0 and col_ptrs[n_cols] is always
n_nonzero, there's some very minor optimization that could be done. But
I've never considered it to be worthwhile.

(2) If the length is n_cols + 1, the content of col_ptrs is 0 1 2 3 0
0 . The column pointers are padded by zeros since three are only three
nonzero elements.

I'm not sure I understand what you mean, are you saying that you run
speye() and the resulting matrix has column pointers [0 1 2 3 0 0]?
I'm traveling right now, so any response to this will unfortunately take
a couple of days---but I will get back to it. :)

and by now we already have some new code to test thanks to@conradsnicta.

I think this has been fixed byconradsnicta/armadillo-code@763ed1a

That would be awesome. I haven't had a chance to look into the diff but I may then apply this to the pending RcppArmadillo release based on 7.960.0 (as well as all of the changes by@binxiangni).

I have applied that patch locally.@eddelbuettel

You mean you added it to this (pending) PR? That's very helpful.

And I presume the unit test fails without it, but now passes?

I add that patch in this PR and it just fixed the problem and passed new tests.