diff --git a/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp b/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp index 64254dbb..2236603c 100644 --- a/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp +++ b/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); } diff --git a/uno/ingredients/hessian_models/ConvexifiedHessian.cpp b/uno/ingredients/hessian_models/ConvexifiedHessian.cpp index 3e981ed9..cd9ce44a 100644 --- a/uno/ingredients/hessian_models/ConvexifiedHessian.cpp +++ b/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& 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 \ No newline at end of file +} // namespace diff --git a/uno/linear_algebra/SymmetricMatrix.hpp b/uno/linear_algebra/SymmetricMatrix.hpp index 5094716e..e1394690 100644 --- a/uno/linear_algebra/SymmetricMatrix.hpp +++ b/uno/linear_algebra/SymmetricMatrix.hpp @@ -4,6 +4,7 @@ #ifndef UNO_SYMMETRICMATRIX_H #define UNO_SYMMETRICMATRIX_H +#include #include #include #include @@ -63,20 +64,19 @@ namespace uno { SymmetricMatrix::SymmetricMatrix(size_t dimension, size_t capacity, bool use_regularization, const std::string& sparse_format) : sparse_storage(SparseStorageFactory::create(sparse_format, dimension, capacity, use_regularization)) { } - + template - // TODO fix. We need to scan through all the columns inline ElementType SymmetricMatrix::smallest_diagonal_entry(size_t max_dimension) const { - ElementType smallest_entry = INF; + // diagonal entries might be at several locations and must be accumulated + // TODO preallocate this vector somewhere + std::vector 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) { - smallest_entry = ElementType(0); - } - return smallest_entry; + return *std::min_element(diagonal_entries.begin(), diagonal_entries.end()); } template