Sparse matrix class with parameterizable "zero"

12
votes

I am doing some computations on a sparse matrix of floats in the log domain, so the "empty" entries are actually -Inf (using -FLT_MAX). I'm using a custom sparse matrix class right now but I am eager to swap in an off-the-shelf replacement.

This is in C++. My inclinations were to look at the compressed column matrices in Eigen and Boost uBlas. However it is not clear that either supports a custom value for "zero" (perhaps provided by a template parameter). Does anyone have a suggestion?

Clarification:

What I want is this: for any cell (i,j) that has not been "set" previously, I would like mat[i,j] to return -Inf ... so this is perhaps better described as a "default" value for the "empty" entries of the sparse matrix.

I am using this to perform HMM recursions (Viterbi, sum-product) with probabilities kept in the log domain to avoid underflow.

I am not doing any matrix operations ... I am just filling in a dynamic programming tableau, essentially. I want to use a sparse matrix class because I am only filling in a band of the matrix and I would like efficient memory use. The compressed band matrices would give good performance since I am filling in the matrix "in order."

2
c++mathboostsparse-matrixeigen
The "zero" elements of a conventional sparse matrix are not stored, nor are they used in any numerical calculations. I understand that log(0) = -Inf but I don't understand what you're trying to do - if you add -FLT_MAX for all zero elements you'll end up with a dense matrix, and at any rate I can't see how doing numerical operations on matrices involving -Inf is meaningful. Maybe I've missed the point? - Darren Engwirda
For example, x * 0 == 0, but x * -Inf != -Inf when x is negative (and likewise x + 0 == x, but x + -Inf != x). So when doing a sparse matrix multiplication you can't just ignore the so-called empty values if those values are -Inf, whereas you can jump straight past a bunch of zeroes. - Steve Jessop
I think what the OP wants is a sparse matrix where the zero cells are in fact some other value and not zero. So instead of a matrix where most of the values are zero (so it's sparse), most of the values are some value X (but still sparse in respect to non-X values). - Skizz
What kind of computations do you intend to do with your matrix? - Nicolas Grebille
@David: the questioner did explain - the values in the matrix are logs of probabilities, and will be anti-logged before multiplication. -Inf is used as an approximation to log(0). So the sparse thing presumably doesn't needn't matrix arithmetic (which is meaningless in the log domain AFAIK), just a sparse 2-dimensional array with a tunable empty value would do. - Steve Jessop

2 Answers

2
votes

How about something like this?

class compressed_matrix_nonzero_default : public boost::numeric::ublas::compressed_matrix<double>
{
    double def;
public:
    compressed_matrix_nonzero_default( int s1, int s2 )
        : boost::numeric::ublas::compressed_matrix<double>(s1,s2)
        , def(0)
    {

    }
    void setDefault( double d ) { def = d; }
    double value( int i, int j )
    {
        typedef boost::numeric::ublas::compressed_matrix<double>::iterator1 it1_t;
        typedef boost::numeric::ublas::compressed_matrix<double>::iterator2 it2_t;
        for (it1_t it1 = begin1(); it1 != end1(); it1++)
        {
            if( it1.index1() <  i )
                continue;
            if( it1.index1() > i ) {
                return def;
            }
            for (it2_t it2 = it1.begin(); it2 != it1.end(); it2++)
            {
                if( it2.index2() < j )
                    continue;
                if( it2.index2() == j )
                    return *it2;
                if( it2.index2() > j )
                    return def;
            }


        }
        return def;
    }

};

Usage

compressed_matrix_nonzero_default MNZ(3,3);
MNZ.setDefault(-100);
MNZ (1,1) = 45;

for( int i = 0; i < 3; i++ ) {
    for( int j = 0; j < 3; j++ ) {
        std::cout << MNZ.value(i,j) << ",";
    }
    std::cout << "\n";
}
1
votes

The solution I have currently cooked up is this. Define a class lfloat:

class lfloat {
  float value;
public:
  lfloat(float f=-FLT_MAX)
  {
    value = f;
  }

  lfloat& operator=(float f)
  {
    value = f;
    return *this;
  }

  operator float()   { return value; }
};

and use it like so:

compressed_matrix<lfloat> M1(3,3);

This way we are not rewriting any of the functionality in the boost matrix classes, but we should get the desired result.