12

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

David Alexander
  • 358
  • 1
  • 12
  • 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 Aug 31 '11 at 14:48
  • 3
    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 Aug 31 '11 at 14:58
  • 1
    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 Aug 31 '11 at 15:10
  • What kind of computations do you intend to do with your matrix? – Nicolas Grebille Aug 31 '11 at 16:34
  • @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 Aug 31 '11 at 17:27
  • I am doing explicit addition in the log domain. 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." – David Alexander Aug 31 '11 at 17:41
  • What you are asking for is a bit of a dependency inversion. I would expect the best practice would be to have the algorithm handle the handling of these values. – Tom Kerr Aug 31 '11 at 18:26
  • You say your matrix is banded. If so use a banded matrix. If you want to use `-Inf` as a sentinel value then you can just do so. – David Heffernan Aug 31 '11 at 19:17
  • Since you are in c++, why not use IEEE754's `-std::numeric_limits::infinity()` (or `double`)? You get `std::isinf` as bonus. – eudoxos Sep 06 '11 at 15:18

2 Answers2

2

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";
}
ravenspoint
  • 19,093
  • 6
  • 57
  • 103
  • Thanks ravenspoint, that's definitely a viable solution. My only fault with it is that it always involves a linear scan. If the underlying boost compressed_matrix has some fast caching scheme to speed up in-order accesses, we cannot take advantage of it in this way. But I do not know whether boost or eigen actually does implement such a scheme---I will have to check. – David Alexander Aug 31 '11 at 21:52
  • Boost compressed_matrix actually has a nice separation of functionality here: find_element locates the (i,j) element and returns a pointer, while operator() either dereferences that pointer or returns zero_. So I think the hook I am looking for is either to set zero_ somehow or inherit and override operator() – David Alexander Aug 31 '11 at 22:24
  • You may want to test the performance - that call to sort made by find_element looks expensive! – ravenspoint Aug 31 '11 at 22:44
  • According to [http://www.guwi17.de/ublas/matrix_sparse_usage.html](http://www.guwi17.de/ublas/matrix_sparse_usage.html), the appropriate way to fill a sparse_matrix "in order" is using push_back(i,j, val), which is then a constant-time operation. – David Alexander Sep 01 '11 at 15:39
  • Sure. But what has that got to do with extracting a value from a random location? – ravenspoint Sep 01 '11 at 19:23
  • The idea is that you are filling in each row up to a point where the entries start getting really small (very negative), and then go to the next row. Since the entries that are skipped are assumed to be really small, there is not much loss in having them be -FLT_MAX (since e^-FLT_MAX won't be that much different from e^(really small)). Thus we can then query any random (i,j) position and get a reasonable answer. – David Alexander Sep 03 '11 at 14:27
1

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.

David Alexander
  • 358
  • 1
  • 12
  • Have you tested this? It seems to me that the compressed matrix will return a zero value for absent values, which is exactly what you do not want. – ravenspoint Sep 01 '11 at 19:26
  • I tested it and it works. I returns value_type() as "zero", which in this case will be lfloat()=-FLT_MAX – David Alexander Sep 01 '11 at 21:15
  • That is neat. May I suggest a small optimization. Since each return of an absent value invokes the default constructor, you can speed it up a little by writing lfoat::lfloat() : value( -FLT_MAX ) {} – ravenspoint Sep 02 '11 at 13:16