Properly implemented finding the smallest diagonal entry of a square …

…matrix

uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp

@@ -210,7 +210,8 @@ namespace uno {
       const double linearized_residual_reduction = current_iterate.progress.infeasibility - linearized_residual;
       const double lowest_linearized_residual_reduction = current_iterate.progress.infeasibility - residual_lowest_violation;
       if (lowest_linearized_residual_reduction < 0.) {
-         WARNING << "The lowest_linearized_residual_reduction quantity is negative. Negative curvature in your problem\n";
+         WARNING << "lowest_linearized_residual_reduction = " << lowest_linearized_residual_reduction <<
+         " is negative. Negative curvature in your problem\n";
       }
       return (linearized_residual_reduction >= this->parameters.epsilon1 * lowest_linearized_residual_reduction);
    }

uno/ingredients/hessian_models/ConvexifiedHessian.cpp

@@ -29,16 +29,16 @@ namespace uno {
       problem.evaluate_lagrangian_hessian(primal_variables, constraint_multipliers, this->hessian);
       this->evaluation_count++;
       // regularize (only on the original variables) to convexify the problem
-      DEBUG2 << "hessian before convexification: " << this->hessian;
       this->regularize(statistics, this->hessian, problem.get_number_original_variables());
    }
 
    // Nocedal and Wright, p51
    void ConvexifiedHessian::regularize(Statistics& statistics, SymmetricMatrix<size_t, double>& hessian, size_t number_original_variables) {
+      DEBUG << "Current Hessian:\n" << hessian << '\n';
       const double smallest_diagonal_entry = hessian.smallest_diagonal_entry(number_original_variables);
       DEBUG << "The minimal diagonal entry of the matrix is " << smallest_diagonal_entry << '\n';
 
-      double regularization_factor = (smallest_diagonal_entry <= 0.) ? this->regularization_initial_value - smallest_diagonal_entry : 0.;
+      double regularization_factor = (smallest_diagonal_entry > 0.) ? 0. : this->regularization_initial_value - smallest_diagonal_entry;
       bool good_inertia = false;
       bool symbolic_factorization_performed = false;
       while (not good_inertia) {
@@ -48,6 +48,8 @@ namespace uno {
                return (variable_index < number_original_variables) ? regularization_factor : 0.;
             });
          }
+         DEBUG << "Current Hessian:\n" << hessian << '\n';
+
          // perform the symbolic factorization only once
          if (not symbolic_factorization_performed) {
             this->linear_solver->do_symbolic_factorization(hessian);
@@ -68,4 +70,4 @@ namespace uno {
       }
       statistics.set("regulariz", regularization_factor);
    }
-} // namespace
+} // namespace

uno/linear_algebra/SymmetricMatrix.hpp

@@ -4,6 +4,7 @@
 #ifndef UNO_SYMMETRICMATRIX_H
 #define UNO_SYMMETRICMATRIX_H
 
+#include <algorithm>
 #include <memory>
 #include <functional>
 #include <cassert>
@@ -63,20 +64,19 @@ namespace uno {
    SymmetricMatrix<IndexType, ElementType>::SymmetricMatrix(size_t dimension, size_t capacity, bool use_regularization, const std::string& sparse_format) :
          sparse_storage(SparseStorageFactory<IndexType, ElementType>::create(sparse_format, dimension, capacity, use_regularization)) {
    }
-   
+
    template <typename IndexType, typename ElementType>
-   // TODO fix. We need to scan through all the columns
    inline ElementType SymmetricMatrix<IndexType, ElementType>::smallest_diagonal_entry(size_t max_dimension) const {
-      ElementType smallest_entry = INF<ElementType>;
+      // diagonal entries might be at several locations and must be accumulated
+      // TODO preallocate this vector somewhere
+      std::vector<ElementType> diagonal_entries(max_dimension, ElementType(0));
+
       for (const auto [row_index, column_index, element]: *this->sparse_storage) {
          if (row_index == column_index && row_index < max_dimension) {
-            smallest_entry = std::min(smallest_entry, element);
+            diagonal_entries[row_index] += element;
          }
       }
-      if (smallest_entry == INF<ElementType>) {
-         smallest_entry = ElementType(0);
-      }
-      return smallest_entry;
+      return *std::min_element(diagonal_entries.begin(), diagonal_entries.end());
    }
 
    template <typename IndexType, typename ElementType>