Matrix Diagonal Fill Bugged for Nonsymmetric #619

Closed
coatless opened this Issue Dec 27, 2016 · 3 comments

Projects

None yet

3 participants

@coatless
Contributor
coatless commented Dec 27, 2016 edited

Consider:

Rcpp::cppFunction("NumericMatrix fdiag(NumericMatrix x) {
    x.fill_diag(0.0);
    return x;
}")

The issue arises with the .fill_diag()'s fill_diag__dispatch() method in Matrix.h primarily with a nonsymmetric matrix under either n > p or n < p.

Case 1: Nonsymmetric (n < p)

m1 <- matrix(1:10, nrow = 2, ncol = 5, byrow = TRUE)
fdiag(m1)
     [,1] [,2] [,3] [,4] [,5]
[1,]    0    2    3    0    5
[2,]    6    7    8    9   10

Expected:

     [,1] [,2] [,3] [,4] [,5]
[1,]    0    2    3    4    5
[2,]    6    0    8    9   10

Case 2: Nonsymmetric (n > p)

m3 <- matrix(1:10, nrow = 5, ncol = 2, byrow = TRUE)
fdiag(m3)
     [,1] [,2]
[1,]    0    2
[2,]    3    4
[3,]    5    6
[4,]    0    8
[5,]    9   10

Expected:

     [,1] [,2]
[1,]    0    2
[2,]    3    0
[3,]    5    6
[4,]    7    8
[5,]    9   10

Case 3: n = p

Though, the square matrix iteration works correctly.

m3 <- matrix(1:9, nrow = 3, ncol = 3, byrow = TRUE)
fdiag(m3)
     [,1] [,2] [,3]
[1,]    0    2    3
[2,]    4    0    6
[3,]    7    8    0
@eddelbuettel
Member
eddelbuettel commented Dec 27, 2016 edited

Hm. Shall we warn, or even stop if we have a non-symmetric matrix? Maybe require an override var?
As in fill_diag(double val, bool nonsym=FALSE) ?

What is fill_diag() being modeled after?

Maybe still easiest to just refuse... Dunno.

Edit OTOH it is clearly named ..._diag() and if uses on something else ... you get to keep the pieces.

@nathan-russell
Contributor

I think it would be pretty safe to go with the main diagonal definition, i.e. diagonal elements are where i == j. This seems to be what R does anyhow:

`diag<-`(matrix(1, 3, 4), 0)
#      [,1] [,2] [,3] [,4]
# [1,]    0    1    1    1
# [2,]    1    0    1    1
# [3,]    1    1    0    1

`diag<-`(matrix(1, 4, 3), 0)
#      [,1] [,2] [,3]
# [1,]    0    1    1
# [2,]    1    0    1
# [3,]    1    1    0
# [4,]    1    1    1

`diag<-`(matrix(1, 3, 3), 0)
#      [,1] [,2] [,3]
# [1,]    0    1    1
# [2,]    1    0    1
# [3,]    1    1    0

I don't know why an iterator is used in the current implementation but this would be trivial with an ordinary loop:

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericMatrix FillDiag(NumericMatrix x) {
    R_xlen_t i = 0, n = x.nrow() < x.ncol() ? x.nrow() : x.ncol();
    for (; i < n; i++) {
        x(i, i) = 0.0;
    }
    return x;
}

/*** R

FillDiag(matrix(1, 3, 4))
#      [,1] [,2] [,3] [,4]
# [1,]    0    1    1    1
# [2,]    1    0    1    1
# [3,]    1    1    0    1

FillDiag(matrix(1, 4, 3))
#      [,1] [,2] [,3]
# [1,]    0    1    1
# [2,]    1    0    1
# [3,]    1    1    0
# [4,]    1    1    1

FillDiag(matrix(1, 3, 3))
#      [,1] [,2] [,3]
# [1,]    0    1    1
# [2,]    1    0    1
# [3,]    1    1    0

*/
@coatless
Contributor

Apologies for the lack of context regarding .fill_diag(). I was a bit surprised when I stumbled across the behavior as I was documenting the matrix features.

To quickly address the points:

  • The .fill_diag() is attempting to mimic diag(x) <- value and, thus, it does have a place.
  • There is no need to add in an extra parameter e.g. nonsym as the diag fill is based off of the major diagonal.
  • As @nathan-russell suggested, the main fix is most likely to move away from the iterator approach.
    • I'm wondering if the error is due to a UB in Cases 1 & 2 since the iterator is extended past the end. Though, given the pattern I have a feeling it is related to the elements being organized in the Vector by column and not row.
    • One last note, it looks as if this function also didn't receive the face lift over to using R_xlen_t from int.
  • Side note: Might be nice to extend the sugar diag() to support extracting the antidiagonal and adding an antidiagonal fill option to the fill_diag() method

Fun trivia factoid, the bug has been present since at least 0.11.0.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment