Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Sparse matrix class with parameterizable "zero"

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."

like image 884
David Alexander Avatar asked Aug 31 '11 14:08

David Alexander


2 Answers

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";
}
like image 117
ravenspoint Avatar answered Sep 28 '22 12:09

ravenspoint


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.

like image 36
David Alexander Avatar answered Sep 28 '22 12:09

David Alexander