Backward Analysis using Preimage Approximation

CHANGELOG.md

@@ -5,6 +5,7 @@
   - Dropped support for parsing Tensorflow network format. Newest Marabou version that supports Tensorflow is at commit 190555573e4702.
   - Fixed bug in the parsing of `transpose` nodes in command line C++ parser.
   - Implemented forward-backward abstract interpretation, symbolic bound tightening, interval arithmetic and simulations for all activation functions.
+  - Implemented backward analysis using preimage-approximation algorithm for `Relu`, `LeakyRelu`, `Sign` and `Bilinear` Layers.
   - Added the BaBSR heuristic as a new branching strategy for ReLU Splitting
 
 ## Version 2.0.0

src/configuration/GlobalConfiguration.cpp

@@ -67,6 +67,13 @@ const unsigned GlobalConfiguration::ROW_BOUND_TIGHTENER_SATURATION_ITERATIONS =
 const double GlobalConfiguration::COST_FUNCTION_ERROR_THRESHOLD = 0.0000000001;
 
 const unsigned GlobalConfiguration::SIMULATION_RANDOM_SEED = 1;
+const unsigned GlobalConfiguration::VOLUME_ESTIMATION_RANDOM_SEED = 1;
+const unsigned GlobalConfiguration::PREIMAGE_APPROXIMATION_OPTIMIZATION_RANDOM_SEED = 1;
+const unsigned GlobalConfiguration::VOLUME_ESTIMATION_ITERATIONS = 25000;
+const unsigned GlobalConfiguration::PREIMAGE_APPROXIMATION_OPTIMIZATION_MAX_ITERATIONS = 25;
+const double GlobalConfiguration::PREIMAGE_APPROXIMATION_OPTIMIZATION_STEP_SIZE = 0.025;
+const double GlobalConfiguration::PREIMAGE_APPROXIMATION_OPTIMIZATION_LEARNING_RATE = 0.25;
+const double GlobalConfiguration::PREIMAGE_APPROXIMATION_OPTIMIZATION_WEIGHT_DECAY = 0;
 
 const bool GlobalConfiguration::USE_HARRIS_RATIO_TEST = true;
 

src/configuration/GlobalConfiguration.h

@@ -151,6 +151,27 @@ class GlobalConfiguration
     // Random seed for generating simulation values.
     static const unsigned SIMULATION_RANDOM_SEED;
 
+    // Random seed for EstimateVolume procedure (PreimageApproximation).
+    static const unsigned VOLUME_ESTIMATION_RANDOM_SEED;
+
+    // Number of iterations for EstimateVolume procedure (PreimageApproximation).
+    static const unsigned VOLUME_ESTIMATION_ITERATIONS;
+
+    // Random seed for PreimageApproximation optimization.
+    static const unsigned PREIMAGE_APPROXIMATION_OPTIMIZATION_RANDOM_SEED;
+
+    // Maximum iterations for PreimageApproximation optimization.
+    static const unsigned PREIMAGE_APPROXIMATION_OPTIMIZATION_MAX_ITERATIONS;
+
+    // Step size for PreimageApproximation optimization.
+    static const double PREIMAGE_APPROXIMATION_OPTIMIZATION_STEP_SIZE;
+
+    // Learning rate for PreimageApproximation optimization.
+    static const double PREIMAGE_APPROXIMATION_OPTIMIZATION_LEARNING_RATE;
+
+    // Weight decay for PreimageApproximation optimization.
+    static const double PREIMAGE_APPROXIMATION_OPTIMIZATION_WEIGHT_DECAY;
+
     // How often should projected steepest edge reset the reference space?
     static const unsigned PSE_ITERATIONS_BEFORE_RESET;
 

src/configuration/OptionParser.cpp

@@ -267,7 +267,8 @@ void OptionParser::initialize()
                 &( ( *_stringOptions )[Options::MILP_SOLVER_BOUND_TIGHTENING_TYPE] ) )
                 ->default_value( ( *_stringOptions )[Options::MILP_SOLVER_BOUND_TIGHTENING_TYPE] ),
             "The MILP solver bound tightening type: "
-            "lp/backward-once/backward-converge/lp-inc/milp/milp-inc/iter-prop/none." )
+            "lp/backward-once/backward-converge/backward-preimage-approx/lp-inc/milp/milp-inc/"
+            "iter-prop/none." )
 #endif
         ;
 

src/configuration/Options.cpp

@@ -209,6 +209,8 @@ MILPSolverBoundTighteningType Options::getMILPSolverBoundTighteningType() const
             return MILPSolverBoundTighteningType::BACKWARD_ANALYSIS_ONCE;
         if ( strategyString == "backward-converge" )
             return MILPSolverBoundTighteningType::BACKWARD_ANALYSIS_CONVERGE;
+        if ( strategyString == "backward-preimage-approx" )
+            return MILPSolverBoundTighteningType::BACKWARD_ANALYSIS_PREIMAGE_APPROX;
         else if ( strategyString == "milp" )
             return MILPSolverBoundTighteningType::MILP_ENCODING;
         else if ( strategyString == "milp-inc" )

src/engine/Engine.cpp

@@ -1594,6 +1594,7 @@ void Engine::performMILPSolverBoundedTightening( Query *inputQuery )
         case MILPSolverBoundTighteningType::LP_RELAXATION_INCREMENTAL:
         case MILPSolverBoundTighteningType::BACKWARD_ANALYSIS_ONCE:
         case MILPSolverBoundTighteningType::BACKWARD_ANALYSIS_CONVERGE:
+        case MILPSolverBoundTighteningType::BACKWARD_ANALYSIS_PREIMAGE_APPROX:
             _networkLevelReasoner->lpRelaxationPropagation();
             break;
         case MILPSolverBoundTighteningType::MILP_ENCODING:
@@ -1658,6 +1659,15 @@ void Engine::performAdditionalBackwardAnalysisIfNeeded()
                 printf( "Backward analysis tightened %u bounds\n", tightened );
         }
     }
+
+    if ( _milpSolverBoundTighteningType ==
+         MILPSolverBoundTighteningType::BACKWARD_ANALYSIS_PREIMAGE_APPROX )
+    {
+        performMILPSolverBoundedTightening( &( *_preprocessedQuery ) );
+        unsigned tightened = performSymbolicBoundTightening( &( *_preprocessedQuery ) );
+        if ( _verbosity > 0 )
+            printf( "Backward analysis tightened %u bounds\n", tightened );
+    }
 }
 
 void Engine::performMILPSolverBoundedTighteningForSingleLayer( unsigned targetIndex )
@@ -1685,6 +1695,7 @@ void Engine::performMILPSolverBoundedTighteningForSingleLayer( unsigned targetIn
             return;
         case MILPSolverBoundTighteningType::BACKWARD_ANALYSIS_ONCE:
         case MILPSolverBoundTighteningType::BACKWARD_ANALYSIS_CONVERGE:
+        case MILPSolverBoundTighteningType::BACKWARD_ANALYSIS_PREIMAGE_APPROX:
         case MILPSolverBoundTighteningType::ITERATIVE_PROPAGATION:
         case MILPSolverBoundTighteningType::NONE:
             return;
@@ -3879,4 +3890,4 @@ void Engine::addPLCLemma( std::shared_ptr<PLCLemma> &explanation )
     ASSERT( explanation && _UNSATCertificate && _UNSATCertificateCurrentPointer )
     _statistics.incUnsignedAttribute( Statistics::NUM_LEMMAS );
     _UNSATCertificateCurrentPointer->get()->addPLCLemma( explanation );
-}
+}

src/engine/MILPSolverBoundTighteningType.h

@@ -34,6 +34,9 @@ enum class MILPSolverBoundTighteningType {
     // Perform backward analysis
     BACKWARD_ANALYSIS_ONCE = 5,
     BACKWARD_ANALYSIS_CONVERGE = 6,
+    // Perform backward analysis using the PreimageApproximation Algorithm (arXiv:2305.03686v4
+    // [cs.SE])
+    BACKWARD_ANALYSIS_PREIMAGE_APPROX = 7,
     // Option to have no MILP bound tightening performed
     NONE = 10,
 };

src/engine/PolygonalTightening.h

@@ -0,0 +1,101 @@
+/*********************                                                        */
+/*! \file PolygonalTightening.h
+ ** \verbatim
+ ** Top contributors (to current version):
+ **   Duligur Ibeling, Guy Katz, Ido Shmuel
+ ** This file is part of the Marabou project.
+ ** Copyright (c) 2017-2024 by the authors listed in the file AUTHORS
+ ** in the top-level source directory) and their institutional affiliations.
+ ** All rights reserved. See the file COPYING in the top-level source
+ ** directory for licensing information.\endverbatim
+ **
+ ** [[ Add lengthier description here ]]
+
+**/
+
+#ifndef __PolygonalTightening_h__
+#define __PolygonalTightening_h__
+
+#include "Map.h"
+#include "NeuronIndex.h"
+
+#include <cstdio>
+
+class PolygonalTightening
+{
+public:
+    enum PolygonalBoundType {
+        LB = 0,
+        UB = 1,
+    };
+
+    PolygonalTightening( Map<NLR::NeuronIndex, double> neuronToCoefficient,
+                         double value,
+                         PolygonalBoundType type )
+        : _neuronToCoefficient( neuronToCoefficient )
+        , _value( value )
+        , _type( type )
+    {
+    }
+
+    /*
+      The coefficient of each neuron.
+    */
+    Map<NLR::NeuronIndex, double> _neuronToCoefficient;
+
+    /*
+      Its new value.
+    */
+    double _value;
+
+    /*
+      Whether the tightening tightens the
+      lower bound or the upper bound.
+    */
+    PolygonalBoundType _type;
+
+    /*
+      Equality operator.
+    */
+    bool operator==( const PolygonalTightening &other ) const
+    {
+        bool allFound = true;
+        for ( const auto &pair : _neuronToCoefficient )
+        {
+            bool currentFound = false;
+            for ( const auto &otherPair : other._neuronToCoefficient )
+            {
+                currentFound |= ( pair.first._layer == otherPair.first._layer &&
+                                  pair.first._neuron == otherPair.first._neuron &&
+                                  pair.second == otherPair.second );
+            }
+            allFound &= currentFound;
+        }
+        bool result = allFound && _value == other._value && _type == other._type &&
+                      _neuronToCoefficient.size() == other._neuronToCoefficient.size();
+        return result;
+    }
+
+    void dump() const
+    {
+        printf( "PolygonalTightening: x %s %.2lf\n", _type == LB ? ">=" : "<=", _value );
+
+        for ( const auto &pair : _neuronToCoefficient )
+        {
+            printf( "NeuronIndex: (layer %u, neuron %u), Coefficient: %.2lf )\n",
+                    pair.first._layer,
+                    pair.first._neuron,
+                    pair.second );
+        }
+    }
+};
+
+#endif // __PolygonalTightening_h__
+
+//
+// Local Variables:
+// compile-command: "make -C ../.. "
+// tags-file-name: "../../TAGS"
+// c-basic-offset: 4
+// End:
+// s

src/engine/SumOfInfeasibilitiesManager.cpp

@@ -195,8 +195,19 @@ void SumOfInfeasibilitiesManager::proposePhasePatternUpdateRandomly()
     } );
 
     // First, pick a pl constraint whose cost component we will update.
-    unsigned index = (unsigned)T::rand() % _plConstraintsInCurrentPhasePattern.size();
-    PiecewiseLinearConstraint *plConstraintToUpdate = _plConstraintsInCurrentPhasePattern[index];
+    bool fixed = true;
+    PiecewiseLinearConstraint *plConstraintToUpdate;
+    while ( fixed )
+    {
+        if ( _plConstraintsInCurrentPhasePattern.empty() )
+            return;
+
+        unsigned index = (unsigned)T::rand() % _plConstraintsInCurrentPhasePattern.size();
+        plConstraintToUpdate = _plConstraintsInCurrentPhasePattern[index];
+        fixed = plConstraintToUpdate->phaseFixed();
+        if ( fixed )
+            removeCostComponentFromHeuristicCost( plConstraintToUpdate );
+    }
 
     // Next, pick an alternative phase.
     PhaseStatus currentPhase = _currentPhasePattern[plConstraintToUpdate];

src/nlr/CMakeLists.txt

@@ -19,6 +19,11 @@ network_level_reasoner_add_unit_test(NetworkLevelReasoner)
 network_level_reasoner_add_unit_test(WsLayerElimination)
 network_level_reasoner_add_unit_test(ParallelSolver)
 
+if (${ENABLE_GUROBI})
+    network_level_reasoner_add_unit_test(LPRelaxation)
+endif()
+
 if (${BUILD_PYTHON})
     target_include_directories(${MARABOU_PY} PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}")
 endif()
+