[Giuseppe Passino] Optimised maximum-residual BP...
authorJoris Mooij <joris.mooij@tuebingen.mpg.de>
Tue, 3 Mar 2009 08:02:48 +0000 (09:02 +0100)
committerJoris Mooij <joris.mooij@tuebingen.mpg.de>
Tue, 3 Mar 2009 08:02:48 +0000 (09:02 +0100)
...by using a reversed ordered set instead of the linear search (which can
yield enormous speedups - a factor 500 has been measured on a binary Ising grid
with 128x128 variables!)

ChangeLog
include/dai/bp.h
src/bp.cpp

index e2b3032..ad1152b 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,14 @@
+* [Giuseppe Passino] Optimised maximum-residual BP by using a reversed ordered
+  set instead of the linear search (which can yield enormous speedups - a factor
+  500 has been measured on a binary Ising grid with 128x128 variables!)
+* Added debug assertions to Var which check for inconsistent dimensions of 
+  variables with the same labels
+* [Giuseppe Passino] Added prefix ++ operator to State (State::operator++())
+* [Giuseppe Passino] Added iterators to FactorGraph (FactorGraph::begin, FactorGraph::end)
+* [Giuseppe Passino] Added iterators to TFactor (TFactor::begin, TFactor::end)
+* [Giuseppe Passino] Added iterators to TProb (TProb::begin, TProb::end)
+* [Giuseppe Passino] Added BP::findMaximum(), which can be used after running
+  max-product BP to construct a global state with maximum probability
 * Improved Makefile (libDAI now also builds out-of-the-box on MacOSX, 
   thanks to Dan Preston; merged Makefile.win and Makefile.shared into Makefile)
 * Fixed bug in calcMarginal, calcPairBeliefs, calcPairBeliefsNew where
index 5ee3b35..361e042 100644 (file)
@@ -40,11 +40,10 @@ namespace dai {
 
 
 /// Approximate inference algorithm "(Loopy) Belief Propagation"
-/** \todo Optimize BP_SEQMAX (it should use a better data structure than a vector for the residuals).
- */
 class BP : public DAIAlgFG {
     private:
         typedef std::vector<size_t> ind_t;
+           typedef std::multimap<double, std::pair<std::size_t, std::size_t> > LutType;
         struct EdgeProp {
             ind_t  index;
             Prob   message;
@@ -52,6 +51,8 @@ class BP : public DAIAlgFG {
             double residual;
         };
         std::vector<std::vector<EdgeProp> > _edges;
+        std::vector<std::vector<LutType::iterator> > _edge2lut;
+        LutType _lut;
         /// Maximum difference encountered so far
         double _maxdiff;
         /// Number of iterations needed
@@ -93,16 +94,22 @@ class BP : public DAIAlgFG {
 
     public:
         /// Default constructor
-        BP() : DAIAlgFG(), _edges(), _maxdiff(0.0), _iters(0U), props() {}
+        BP() : DAIAlgFG(), _edges(), _edge2lut(), _lut(), _maxdiff(0.0), _iters(0U), props() {}
 
         /// Copy constructor
-        BP( const BP &x ) : DAIAlgFG(x), _edges(x._edges), _maxdiff(x._maxdiff), _iters(x._iters), props(x.props) {}
+        BP( const BP &x ) : DAIAlgFG(x), _edges(x._edges), _edge2lut(x._edge2lut), _lut(x._lut), _maxdiff(x._maxdiff), _iters(x._iters), props(x.props) {
+            for( LutType::iterator l = _lut.begin(); l != _lut.end(); ++l )
+                _edge2lut[l->second.first][l->second.second] = l;
+        }
 
         /// Assignment operator
         BP& operator=( const BP &x ) {
             if( this != &x ) {
                 DAIAlgFG::operator=( x );
                 _edges = x._edges;
+                _lut = x._lut;
+                for( LutType::iterator l = _lut.begin(); l != _lut.end(); ++l )
+                    _edge2lut[l->second.first][l->second.second] = l;
                 _maxdiff = x._maxdiff;
                 _iters = x._iters;
                 props = x.props;
@@ -156,15 +163,8 @@ class BP : public DAIAlgFG {
         const double & residual(size_t i, size_t _I) const { return _edges[i][_I].residual; }
 
         void calcNewMessage( size_t i, size_t _I );
-        void updateMessage( size_t i, size_t _I ) {
-            if( props.damping == 0.0 ) {
-                message(i,_I) = newMessage(i,_I);
-                residual(i,_I) = 0.0;
-            } else {
-                message(i,_I) = (message(i,_I) ^ props.damping) * (newMessage(i,_I) ^ (1.0 - props.damping));
-                residual(i,_I) = dist( newMessage(i,_I), message(i,_I), Prob::DISTLINF );
-            }
-        }
+        void updateMessage( size_t i, size_t _I );
+        void updateResidual( size_t i, size_t _I, double r );
         void findMaxResidual( size_t &i, size_t &_I );
         /// Calculates unnormalized belief of variable
         void calcBeliefV( size_t i, Prob &p ) const;
index 4581cfa..fb30e9a 100644 (file)
@@ -101,9 +101,16 @@ void BP::construct() {
     // create edge properties
     _edges.clear();
     _edges.reserve( nrVars() );
+    _edge2lut.clear();
+    if( props.updates == Properties::UpdateType::SEQMAX )
+        _edge2lut.reserve( nrVars() );
     for( size_t i = 0; i < nrVars(); ++i ) {
         _edges.push_back( vector<EdgeProp>() );
         _edges[i].reserve( nbV(i).size() ); 
+        if( props.updates == Properties::UpdateType::SEQMAX ) {
+            _edge2lut.push_back( vector<LutType::iterator>() );
+            _edge2lut[i].reserve( nbV(i).size() );
+        }
         foreach( const Neighbor &I, nbV(i) ) {
             EdgeProp newEP;
             newEP.message = Prob( var(i).states() );
@@ -117,6 +124,8 @@ void BP::construct() {
 
             newEP.residual = 0.0;
             _edges[i].push_back( newEP );
+            if( props.updates == Properties::UpdateType::SEQMAX )
+                _edge2lut[i].push_back( _lut.insert( std::make_pair( newEP.residual, std::make_pair( i, _edges[i].size() - 1 ))) );
         }
     }
 }
@@ -128,12 +137,15 @@ void BP::init() {
         foreach( const Neighbor &I, nbV(i) ) {
             message( i, I.iter ).fill( c );
             newMessage( i, I.iter ).fill( c );
+            if( props.updates == Properties::UpdateType::SEQMAX )
+                               updateResidual( i, I.iter, 0.0 );
         }
     }
 }
 
 
 void BP::findMaxResidual( size_t &i, size_t &_I ) {
+/*
     i = 0;
     _I = 0;
     double maxres = residual( i, _I );
@@ -144,6 +156,12 @@ void BP::findMaxResidual( size_t &i, size_t &_I ) {
                 _I = I.iter;
                 maxres = residual( i, _I );
             }
+*/
+    assert( !_lut.empty() );
+    LutType::const_iterator largestEl = _lut.end();
+    --largestEl;
+    i  = largestEl->second.first;
+    _I = largestEl->second.second;
 }
 
 
@@ -217,6 +235,10 @@ void BP::calcNewMessage( size_t i, size_t _I ) {
         else
             newMessage(i,_I) = marg;
     }
+
+    // Update the residual if necessary
+    if( props.updates == Properties::UpdateType::SEQMAX )
+        updateResidual( i, _I , dist( newMessage( i, _I ), message( i, _I ), Prob::DISTLINF ) );
 }
 
 
@@ -245,8 +267,6 @@ double BP::run() {
         for( size_t i = 0; i < nrVars(); ++i )
             foreach( const Neighbor &I, nbV(i) ) {
                 calcNewMessage( i, I.iter );
-                // calculate initial residuals
-                residual( i, I.iter ) = dist( newMessage( i, I.iter ), message( i, I.iter ), Prob::DISTLINF );
             }
     } else {
         update_seq.reserve( nredges );
@@ -273,10 +293,8 @@ double BP::run() {
                     if( J.iter != _I ) {
                         foreach( const Neighbor &j, nbF(J) ) {
                             size_t _J = j.dual;
-                            if( j != i ) {
+                            if( j != i )
                                 calcNewMessage( j, _J );
-                                residual( j, _J ) = dist( newMessage( j, _J ), message( j, _J ), Prob::DISTLINF );
-                            }
                         }
                     }
                 }
@@ -465,12 +483,40 @@ string BP::identify() const {
 void BP::init( const VarSet &ns ) {
     for( VarSet::const_iterator n = ns.begin(); n != ns.end(); ++n ) {
         size_t ni = findVar( *n );
-        foreach( const Neighbor &I, nbV( ni ) )
-            message( ni, I.iter ).fill( props.logdomain ? 0.0 : 1.0 );
+        foreach( const Neighbor &I, nbV( ni ) ) {
+            double val = props.logdomain ? 0.0 : 1.0;
+            message( ni, I.iter ).fill( val );
+            newMessage( ni, I.iter ).fill( val );
+            if( props.updates == Properties::UpdateType::SEQMAX )
+                updateResidual( ni, I.iter, 0.0 );
+        }
     }
 }
 
 
+void BP::updateMessage( size_t i, size_t _I ) {
+    if( props.damping == 0.0 ) {
+        message(i,_I) = newMessage(i,_I);
+        if( props.updates == Properties::UpdateType::SEQMAX )
+            updateResidual( i, _I, 0.0 );
+    } else {
+        message(i,_I) = (message(i,_I) ^ props.damping) * (newMessage(i,_I) ^ (1.0 - props.damping));
+        if( props.updates == Properties::UpdateType::SEQMAX )
+            updateResidual( i, _I, dist( newMessage(i,_I), message(i,_I), Prob::DISTLINF ) );
+    }
+}
+
+
+void BP::updateResidual( size_t i, size_t _I, double r ) {
+       EdgeProp* pEdge = &_edges[i][_I];
+       pEdge->residual = r;
+       
+       // rearrange look-up table (delete and reinsert new key)
+       _lut.erase( _edge2lut[i][_I] );
+       _edge2lut[i][_I] = _lut.insert( std::make_pair( r, std::make_pair(i, _I) ) );
+}
+
+
 std::vector<size_t> BP::findMaximum() const {
     std::vector<size_t> maximum( nrVars() );
     std::vector<bool> visitedVars( nrVars(), false );