Implement one-phase IPM

docs/algorithms.md

@@ -181,8 +181,10 @@ Section 6 of [^3] describes how to check for local infeasibility.
 
 ## Works cited
 
-[^1]: Nocedal, J. and Wright, S. "Numerical Optimization", 2nd. ed., Ch. 19. Springer, 2006.
+[^1]: Hinder, O. and Ye, Y. "A one-phase interior point method for nonconvex optimization", 2018. [https://arxiv.org/pdf/1801.03072.pdf](https://arxiv.org/pdf/1801.03072.pdf)
 
-[^2]: Wächter, A. and Biegler, L. "On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming", 2005. [http://cepac.cheme.cmu.edu/pasilectures/biegler/ipopt.pdf](http://cepac.cheme.cmu.edu/pasilectures/biegler/ipopt.pdf)
+[^2]: Nocedal, J. and Wright, S. "Numerical Optimization", 2nd. ed., Ch. 19. Springer, 2006.
 
-[^3]: Byrd, R. and Nocedal J. and Waltz R. "KNITRO: An Integrated Package for Nonlinear Optimization", 2005. [https://users.iems.northwestern.edu/~nocedal/PDFfiles/integrated.pdf](https://users.iems.northwestern.edu/~nocedal/PDFfiles/integrated.pdf)
+[^3]: Wächter, A. and Biegler, L. "On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming", 2005. [http://cepac.cheme.cmu.edu/pasilectures/biegler/ipopt.pdf](http://cepac.cheme.cmu.edu/pasilectures/biegler/ipopt.pdf)
+
+[^4]: Byrd, R. and Nocedal J. and Waltz R. "KNITRO: An Integrated Package for Nonlinear Optimization", 2005. [https://users.iems.northwestern.edu/~nocedal/PDFfiles/integrated.pdf](https://users.iems.northwestern.edu/~nocedal/PDFfiles/integrated.pdf)

include/sleipnir/optimization/OptimizationProblem.hpp

@@ -189,10 +189,12 @@ class SLEIPNIR_DLLEXPORT OptimizationProblem {
           std::max(status.equalityConstraintType, c.Type());
     }
 
-    m_equalityConstraints.reserve(m_equalityConstraints.size() +
-                                  constraint.constraints.size());
-    std::copy(constraint.constraints.begin(), constraint.constraints.end(),
-              std::back_inserter(m_equalityConstraints));
+    m_inequalityConstraints.reserve(m_inequalityConstraints.size() +
+                                    2 * constraint.constraints.size());
+    for (const auto& c : constraint.constraints) {
+      m_inequalityConstraints.emplace_back(c);
+      m_inequalityConstraints.emplace_back(-c);
+    }
   }
 
   /**
@@ -208,10 +210,12 @@ class SLEIPNIR_DLLEXPORT OptimizationProblem {
           std::max(status.equalityConstraintType, c.Type());
     }
 
-    m_equalityConstraints.reserve(m_equalityConstraints.size() +
-                                  constraint.constraints.size());
-    std::copy(constraint.constraints.begin(), constraint.constraints.end(),
-              std::back_inserter(m_equalityConstraints));
+    m_inequalityConstraints.reserve(m_inequalityConstraints.size() +
+                                    2 * constraint.constraints.size());
+    for (const auto& c : constraint.constraints) {
+      m_inequalityConstraints.emplace_back(c);
+      m_inequalityConstraints.emplace_back(-c);
+    }
   }
 
   /**
@@ -304,11 +308,21 @@ class SLEIPNIR_DLLEXPORT OptimizationProblem {
       return status;
     }
 
+    // Turn each equality constraint into two inequality constraints:
+    //
+    //   cₑ(x) = 0
+    //
+    //   cₑ ≥ 0
+    //   cₑ ≤ 0
+    //
+    //   cₑ ≥ 0
+    //   −cₑ ≥ 0
+
     // Solve the optimization problem
     Eigen::VectorXd s = Eigen::VectorXd::Ones(m_inequalityConstraints.size());
     InteriorPoint(m_decisionVariables, m_equalityConstraints,
-                  m_inequalityConstraints, m_f.value(), m_callback, config,
-                  false, x, s, &status);
+                  m_inequalityConstraints, m_f.value(), m_callback, config, x,
+                  s, &status);
 
     if (config.diagnostics) {
       sleipnir::println("Exit condition: {}", ToMessage(status.exitCondition));

include/sleipnir/optimization/SolverIterationInfo.hpp

@@ -26,9 +26,6 @@ struct SolverIterationInfo {
   /// The Hessian of the Lagrangian.
   const Eigen::SparseMatrix<double>& H;
 
-  /// The equality constraint Jacobian.
-  const Eigen::SparseMatrix<double>& A_e;
-
   /// The inequality constraint Jacobian.
   const Eigen::SparseMatrix<double>& A_i;
 };

include/sleipnir/optimization/solver/InteriorPoint.hpp

@@ -36,8 +36,6 @@ are the inequality constraints.
 @param[in] f The cost function.
 @param[in] callback The user callback.
 @param[in] config Configuration options for the solver.
-@param[in] feasibilityRestoration Whether to use feasibility restoration instead
-  of the normal algorithm.
 @param[in,out] x The initial guess and output location for the decision
   variables.
 @param[in,out] s The initial guess and output location for the inequality
@@ -49,7 +47,7 @@ SLEIPNIR_DLLEXPORT void InteriorPoint(
     std::span<Variable> equalityConstraints,
     std::span<Variable> inequalityConstraints, Variable& f,
     function_ref<bool(const SolverIterationInfo& info)> callback,
-    const SolverConfig& config, bool feasibilityRestoration, Eigen::VectorXd& x,
-    Eigen::VectorXd& s, SolverStatus* status);
+    const SolverConfig& config, Eigen::VectorXd& x, Eigen::VectorXd& s,
+    SolverStatus* status);
 
 }  // namespace sleipnir

jormungandr/cpp/Docstrings.hpp

@@ -265,10 +265,6 @@ Parameter ``callback``:
 Parameter ``config``:
     Configuration options for the solver.
 
-Parameter ``feasibility_restoration``:
-    Whether to use feasibility restoration instead of the normal
-    algorithm.
-
 Parameter ``x``:
     The initial guess and output location for the decision variables.
 
@@ -907,8 +903,6 @@ static const char *__doc_sleipnir_SolverExitCondition_kTooFewDOFs = R"doc(The so
 
 static const char *__doc_sleipnir_SolverIterationInfo = R"doc(Solver iteration information exposed to a user callback.)doc";
 
-static const char *__doc_sleipnir_SolverIterationInfo_A_e = R"doc(The equality constraint Jacobian.)doc";
-
 static const char *__doc_sleipnir_SolverIterationInfo_A_i = R"doc(The inequality constraint Jacobian.)doc";
 
 static const char *__doc_sleipnir_SolverIterationInfo_H = R"doc(The Hessian of the Lagrangian.)doc";

jormungandr/cpp/optimization/BindSolverIterationInfo.cpp

@@ -26,9 +26,6 @@ void BindSolverIterationInfo(nb::class_<SolverIterationInfo>& cls) {
   cls.def_prop_ro(
       "H", [](const SolverIterationInfo& self) { return self.H; },
       DOC(sleipnir, SolverIterationInfo, H));
-  cls.def_prop_ro(
-      "A_e", [](const SolverIterationInfo& self) { return self.A_e; },
-      DOC(sleipnir, SolverIterationInfo, A_e));
   cls.def_prop_ro(
       "A_i", [](const SolverIterationInfo& self) { return self.A_i; },
       DOC(sleipnir, SolverIterationInfo, A_i));

src/optimization/RegularizedLDLT.hpp

@@ -2,7 +2,6 @@
 
 #pragma once
 
-#include <cmath>
 #include <cstddef>
 
 #include <Eigen/Core>
@@ -37,22 +36,15 @@ class RegularizedLDLT {
    * Computes the regularized LDLT factorization of a matrix.
    *
    * @param lhs Left-hand side of the system.
-   * @param numEqualityConstraints The number of equality constraints in the
-   *   system.
-   * @param μ The barrier parameter for the current interior-point iteration.
    */
-  void Compute(const Eigen::SparseMatrix<double>& lhs,
-               size_t numEqualityConstraints, double μ) {
+  void Compute(const Eigen::SparseMatrix<double>& lhs) {
     // The regularization procedure is based on algorithm B.1 of [1]
-    m_numDecisionVariables = lhs.rows() - numEqualityConstraints;
-    m_numEqualityConstraints = numEqualityConstraints;
+    m_numDecisionVariables = lhs.rows();
 
-    const Inertia idealInertia{m_numDecisionVariables, m_numEqualityConstraints,
-                               0};
+    const Inertia idealInertia{m_numDecisionVariables, 0, 0};
     Inertia inertia;
 
     double δ = 0.0;
-    double γ = 0.0;
 
     AnalyzePattern(lhs);
     m_solver.factorize(lhs);
@@ -67,13 +59,6 @@ class RegularizedLDLT {
       }
     }
 
-    // If the decomposition succeeded and the inertia has some zero eigenvalues,
-    // or the decomposition failed, regularize the equality constraints
-    if ((m_solver.info() == Eigen::Success && inertia.zero > 0) ||
-        m_solver.info() != Eigen::Success) {
-      γ = 1e-8 * std::pow(μ, 0.25);
-    }
-
     // Also regularize the Hessian. If the Hessian wasn't regularized in a
     // previous run of Compute(), start at a small value of δ. Otherwise,
     // attempt a δ half as big as the previous run so δ can trend downwards over
@@ -87,9 +72,8 @@ class RegularizedLDLT {
     while (true) {
       // Regularize lhs by adding a multiple of the identity matrix
       //
-      // lhs = [H + AᵢᵀΣAᵢ + δI   Aₑᵀ]
-      //       [       Aₑ        −γI ]
-      Eigen::SparseMatrix<double> lhsReg = lhs + Regularization(δ, γ);
+      // lhs = [H + AᵢᵀΣAᵢ + δI]
+      Eigen::SparseMatrix<double> lhsReg = lhs + Regularization(δ);
       AnalyzePattern(lhsReg);
       m_solver.factorize(lhsReg);
       inertia = Inertia{m_solver};
@@ -132,9 +116,6 @@ class RegularizedLDLT {
   /// The number of decision variables in the system.
   size_t m_numDecisionVariables = 0;
 
-  /// The number of equality constraints in the system.
-  size_t m_numEqualityConstraints = 0;
-
   /// The value of δ from the previous run of Compute().
   double m_δOld = 0.0;
 
@@ -158,19 +139,14 @@ class RegularizedLDLT {
    * Returns regularization matrix.
    *
    * @param δ The Hessian regularization factor.
-   * @param γ The equality constraint Jacobian regularization factor.
    */
-  Eigen::SparseMatrix<double> Regularization(double δ, double γ) {
-    Eigen::VectorXd vec{m_numDecisionVariables + m_numEqualityConstraints};
+  Eigen::SparseMatrix<double> Regularization(double δ) {
+    Eigen::VectorXd vec{m_numDecisionVariables};
     size_t row = 0;
     while (row < m_numDecisionVariables) {
       vec(row) = δ;
       ++row;
     }
-    while (row < m_numDecisionVariables + m_numEqualityConstraints) {
-      vec(row) = -γ;
-      ++row;
-    }
 
     return Eigen::SparseMatrix<double>{vec.asDiagonal()};
   }