Paper deep dive
ODYN: An All-Shifted Non-Interior-Point Method for Quadratic Programming in Robotics and AI
Jose Rojas, Aristotelis Papatheodorou, Sergi Martinez, Andrea Patrizi, Ioannis Havoutis, Carlos Mastalli
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 90%
Last extracted: 7/21/2026, 2:31:54 AM
Summary
The paper introduces ODYN, a novel all-shifted primal-dual non-interior-point quadratic programming (QP) solver designed for robotics and AI applications. It combines all-shifted nonlinear complementarity problem (NCP) functions with the proximal method of multipliers to handle ill-conditioned and degenerate problems without requiring linear independence of constraints. ODYN demonstrates superior warm-start performance and is benchmarked against state-of-the-art solvers on the Maros-Mészáros test set. The paper details its implementation as a backend for SQP-based predictive control (OdynSQP), a differentiable optimization layer for deep learning (ODYNLayer), and a contact-dynamics simulation optimizer (ODYNSim).
Entities (16)
Relation Signals (11)
ODYN → benchmarkedon → Maros-Mészáros test set
confidence 95% · benchmark ODYN on the Maros-Mészáros test set, demonstrating state-of-the-art convergence performance
ODYN → implements → Non-interior-point methods
confidence 95% · ODYN is built upon these principles... Odyn belongs to the class of path-following methods.
ODYN → uses → Proximal method of multipliers
confidence 95% · ODYN combines all-shifted nonlinear complementarity problem (NCP) functions with proximal method of multipliers
ODYN → isbackendfor → OdynSQP
confidence 92% · deploying ODYN as the backend of an SQP-based predictive control framework (OdynSQP)
ODYN → isbackendfor → ODYNLayer
confidence 92% · as the implicitly differentiable optimization layer for deep learning (ODYNLayer)
ODYN → isbackendfor → ODYNSim
confidence 92% · and the optimizer of a contact-dynamics simulation (ODYNSim)
OdynSQP → builton → Crocoddyl
confidence 90% · OdynSQP is developed on top of Crocoddyl, an open-source software framework for model predictive control in robotics.
ODYN → handles → Ill-conditioned problems
confidence 90% · robustly address ill-conditioned and degenerate problems
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:We introduce ODYN, a novel all-shifted primal-dual non-interior-point quadratic programming (QP) solver designed to efficiently handle challenging dense and sparse QPs. ODYN combines all-shifted nonlinear complementarity problem (NCP) functions with proximal method of multipliers to robustly address ill-conditioned and degenerate problems, without requiring linear independence of the constraints. It exhibits strong warm-start performance and is well suited to both general-purpose optimization, and robotics and AI applications, including model-based control, estimation, and kernel-based learning methods. We provide an open-source implementation and benchmark ODYN on the Maros-Mészáros test set, demonstrating state-of-the-art convergence performance in small-to-high-scale problems. The results highlight ODYN's superior warm-starting capabilities, which are critical in sequential and real-time settings common in robotics and AI. These advantages are further demonstrated by deploying ODYN as the backend of an SQP-based predictive control framework (OdynSQP), as the implicitly differentiable optimization layer for deep learning (ODYNLayer), and the optimizer of a contact-dynamics simulation (ODYNSim).
Tags
Links
- Source: https://arxiv.org/abs/2602.16005v2
- Canonical: https://arxiv.org/abs/2602.16005v2
Trouble viewing inline? Open PDF directly →
Full Text
120,020 characters extracted from source content.
Expand or collapse full text
ODYN: An All-Shifted Non-Interior-Point Method for Quadratic Programming in Robotics and AI Jose Rojas† Aristotelis Papatheodorou Sergi Martinez† Andrea Patrizi∗ Ioannis Havoutis Carlos Mastalli† This research was conducted as part of the Advancing MANipulation skills in Legged Robots (AMAN) project, a collaborative project supported by Tata Consultancy Services. (Corresponding author: Carlos Mastalli) Jose Rojas, Sergi Martinez and Carlos Mastalli are part of the Robot Motor Intelligence (RoMI) Lab, Heriot-Watt University, U.K. Aristotelis Papatheodorou and Ioannis Havoutis are part of the Oxford Robotics Institute (ORI), University of Oxford, U.K. Andrea Patrizi is part of Humanoids and Human Centered Mechatronics, Italian Institute of Technology, Italy. Abstract We introduce Odyn, a novel all-shifted primal–dual non-interior-point quadratic programming (QP) solver designed to efficiently handle challenging dense and sparse QPs. Odyn combines all-shifted nonlinear complementarity problem (NCP) functions with the proximal method of multipliers to robustly address ill-conditioned and degenerate problems, without requiring linear independence of the constraints. It exhibits strong warm-start performance and is well-suited to both general-purpose optimization and robotics and AI applications, including model-based control, estimation, and kernel-based learning methods. We provide an open-source implementation and benchmark Odyn on the Maros–Mészáros test set, demonstrating state-of-the-art convergence performance in small-to-high-scale problems. The results highlight Odyn’s superior warm-starting capabilities, which are critical in sequential and real-time settings common in robotics and AI. These advantages are further demonstrated by deploying Odyn as the backend of an SQP-based predictive control framework (OdynSQP), as the implicitly differentiable optimization layer for deep learning (ODYNLayer), and the optimizer of a contact-dynamics simulation (ODYNSim). I Introduction Quadratic programming is a foundational optimization framework in robotics, AI, and a broad range of scientific and engineering disciplines, including operations research. In robotics, quadratic programming (QP) formulations commonly arise in applications such as contact simulation, model predictive control (MPC), whole-body control, state-estimation, and simultaneous localization and mapping (SLAM). Beyond these direct uses, QP solvers also underpin more advanced optimization paradigms, serving as the building blocks of sequential quadratic programming (SQP) methods for nonlinear programs and mixed-integer quadratic programming (MIQP) approaches for problems with discrete decision variables, as detailed in [48] and [31], respectively. Modern numerical QP solvers can be broadly categorized into two main classes: active-set methods and relaxation-based methods. Active-set methods explicitly maintain and update a working set of active constraints, making them well suited to warm-started and small-to-medium-scale problems in which the active set evolves smoothly. Relaxation-based methods, including interior-point methods (IPMs) and augmented Lagrangian methods (ALMs), instead enforce feasibility through barrier or penalty mechanisms, respectively, and are widely adopted for their robustness and scalability. The standard formulation supported by many off-the-shelf QP solvers is defined as follows: min∈ℝn _x ^n 12⊺+⊺ 12\,x Q\,x+c x (1) subject to =, \,x=b, ≤, \,x , where ∈ℝnx ^n denotes the decision vector, ∈ℝn×nQ ^n× n is a symmetric and positive semidefinite quadratic cost matrix (⪰Q 0), ∈ℝnc ^n is a vector of linear cost terms, ∈ℝm×nA ^m× n is the matrix of equality constraints, ∈ℝmb ^m is the right-hand side vector for the equality constraints, ∈ℝp×nG ^p× n is the matrix of inequality constraints, and ∈ℝph ^p is the right-hand side vector for the inequality constraints. For large-scale problems (e.g., n≫100n 100), QP solvers typically leverage sparse backends, relying on sparse linear solvers [23] and representing Q, A, and G as sparse matrices. In contrast, for small to moderately sized problems—frequently encountered in robotics and AI—dense linear solvers and matrix representations tend to be more efficient and convenient. Robotics and AI applications often demand optimizers that support warm-starting, high efficiency, robustness, and scalability. For example, MPC algorithms continuously update motion plans and control actions by solving large, nonlinear optimization problems, typically warm-started from previous solutions. Similarly, warm-start strategies are critical for real-time localization and system identification [39, 40]. Moreover, contact simulators rely on optimizers that must handle rank-deficient contact Jacobians or positive semidefinite Delassus matrices [62], while still benefiting from warm-start information. Yet state-of-the-art solvers such as Gurobi [26], Mosek [46], PiQP [55], ProxQP [3], and OsQP [58] often struggle to satisfy all of these requirements simultaneously. A promising class of algorithms, known as non-interior-point methods (NIPMs), offers an alternative to classical interior-point methods. They are motivated in part by the desire to avoid the strict interior-feasibility requirements imposed by its log-barrier functions. Although these methods scale well to problems with many inequality constraints, they often make robust warm-starting difficult when solving sequences of closely related optimization problems. Odyn is built upon these principles and demonstrates strong capabilities that advance the state of the art in quadratic programming. Figure 1: Overview of Odyn applications in robotics and AI. Odyn serves as the computational core for constrained nonlinear trajectory optimization (OdynSQP), contact-dynamics simulation (ODYNSim), and differentiable optimization layers (ODYNLayer), providing a common optimization backbone across control, simulation, and learning. I-A Contribution We introduce Odyn, a warm-startable QP solver that we characterize as an all-shifted NIPMs—a class of algorithms that remains relatively underexplored in the numerical optimization community, yet holds significant potential for robotics and AI workloads. Our main contributions are: (i) A new NIPM-based QP formulation. We develop a novel solver that blends elements of primal–dual interior-point methods with augmented-Lagrangian techniques. Odyn incorporates three key novelties: a centering-weighted barrier, an all-shifted nonlinear complementarity problem (NCP) formulation, and proximal Lagrangian penalties for equality and inequality constraints. They are critical for robust performance under rank deficiencies, degeneracies, and poor scaling. (i) A complete open-source implementation. We release an efficient C++ implementation with Python bindings.111The code will be available after acceptance. Odyn currently supports dense and sparse backends, multiple floating-point numbers, code generation, and nullspace handling of equality constraints. Moreover, we incorporate ODYNLayer, a framework for training QP layers within neural-network architectures in Pytorch. Finally, OdynSQP is developed on top of Crocoddyl, an open-source software framework for model predictive control in robotics. (i) Benchmarking and robotics/AI evaluation. We benchmark Odyn against state-of-the-art interior-point and augmented-Lagrangian QP solvers using the Maros–Mészáros test set. We further assess its practical utility across robotics and AI use cases: (i) an efficient SQP framework for real-time MPC named OdynSQP, (i) contact simulation solved via ODYNSim, and (i) training differentiable QP layers through the ODYNLayer. These characteristics make Odyn particularly well suited for robotics and AI applications (Fig. 1). More broadly, its strong convergence properties also make it relevant to a wide range of engineering and scientific problems. Below, we review the main classes of available QP solvers, the underlying theory, and their salient characteristics. I Related Work Quadratic programming (QP) plays a central role in robotics, control, and numerical optimization, motivating the development of numerous solvers built on distinct algorithmic paradigms. Modern QP methods can be broadly grouped into three families: active-set methods, augmented Lagrangian methods (penalty-based), and interior-point methods (barrier-based). Each family offers different trade-offs in scalability, warm-start capability, and robustness to ill-conditioning. ALMs, for example, are traditionally designed for equality-constrained problems, although solvers such as Algencan [2] and Lancelot [12] have introduced mechanisms to handle inequality constraints. Their appeal lies in their ability to cope with rank deficiencies in equality constraints while mitigating ill-conditioning. In contrast, IPMs excel at handling large numbers of inequality constraints through barrier functions and typically incorporate equality constraints via primal–dual Karush–Kuhn–Tucker (KKT) systems, as in Ipopt [63]. By avoiding the combinatorial complexity of explicitly tracking active constraints, they offer better scalability than active-set or penalty-based methods. However, their reliance on strict interior feasibility makes warm-starting inherently difficult. Section I-B reviews recent efforts to mitigate this limitation through mechanisms such as exact-penalty formulations and all-shifted complementarity strategies. Motivated by the warm-starting limitations of IPMs, non-interior-point methods (NIPMs) have gained renewed attention. These approaches relax the strict feasibility requirements of traditional IPMs by reformulating complementarity constraints using NCP functions. Section I-C surveys the development of NIPMs, including nonsmooth and path-following variants, and highlights recent hybrid strategies that blend IPMs with ALM-style updates. We begin by reviewing the most prominent QP solvers in the literature. I-A Available QP solvers and their features qpOASES is a classical example of an active-set solver, introduced in [17]. It is tailored for online and real-time applications, such as MPC, and is particularly effective for small- to medium-sized dense QPs. Designed with warm-starting in mind, qpOASES performs well when the problem structure evolves gradually over time. Another widely used active-set solver in robotics is eiquadprog [57], which implements the Goldfarb–Idnani dual active-set algorithm [22]. Originally developed within the Stack-of-Tasks framework, eiquadprog is commonly employed in whole-body control and legged locomotion. Like qpOASES, it is most effective for dense problems and benefits from warm-starting. However, active-set solvers such as qpOASES and eiquadprog tend to struggle with large-scale or highly degenerate problems. They may require a large number of iterations when the active set changes substantially between successive iterates. This behavior stems from the combinatorial complexity of identifying the optimal set of active constraints. OsQP belongs to the family of operator-splitting methods and leverages the alternating direction method of multipliers (ADMM) algorithm [58]. This places OsQP within the class of ALMs. It is known for its robustness and scalability to large, sparse problems and supports both structure-exploiting updates and warm-starting. OsQP also offers code generation capabilities for maximal runtime efficiency in CPUs [5], and can be executed on GPUs [54]. However, as a first-order method, its key limitation lies in the fact that typically achieves only moderate solution accuracy. It may also struggle with problems that exhibit strong ill-conditioning or involve tight constraints. ProxQP is another ALM-based solver, with notable distinctions inspired from Lancelot as highlighted in [3]. It employs a proximal augmented Lagrangian method and incorporates second-order information, enabling more accurate and robust solutions compared to first-order methods such as OsQP. ProxQP is designed with a modular architecture and offers a differentiable QP layer for integration into deep learning pipelines [4]. It supports warm-starting and performs well for medium-scale problems. However, its has limited robustness to solve severe ill-conditioned problems and does not yet provide strong guarantees for exact infeasibility certification. PiQP was recently introduced in [55]. It adopts a hybrid approach by combining an infeasible primal–dual IPM with the proximal method of multipliers. This combination enhances numerical robustness and enables PiQP to handle ill-conditioned convex QP problems without requiring the linear independence of the constraints—an advantage over traditional IPMs. While this hybridization results in a more robust solver, PiQP inherits the typical limitations of IPMs in lacking efficient warm-starting and may require repeated factorization across iterations. Nonetheless, PiQP demonstrates excellent performance on challenging QP problems, particularly when benchmark it against the Maros–Mészáros test set. The Maros–Mészáros benchmark comprises both separable and non-separable QP problems [38]. Separable problems are characterized by a diagonal Q matrix, containing only squared terms, while non-separable problems include off-diagonal entries, with Q remaining symmetric. In addition to the aforementioned QP solvers, several other notable approaches merit consideration. qpSWIFT is an IPM-based solver tailored for embedded QPs presented in [50]. Commercial solvers such as Gurobi and Mosek provide highly optimized interior-point algorithms and support a broad class of problems, though they are closed-source and not specifically designed for real-time robotics or AI. Hpipm offers a dense, Riccati-based formulation particularly suited for MPC applications [19]. However, as an IPM-type method, it inherits the typical limitations related to warm-starting. I-B Warm-starting IPMs Interior-point (IP) methods face inherent challenges with warm-starting, as their iterates must remain strictly within the inequality-feasible set (or the positive orthant when slack variables are introduced). A key difficulty is that reusing the optimal solution of a previous problem (e.g., in SQP solvers) as an initial point for a related problem often leads to blocked search directions and severe ill-conditioning, making it hard for the solver to identify changes in the active set [24]. A common workaround is to retain one or several unconverged iterates and modify them to satisfy the feasibility requirements of the new problem. Another line of research includes the exact-penalty reformulation of [7], although its reliance on adaptive penalty tuning can introduce numerical instabilities. However, a related approach introduced in [15] avoids penalty parameters but still requires explicit identification of the active constraints. Alternatively, the approach of [56] uses a convex combination of a well-centered reference point and the previous solution to warm-start linear programs (LPs) and second-order cone programs (SOCPs) within the homogeneous self-dual framework. Despite these efforts, progress on robust and general warm-starting mechanisms for IPMs has remained limited. More recently, the method of [21] leverages all-shifted complementarity constraints to relax non-negativity requirements through a projected line search mechanism. This strategy offers a promising alternative to classical warm-starting techniques by allowing iterates to move outside the strict interior region while still preserving the structure required by primal–dual interior-point algorithms. I-C Non-interior-point methods Non-interior-point (NIP) methods were developed as an alternative to IPMs, with the goal of solving constrained optimization problems, including linear and nonlinear complementarity problems, by leveraging NCP functions. NCP functions avoid the need for iterates to remain strictly interior to the feasible set. NIP methods can be broadly classified into two main categories: semismooth and path-following methods. Semismooth NIP algorithms utilize semismooth NCP functions, such as the minimum map and the Fischer–Burmeister function [52]. These approaches have been successfully applied in diverse domains including physics simulation [61, 36], quadratic programming [33], nonlinear complementarity problems [59], and second-order cone complementarity problems [35]. In contrast, path-following NIP algorithms solve a sequence of nonlinear systems defined by smoothed NCP functions, parameterized by a smoothing parameter. The trajectory of optimal solutions for these systems is known as the central path. These algorithms have been explored for solving the linear complementarity problem (LCP) [29, 10], optimal control problems with complementarity constraints [34], and linear programming using predictor-corrector strategies [9, 16]. More recently, the relaxed interior-point method has been introduced, which combines interior-point method barrier terms with augmented Lagrangian penalties. This leads to a path-following NIP algorithm that incorporates a smoothed minimum NCP function, and has been applied to quadratic programming [66] and semidefinite programming [67]. For context, we provide below a summary of representative state-of-the-art QP solvers discussed above. Table I categorizes these solvers according to their underlying algorithmic principles and highlights their main practical features, including backend support, warm-start capabilities, robustness to degeneracy, achievable solution accuracy, and suitability for differentiation. Odyn belongs to the class of path-following methods. Before detailing the Odyn’s NIP algorithm, we first introduce the mathematical foundations underlying the construction of Odyn, some of which are novel. TABLE I: Comparison of representative QP solvers grouped by algorithmic method. Method Solver Algorithm Backend Warm-start Degeneracy Accuracy Differentiable Non-interior Point Odyn All-shifted NIPM Dense / Sparse ✓ ✓ High ✓ Fbstab [33] Semismooth NIPM Sparse ✓ ✓ High ✗ Interior-Point PiQP [55] Proximal IPM Sparse ✗ ✓ High ✗ Mosek [46] Homogeneous self-dual IPM Dense / Sparse ✗ ✓ High ✗ Gurobi [26] Primal–dual IPM Dense / Sparse ✗ ✓ High ✗ qpSWIFT [50] Primal–dual IPM Sparse ✗ ✓ High ✗ Augmented Lagrangian ProxQP [3] Proximal ALM Dense / Sparse ✓ Partial High ✓ OsQP [58] ADMM-based ALM Sparse ✓ Limited Medium ✗ Active-set qpOASES [17] Online active-set Dense ✓ Limited High ✗ eiquadprog [57] Goldfarb–Idnani active-set Dense ✓ Limited High ✗ I Algorithmic Foundations We begin by deriving the optimality conditions and establishing the connections that enable us to construct a non-interior approach from interior-point principles. We then discuss mechanisms for introducing regularity in degenerate optimization problems. Both components draw inspiration from the literature on interior-point and augmented-Lagrangian methods. I-A Overview of numerical challenges We start by simplifying inequality constraints complexity of Eq. 1 as follows min, _x,s 12⊺+⊺ 12\,x Q\,x+c x (2) subject to =, \,x=b, +=, \,x+s=h, ≥. 0. where ∈ℝps ^p is a vector of slack variables. To establish the optimality conditions, we define the Lagrangian of Eq. 2 as ℒ(,,,) (x,s,y,z) =12⊺+⊺+⊺(−) = 12x Qx+c x+y (Ax-b) +⊺(+−), +z (Gx+s-h), (3) where, as introduced earlier, ∈ℝnx ^n represent the primal decision variables, and ∈ℝmy ^m, ∈ℝpz ^p corresponds to the Lagrange multipliers of equality and inequality constraints, respectively. The necessary and sufficient conditions222In convex QP problems, sufficiency follows from convexity of the objective and constraints, ensuring that the KKT conditions characterize optimality. for optimality are given by the saddle point equation ∇ℒ(,,,)= (x,s,y,z)=0. They are also known as first-order necessary conditions: ++⊺+⊺ +c+A y+G z =, =0, (4a) − -b =, =0, (4b) +− +s-h =, =0, (4c) ∘=,(,) =0,(s,z) ≥, 0, (4d) where ∘ denotes an element-wise operation, which in this case represents the Hadamard product. A central challenge in numerical optimization arises from the complementarity constraints described in Eq. 4d. This is because they require the solution to remain in the positive orthant. Moreover, QP problems can become degenerate, ill-posed, or large-scale, each of which introduces significant numerical difficulties. Common sources of degeneration include: (i) Rank-deficient constraint matrices: linear dependencies in A or G, (i) Non-uniqueness of dual solutions: multiple valid values for y and z, (i) Singular or semidefinite Hessians: lack of strict convexity in Q, (iv) Ill-conditioning from small constraint violations or near-zero multipliers: e.g., ≈z 0. These issues make QP problems difficult to solve reliably and efficiently with standard solvers. Complicating matters further, QP problems may also be infeasible. Given the impact of these numerical limitations on algorithms in robotics and AI, we design Odyn to tackle hard QP problems from degenerate systems to warm-started sequences of related problems. This makes it especially suitable for real-time applications in robotics and AI, where robustness and efficiency are critical. Below, we delve into the key insights used in the design of the Odyn QP solver. I-B From interior to non-interior point In this section, we highlight the connections between IPMs and ALMs, which motivate the development of our path-following NIP methods. For clarity of exposition, we omit equality constraints; however, they can be easily integrated. IPMs relax the complementarity constraints by employing the log-barrier functions. Inequality constraints are incorporated into the objective as barrier terms, yielding min, _x,s 12⊺+⊺−μ∑j=1plog(j) 12\,x Q\,x+c x-μ _j=1^p (s_j) (5) subject to +=, \,x+s=h, where μ∈ℝ+μ _+ is the barrier parameter. The associated KKT conditions can be written compactly as (;μ)=[++⊺+−∘−μ 1]=,(,)≥,r(U;μ)= bmatrixQx+c+G z\\ Gx+s-h\\ s -μ\,1 bmatrix=0, (s,z) 0, (6) where =(,,)U=(x,s,z) denotes the collection of primal and dual variables. To find the root of the above equation (i.e., KKT point), we apply the Newton’s method to Eq. 6, which leads to the following linear system [⊺][ΔΔΔ]=−[++⊺+−∘−μ 1] bmatrixQ&G &\\ G&&I\\ &S&Z bmatrix bmatrix \\ \\ bmatrix=-\, bmatrixQx+c+G z\\ Gx+s-h\\ s -μ\,1 bmatrix (7) with =diag()S=diag(s) and =diag()Z=diag(z) denoting the diagonal matrices of slack and dual variables, respectively, and 1 as a vector with ones. Eq. 6 satisfies a relaxed complementarity slackness condition, i.e., ∘=μs z= 1. When solving the parametrized QP problem, each value of μ yields a primal–dual point along the central path. In practice, we recover the solution of Eq. 2 by solving a sequence of unconstrained problems with progressively decreasing barrier parameters, letting μ→0μ→ 0. Smoothing the complementarity conditions offers a more “natural” alternative to the projection mechanisms commonly used in ALMs, as it introduces a continuation procedure that gradually drives μ→0μ→ 0 (see Fig. 2). Figure 2: Log-barrier function as μ approaches 0. The approximation becomes closer to the indicator function. However, logarithmic barriers force iterates to remain strictly within the interior of the feasible region. In addition to restricting the search direction, this constraint complicates warm-starting. These limitations have motivated renewed interest in augmented-Lagrangian approaches for robotics and AI [3]. Nevertheless, we argue that their reliance on projection mechanisms can be restrictive. To address the limitations inherent in both IPMs and ALMs, we begin by introducing consensus variables ∈ℝp ξ ^p: min,, _x,s, ξ 12⊺+⊺ 12\,x Q\,x+c x (8) subject to +=, \,x+s=h, =, = ξ, ≥. ξ 0. We then construct an equality-constrained problem by applying log-barrier functions to the inequality constraints (as in IPMs), i.e., min,, _x,s, ξ 12⊺+⊺−μ∑j=1plog(j) 12\,x Q\,x+c x-μ _j=1^p ( ξ_j) (9) subject to +=, \,x+s=h, =. = ξ. We finally relax the strict interior-point condition by introducing a quadratic penalty to the consensus constraint. This yields to a penalty–barrier function that blends ideas from IPMs and ALMs: min,, _x,s, ξ 12⊺+⊺−μ∑j=1plog(j)+12‖−‖22 12\,x Q\,x+c x-μ _j=1^p ( ξ_j)+ 12\,\|s- ξ\|_2^2 (10) subject to +=, \,x+s=h, =. = ξ. The Lagrangian associated with this penalty–barrier formulation is given by (~;μ) ( U;μ) =12⊺+⊺−μ∑j=1plog(j)+12‖−‖22 = 12x Qx+c x-μ _j=1^p ( ξ_j)+ 12 - ξ _2^2 −⊺(−)+⊺(+−), -w (s- ξ)+z (Gx+s-h), (11) where ∈ℝpw ^p is the Lagrange multipliers associated with the consensus constraint and ~=(,,,,) U=(x,s, ξ,w,z) collects primal and dual decision variables. Therefore, the first-order neccesary conditions are given by: ~(~;μ)=[∇]=[++⊺+−μ∘−1−(−)+−+]=. r( U;μ)= bmatrix _xP\\[3.00003pt] _zP\\[3.00003pt] _wP\\[3.00003pt] _ ξP\\[3.00003pt] _sP bmatrix= bmatrixQx+c+G z\\[3.00003pt] Gx+s-h\\[3.00003pt] s- ξ\\[3.00003pt] -μ\, ξ -1-(s- ξ)+w\\[3.00003pt] s- ξ-w+z bmatrix=0. (12) From these conditions, we observe that we can find a closed-form solution for the consensus variable ξ, as ∇ℒ= _ ξL=0 is a second-order polynomial, i.e., ∘2−(−)∘−μ 1=. ξ 2\;-\;(s-w) ξ\;-\;μ\,1=0. (13) To satisfy the condition > ξ>0, we must pick the positive solution: =(−)+(−)∘2+4μ 12. ξ= (s-w)\;+\; (s-w) 2+4μ\,12. (14) We also observe that =w=z from equations ∇ℒ= _wL=0 and ∇ℒ= _sL=0. This leads to: =(−)+(−)∘2+4μ 12. ξ= (s-z)\;+\; (s-z) 2+4μ\,12. (15) Finally, by injecting the above result into ∇ℒ= _wL=0 we obtain the equivalent KKT conditions: (;μ)=[++⊺+−+−(−)∘2+4μ 1]=.r(U;μ)= bmatrixQx+c+G z\\[3.99994pt] Gx+s-h\\[3.99994pt] s+z- (s-z) 2+4μ\,1 bmatrix=0. (16) Eq. 16 can be rewritten as follows: (;μ)=[++⊺+−ϕ(,;μ)]=,r(U;μ)= bmatrixQx+c+G z\\[3.00003pt] Gx+s-h\\[3.00003pt] φ(s,z;μ) bmatrix=0, (17) where ϕ:ℝp×ℝp→ℝp φ:R^p×R^p ^p is a general NCP function, representing ϕMIN φ^MIN or ϕFB φ^FB. This is because we could recognize the smoothed minimum NCP function ϕMIN:ℝp×ℝp→ℝp φ^MIN:R^p×R^p ^p, i.e., ϕMIN(,;μ)=+−(−)∘2+4μ 1. φ^MIN(s,z;μ)=s+z- (s-z) 2+4μ\,1. (18) Moreover, smoothed NCP functions satisfy the following fundamental property [29]: ϕ(,;μ)=⟺≥,≥,∘=μ 1. φ(s,z;μ)=0\; \;s 0, 0, =μ\,1. (19) where, by injecting the centrality constraint ∘=μs z= 1, we transform a smoothed minimum NCP function as a regularized Fischer–Burmeister NCP, i.e., ϕMIN(,;μ) φ^MIN(s,z;μ) =+−∘2−2∘+∘2+4μ 1 =s+z- \,s 2-2\,s +z 2+4μ\,1\, =+−∘2+∘2+2μ 1 =s+z- \,s 2+z 2+2μ\,1\, =ϕFB(,;μ). = φ^FB(s,z;μ). (20) Finally, the nonlinear system of equations in Eq. 17 is solved by applying Newton’s method: [⊺][ΔΔΔ]=−[++⊺+−ϕ(,;μ)], bmatrixQ&G &\\ G&&I\\ & _z& _s bmatrix bmatrix \\ \\ bmatrix=- bmatrixQx+c+G z\\ Gx+s-h\\ φ(s,z;μ) bmatrix, (21) where ∈ℝp×p _z ^p× p and ∈ℝp×p _s ^p× p are diagonal matrices denoting the NCP Jacobians. To promote convergence, the iterates of NIP methods are typically constrained to a neighborhood of the central path. This is defined (see [27]) as follows: (β)=∈ℝn+2p:‖(;μ)‖≤βμ,μ>0,N\! (β )= \U ^\,n+2p:\; (U;μ) \;≤\;β\,μ,\;μ>0 \, (22) where β determines the radius of the neighborhood, ensuring that the iterates remain sufficiently close to the central path as μ→0μ→ 0. Compared to IPMs, this approach does not impose restrictions on s and z, since the domain of ϕ(,;μ) φ(s,z;μ) includes the positive orthant of these variables. Next, we delve into another insight on which Odyn is based. I-C From degeneracy to regularity We begin by understanding what degeneracy is. To do so, we focus on equality-constrained QP problems: min∈ℝn _x ^n 12⊺+⊺ 12\,x Q\,x+c x (23) subject to =. \,x=b. where its KKT system is defined as: [⊺][ΔΔ]=−[++⊺−]. bmatrixQ&A \\ A& bmatrix bmatrix \\ bmatrix=- bmatrixQx+c+A y\\ Ax-b bmatrix. (24) To be able to factorize this KKT system, we require (i) strong convexity (≻Q 0) and (i) linear independence in the constraint matrix A (A is full rank). Handling this becomes more challenging if the KKT matrix is ill-conditioned. For instance, apart of the extra computational demand, pseudo-inverses are not an option because: (i) Numerical instability occurs when Q or A are ill-conditioned or rank deficient. (i) Their one-shot projection (i.e., poor globalization) does not provides a convergence mechanism or guarantee. Proximal point algorithms introduce regularity by relying on quadratic penalties [53]. They can be used to handle QP problems that are not strict convex or with ill-conditioning by casting the problem as follows maxmin _x _y ℒ(,)+ρd2‖−E‖22−ρe2‖−E‖22, \;\;L(x,y)+ _d2 -x_E _2^2- _e2 -y_E _2^2, (25) where ℒ(,)=12⊺+⊺+⊺(−)L(x,y)= 12\,x Qx+c x+y (Ax-b), and E∈ℝnx_E ^n and E∈ℝmy_E ^m are the primal and dual estimates, respectively, and ρd,ρe∈ℝ+ _d, _e _+ are the proximal penalties. Eq. 25 is equivalent to shifting the constraints by (Ex_E, Ey_E). This is because Eq. 25 can be rewritten as: min, _x,y\;\; 12⊺+⊺+ρd2‖−E‖22+(−)⊺E 12x Qx+c x+ _d2 -x_E _2^2+(Ax-b) y_E +12ρe‖−‖22−12ρe‖−ρe(−E)‖22. + 12 _e -b _2^2- 12 _e -b- _e(y-y_E) _2^2. (26) From the problem in Section I-C, we have the following perturbed KKT conditions: (^;^,^E)=[++⊺+ρd(−E)−ρe(−E)],r( U;\, ρ,\, U_E)= bmatrixQx+c+A y+ _d(x-x_E)\\[3.99994pt] Ax-b- _e(y-y_E) bmatrix, (27) where ^=(,) U=(x,y), ^E=(E,E) U_E=(x_E,y_E), ^=(ρd,ρe) ρ=( _d, _e) denote the stacked primal–dual variables, their estimates, and the regularization parameters, respectively. The last couple of terms of Section I-C can be seen as the Forsgren–Gill primal–dual penalty function associated with the shifted constraints: −ρeE=.Ax-b- _e\,y_E=0. (28) This is also know as proximal point Lagrangian function. Additionally, the term ρd2‖−E‖22 _d2 -x_E _2^2 in Section I-C can be recognized as the proximal point penalty function. Therefore, combining both proximal point functions together with the term E⊺(−)y_E (Ax-b) lead to a shifted primal–dual penalty function, or a primal–dual proximal point penalty and Lagrangian function. The parameters ρd _d, ρe _e, Ex_E and Ey_E provide a central path towards the solution. This defines a convergence mechanism in contrast to pseudo-inverses. Finally, we now need to factorize the following perturbed KKT system: [~⊺−ρe][ΔΔ]=−[++⊺+ρd(−E)−ρe(−E)] bmatrix Q&A \\[3.00003pt] A&- _eI bmatrix bmatrix \\[3.00003pt] bmatrix=- bmatrixQx+c+A y+ _d(x-x_E)\\[3.99994pt] Ax-b- _e(y-y_E) bmatrix (29) with ~=+ρd Q=Q+ _dI denoting the regularized Hessian matrix, which ensures that the sufficient optimality conditions are satisfied even for linear programs and quasi-convex problems. We are now ready to describe how degeneracy in complementarity conditions is handled via an all-shifted mechanism, and to introduce the Odyn QP algorithm. IV All-shifted non-interior point QP (Odyn) This section presents the mathematical formulation and algorithmic design of Odyn. We start by describing our novel all-shifted NCP functions. IV-A All-shifted NCP functions To provide regularity to the nonlinear complementarity conditions and effectively handle degeneracy, we formulate Eq. 10 as: max,,min, _x,s, ξ _w,z (~;μ)+ρn4‖−E‖22 ( U;μ)+ _n4 -s_E _2^2 (30) −ρn4‖−E‖22−ρn2(−E)⊺(−E), - _n4 -w_E _2^2- _n2(s-s_E) (w-w_E), where E∈ℝps_E ^p and E∈ℝpw_E ^p are the estimates of the slack variables and consensus multipliers, respectively. Note that we recover the optimality conditions of Eq. 10 when E=∗s_E=s^* and E=∗w_E=w^* and, for the sake of simplicity, we restrict our attention to inequality constraints. Eq. 30 yields to the following perturbed NCP Lagrangian: ~(~;μ,ρn,E,E) P( U;μ, _n,s_E,w_E) =(~;μ)+ρn4‖−E‖22 =P( U;μ)+ _n4\|s-s_E\|_2^2 (31) −ρn4‖−E‖22−ρn2(−E)⊺(−E) -50.00008pt- _n4\|w-w_E\|_2^2- _n2(s-s_E) (w-w_E) in which its first-order necessary conditions are obtained by taking the following derivatives: ∇~ _x\, P =++⊺, =Qx+c+G z, (32) ∇~ _y\, P =+−, =Gx+s-h, ∇~ _w\, P =−+ρn2(−E)+ρn2(−E), =s- ξ+ _n2(s-s_E)+ _n2(w-w_E), ∇~ _ ξ\, P =−μ−1−(−)+, =-\,μ\, ξ^-1-(s- ξ)+w, ∇~ _s\, P =−+ρn2(−E)−ρn2(−E)−+. =s- ξ+ _n2(s-s_E)- _n2(w-w_E)-w+z. Similarly to Section I-B, we find a closed-form solution for the consensus variable ξ. Additionally, we observe that Eq. 32 is the perturbed version of Eq. 12. From ∇~= _w P=0 and ∇~= _s P=0, we obtain a perturbed dual consensus constraint −ρn(−E)z-w- _n(w-w_E). Therefore, the first-order necessary conditions in Eq. 32 can be equivalently rewritten as: (;μ)=[++⊺+−ρn(−E)ϕ(,;μ)+ρn(−E)+ρn(−E)],r(U;μ)= bmatrixQx+c+G z\\[3.99994pt] Gx+s-h\\[3.99994pt] z-w- _n(w-w_E)\\[3.99994pt] φ(s,w;μ)+ _n(s-s_E)+ _n(w-w_E) bmatrix, (33) in which, by substituting the solution ξ in ∇ℒ= _wL=0, the all-shifted NCP functions is obtained as follows ϕ(,;μ)+ρn(−E)+ρn(−E). φ(s,w;μ)+ _n(s-s_E)+ _n(w-w_E). (34) Including this regularity helps dealing with degenerated complementarity conditions, i.e., for some j we may have j=j=0s_j=z_j=0. For instance, if we do not include this regularization, the NCP becomes flat nears convergence (i.e, μ=0μ=0). Additionally, similar to Eq. 19 the all-shifted NCP functions satisfy the following property: ϕ(,;μ)=⟺≥,≥,∘=μ 1+¯n∘(E−+E−), φ(s,w;μ)=0\; \; casess 0, 0,\\[4.0pt] s =μ\,1+ ρ_n (s_E-s+w_E-w ), cases (35) where ¯n=ρn2 ρ_n= _n2s. This property is demonstrated in A. All-shifted complementarity constraints. We describe a complete algorithm for solving degenerated QP problems in the next section. IV-B All-shifted non-interior quadratic programming We combine the all-shifted NCP functions with proximal penalties and Lagrangian terms to develop our path-following NIP method. In particular, Odyn solves the following perturbed saddle-point problem: max,,min,, _x,s, ξ\; _w,y,z ~(~;μ,ρn,E,E)+ρd2‖−E‖22 P( U;μ, _n,s_E,w_E)+ _d2 -x_E _2^2 (36) −ρe2‖−E‖22−ρi2‖−E‖22, - _e2 -y_E _2^2- _i2 -z_E _2^2, where ρi∈ℝ+ _i _+ and E∈ℝpz_E ^p are the inequality regularization parameter and estimates, respectively, and including inequality constraints in ~(~;μ,ρn,E,E) P( U;μ, _n,s_E,w_E) with ~=(,,,,,) U=(x,s, ξ,w,y,z). Its KKT conditions are given by: [~d~e~i~c~g]=[++⊺+⊺+ρd(−E)−ρe(−E)+−ρi(−E)−ρn(−E)ϕ(,;μ)+ρn(−E)+ρn(−E)], bmatrix r_d\\[5.0pt] r_e\\[5.0pt] r_i\\[5.0pt] r_c\\[5.0pt] r_g bmatrix= bmatrixQx+c+A y+G z+ _d(x-x_E)\\[5.0pt] Ax-b- _e(y-y_E)\\[5.0pt] Gx+s-h- _i(z-z_E)\\[5.0pt] z-w- _n(w-w_E)\\[5.0pt] φ(s,w;μ)\;+\; _n(s-s_E)+ _n(w-w_E) bmatrix, (37) where ~d∈ℝn r_d ^n is the perturbed dual residual, ~e∈ℝm r_e ^m, ~i∈ℝp r_i ^p and ~c∈ℝp r_c ^p are the perturbed residual vectors associated to equality, inequality and consensus constraints, respectively. In practice, we observe that the proximal term associated with the consensus constraint ~c r_c has a negligible effect. Therefore, we further condense the KKT conditions by substituting =w=z into Eq. 38, yielding: ~(;)=[~d~e~i~g]=[++⊺+⊺+ρd(−E)−ρe(−E)+−ρi(−E)ϕ(,;μ)+ρn(−E)+ρn(−E)] r(U; )= bmatrix r_d\\ r_e\\ r_i\\ r_g bmatrix= bmatrixQx+c+A y+G z+ _d(x-x_E)\\ Ax-b- _e(y-y_E)\\ Gx+s-h- _i(z-z_E)\\ φ(s,z;μ)+ _n(s-s_E)+ _n(z-z_E) bmatrix (38) IV-C Search direction To compute the KKT point in Eq. 38, we apply the Newton’s method. As the Newton’s method iterates over a linear approximation of the NCP functions, this produces a search direction at each iteration obtained by factorizing the following saddle-point system: [~⊺−ρe−ρi~~][ΔΔΔΔ]=−[~d~e~i~g], bmatrix Q&A &G &\\ A&- _eI&&\\ G&&- _iI&I\\ && _z& _s bmatrix bmatrix \\ \\ \\ bmatrix=- bmatrix r_d\\ r_e\\ r_i\\ r_g bmatrix, (39) where ~=+ρn _s= _s+ _nI and ~=+ρn _z= _z+ _nI are the regularized NCP Jacobians with respect to the slack variables and inequality multipliers, respectively. To increase computational efficiency, we could condense Eq. 39 as follows: [~⊺−][ΔΔ]=−[~d~p], bmatrix Q&J \\[3.00003pt] J&-\,D bmatrix bmatrix \\[3.00003pt] λ bmatrix=- bmatrix r_d\\[3.00003pt] r_p bmatrix, (40) where Δ=[Δ⊺Δ⊺]⊺∈ℝm+p λ= bmatrix & bmatrix^\! ^m+p is the stacked search direction for equality and inequality multipliers, ~p=[~e⊺¯i⊺]⊺∈ℝm+p r_p= bmatrix r_e^\! & r_i^\! bmatrix^\! ^m+p is the right-hand-side term associated to the constraint residuals (i.e., the perturbed residuals), ¯i=~i− r_i= r_i-q is the condensed right-side term used to compute Δ , =~−1~P= _s^-1\, _z is a diagonal matrix, =~−1~gq= _s^-1 r_g, ⊺=[⊺]⊺J = bmatrixA &G bmatrix stacks the Jacobian matrices of the constraints, and the diagonal matrix D is defined as: =[ρe~],D= bmatrix _e\,I&\\[3.00003pt] & P bmatrix, (41) with ~=+ρi P=P+ _iI. The system in Eq. 40 can be factorized using a condensed Hessian matrix =~+⊺−1≻H= Q+J D^-1J 0. Since H is positive definite, we employ a Cholesky factorization, exploiting the trivial inversion of the diagonal matrix D. This yields the following search directions: Δ =−1(~d+⊺−1~p), =-\,H^-1\! ( r_d+J^\! D^-1 r_p ), (42a) Δ λ =−1(Δ+~p), =D^-1\! (J\, + r_p ), (42b) Δ =−(Δ+). =- (P\, +q ). (42c) We refer to this approach as the condensed-KKT formulation. In contrast, the full-KKT system in Eq. 39 can be factorized using an LDL⊤ decomposition or a quasi-definite LDT scheme [23]. We now describe the step-acceptance criteria. IV-D Step acceptance A candidate solution is obtained from Δ=(Δ,Δ,Δ,Δ) =( , , , ) via a line search procedure, i.e, +=+αΔU^+=U+α U, where α∈(0,1]α∈(0,1] is the step length parameter. To assess a candidate solution, we employ a natural merit function defined as follows: ℳ()=12(~d⊺~d+~p⊺~p+~g⊺~g),M(U)= 12( r_d r_d+ r_p r_p+ r_g r_g), (43) where ~p=[~e⊺~i⊺]⊺ r_p= bmatrix r_e & r_i bmatrix is the concatenation of the perturbed equality and inequality constraint residuals. Although NCP functions temporarily allow iterates to leave the positive orthant, the feasible neighborhood they define is typically narrow, limiting the step lengths accepted during line search. To mitigate this limitation, we develop a flexible backtracking line search based on a relaxed Armijo condition. Specifically, for a candidate step length α, we require: ℳ(+αΔ)≤γℳ()+αηDℳ()⊺Δ,M\! (U+α\, )\;≤\;γ\,M(U)\;+\;α\,η\;DM(U)^\! , (44) where γ>0γ>0 is a scaling parameter that controls the degree of non-monotonicity and η>0η>0 is a small positive constant. This condition draws inspiration from the non-monotone line search in [65], but differs in that we use a scaled version of the current merit function as the reference value rather than a moving average. Additionally, the natural merit function ℳM is evaluated using the perturbed KKT residuals. IV-E Centering algorithm The barrier or smoothing parameter for the NCP functions is computed as proposed in [29], i.e., μ=‖min(,)‖22p,μ= \| (s,z)\|_2^2p, (45) where p is the number of inequality constraints. This formulation ensures that μ reflects the degree of complementarity violation. The smoothing parameter is updated using: μ+=max(μmin,σμ),μ^+= ( _ ,σμ), (46) where σ∈(0,1)σ∈(0,1) is a reduction factor, often referred to as the centering parameter in the IPM literature [48]. /* compute predictor for the complementarity measure */ 1 μ←‖min(,)‖2/pμ←\| (s,\,z)\|^2/p μ←max(μmin,σμ)μ← ( _ ,\,σμ) /* update centering parameter */ 2 if ‖(+;)‖≤θu‖(;)‖\|r(U^+; )\|≤ _u\|r(U; )\| then 3 μ←δμ+μ← _μ^+\,μ 4 σ←max(σmin,σ−δ−σ)σ← ( _ ,\,σ-δ^-σ) 5 end if 6else 7 σ←min(σmax,σ+δ+(1−σ))σ← ( _ ,\,σ+δ^+(1-σ)) 8 end if /* safeguard against growth */ 9 μ←min(μ,μold)μ← (μ,\, _old) Algorithm 1 Centering algorithm The centering parameter σ is adjusted dynamically according to a trust-region–inspired rule, as detailed in Algorithm 1. Our strategy resembles the heuristic used in PiQP [55] for updating proximal regularization parameters. However, our update is less aggressive: we apply fixed decrease and increase factors based on the observed progress. Specifically, if the KKT residual +r^+ is sufficiently reduced relative to the previous residual r, we decrease the centering parameter as σ+=max(σmin,σ−δ−σ),σ^+= \! ( _ ,\ σ-δ^-σ ), (47) or increase the centering parameters to promote larger steps as σ+=min(σmax,σ+δ+(1−σ)),σ^+= ( _max,\ σ+δ^+(1-σ)), (48) where σmin _ and σmax _ are the minimum and maximum allowed value for the centering parameter, respectively, with δ−δ^- and δ+δ^+ as the decrease and increase factors. In Algorithm 1, the barrier parameter μ is enforced to decrease monotonically. If an increase in the complementarity residual is detected, μ is clamped to its previous value μold _ old to prevent growth in the barrier parameter. Note that =(μ,,E) =(μ, ρ,U_E) denotes the collection of Odyn parameters IV-F Proximal primal–dual estimates Proximal methods typically use either the previous or the current primal–dual iterate as the reference point for the proximal terms [3, 55]. In contrast, we introduce an interpolated estimate between these two choices. Specifically, the updated primal–dual references are computed as: E _E =++θd(−+), =x^++ _d(x-x^+), E _E =++θe(−+), =y^++ _e(y-y^+), E _E =++θi(−+), =z^++ _i(z-z^+), E _E =++θi(−+), =s^++ _i(s-s^+), (49) where θd,θe,θi∈[0,1] _d, _e, _i∈[0,1] are interpolation parameters for dual, equality, inequality, and NCP feasibility, respectively. These interpolation parameters are updated dynamically using a trust-region-inspired strategy. When a sufficient reduction in relative feasibility is observed (i.e., measured against a threshold factor θl _l), the corresponding interpolation parameter is decreased exponentially. The general update rule for each θl _l is given by: θl=(θl−θ−θl;θmin), _l=T( _l-θ^- _l; _ ), (50) where θ−>0θ^->0 is the decay factor and (x;t)T(x;t) is a hard threshold function that sets its input x to zero if it falls below a minimum threshold t. Conversely, if the relative feasibility degrades beyond an upper threshold factor θu _u, the interpolation parameter is increased: θl=min(1,θl+θ+(1−θl)), _l= (1, _l+θ^+(1- _l)), (51) where θ+>0θ^+>0 is the growth factor. The update logic is summarized in Algorithm 2. /* dual step length and regularization update */ 1 if ‖d+‖∞≤θl‖d‖∞\|r_d^+\|_∞≤ _l\|r_d\|_∞ then 2 θd←(max(0,θd−θ−θd);θmin) _d ( (0,\, _d-θ^- _d); _ ) 3 ρd←min(ρd⋅δ−1,ρd,min) _d← ( _d·δ^-1,\, _d, ) 4 end if 5else if ‖d+‖∞≥θu‖d‖∞\|r_d^+\|_∞≥ _u\|r_d\|_∞ then 6 θd←min(1,θd+θ+(1−θd)) _d← (1,\, _d+θ^+(1- _d)) 7 end if /* primal step length and regularization update */ 8 if ‖p+‖∞≤θl‖p‖∞\|r_p^+\|_∞≤ _l\|r_p\|_∞ then 9 θp←(max(0,θp−θ−θp);θmin) _p ( (0,\, _p-θ^- _p); _ ) 10 ρp←min(ρp⋅δ−1,ρp,min) _p← ( _p·δ^-1,\, _p, ) 11 end if 12else if ‖p+‖∞≥θu‖p‖∞\|r_p^+\|_∞≥ _u\|r_p\|_∞ then 13 θp←min(1,θp+θ+(1−θp)) _p← (1,\, _p+θ^+(1- _p)) 14 end if /* NCP regularization update */ 15 if ‖g+‖∞≤θl‖g‖∞\|r_g^+\|_∞≤ _l\|r_g\|_∞ then 16 ρn←min(ρn⋅δ−1,ρn,min) _n← ( _n·δ^-1,\, _n, ) 17 end if Algorithm 2 Compute step length estimates IV-G Neighborhood termination criteria The neighborhood condition typically used in NIP algorithms (see Eq. 22) may require several Newton iterations to be satisfied before the smoothing parameter can be reduced. With hard QP problems, however, this condition can become difficult to meet. This is because it does not account for the scale of the problem and effectively imposes an absolute tolerance on a theoretical subproblem. To address this, we solve the perturbed primal–dual KKT system using fixed values of the smoothing parameter μ and the regularization terms introduced in Eq. 40. Unlike primal–dual ALMs such as ProxQP [3], which rely on an absolute tolerance to terminate the parameterized subproblem, we instead employ a relaxed relative neighborhood condition defined as: (τ)=:∥~(+α;) (τ)= \U:\| r(U^+α; ) ∥∞≤θ∥~(;)∥∞+βμ, \|_∞≤θ\| r(U; )\|_∞+βμ \, (52) where =(β,θ) τ=(β,θ), =(μ,,E) =(μ, ρ,U_E) denotes the collection of Odyn parameters, β>0β>0 is a constant, and θ∈(0,1)θ∈(0,1) specifies the required relative decrease in the perturbed residual norm. When the condition in Eq. 52 is satisfied, we update the smoothing and regularization parameters as well as the proximal primal–dual estimates. IV-H Stopping criteria Inspired by PiQP [55], we employ a termination criterion that combines relative and absolute tolerances. The Odyn QP solver terminates when the primal and dual residuals fall below a threshold defined as the sum of an absolute tolerance ϵa>0 _a>0 and a relative tolerance ϵr≥0 _r≥ 0 scaled by the problem data. Specifically, we require: ‖d‖∞ \|r_d\|_∞ ≤ϵa+ϵrmax(‖∞,‖∞,‖⊤‖∞,‖⊤‖∞), ≤ _a+ _r \! (\|Qx\|_∞,\|c\|_∞,\|A y\|_∞,\|G z\|_∞ ), (53) ‖e‖∞ \|r_e\|_∞ ≤ϵa+ϵrmax(‖∞,‖∞), ≤ _a+ _r \! (\|Ax\|_∞,\|b\|_∞ ), ‖i‖∞ \|r_i\|_∞ ≤ϵa+ϵrmax(‖∞,‖∞,‖∞), ≤ _a+ _r \! (\|Gx\|_∞,\|s\|_∞,\|h\|_∞ ), ‖n‖∞ \|r_n\|_∞ ≤ϵa+ϵrmax(‖∞,‖∞,‖g‖∞), ≤ _a+ _r \! (\|s\|_∞,\|z\|_∞,\|r_g\|_∞ ), where n=min(,)r_n= (s,z) denotes the minimum function, and ϵa>0 _a>0 and ϵr≥0 _r≥ 0 are the absolute and relative tolerances used in the stopping criterion. We use the problem residuals dr_d, er_e, and ir_i, together with gr_g, to scale the NCP-based feasibility measure. IV-I Infeasibility detection Primal and dual infeasibility are certified using Farkas-type conditions, following the approaches in [6, 32]. In Odyn, a QP is declared dual infeasible when the following termination criterion is satisfied: ‖Δ‖∞ \|Q \|_∞ ≤ϵd‖Δ‖∞, ≤ _d\| \|_∞, ⊺ x <−ϵd, <- _d, ‖Δ‖∞ \|A \|_∞ ≤ϵd‖Δ‖∞, ≤ _d\| \|_∞, ‖Δ+Δ‖∞ \|G + \|_∞ ≤ϵd‖Δ‖∞, ≤ _d\| \|_∞, (54) where ϵp>0 _p>0 is a dual infeasibility certification tolerance. Moreover, the primal infeasible is certified when the following condition is met: ‖⊺Δ+⊺Δ‖∞≤ϵp(‖Δ‖∞+‖Δ‖∞),\|A +G \|_∞≤ _p(\| \|_∞+\| \|_∞), ⊺Δ+⊺Δ<−ϵp,b +h <- _p, (55) where ϵp>0 _p>0 is a primal infeasibility certification tolerance. IV-J Overall algorithm and software The overall algorithm of Odyn QP solver is presented in Algorithm 3. /* Odyn QP iterations */ 1 for i←0i← 0 to max_itersmax\_iters do 2 search direction: Δ Eq. 42 3 flexible Armijo line search: α Eq. 44 4 if ‖~(+α;)‖∞≤θ‖~(;)‖∞+βμ\| r(U^+α; )\|_∞≤θ\| r(U; )\|_∞+βμ then 5 update centering Algorithm 1 6 update step-length estimates Algorithm 2 7 update proximal primal–dual estimates Section IV-F 8 end if 9 test stopping criteria Eq. 53 10 end for Algorithm 3 Odyn QP solver Odyn is distributed as an open-source implementation.333The code will be available after acceptance. The solver is written in C++ with Python bindings and relies on the Eigen library [25] for efficient linear-algebra operations in both its dense and sparse backends. Formulating and solving a QP with Odyn involves four main components: Model, Data, Params, and Solver. All components, except for the Params structure, are specialized for both dense and sparse backends. Below, we illustrate how to set up and solve QP problems for each backend in Python. Additionally, the default hyperparameters employed by Odyn (using double precision numbers) are outlined in Table I. For single precision numbers, we increase ρ and μmin _ by 10410^4. Next, we introduce applications of Odyn in the contexts of predictive control, deep learning, and simulation. ⬇ 1import odyn 2 3# Formulate and solve a dense QP problem. 4# We start creating a random QP model, data, as well as a QP solver and parameters. 5# The ’print_level’ parameter controls the verbosity of the output, with ’odyn.VerboseLevel.High’ being the most verbose. 6model = odyn.DenseModel.Random(10, 5, 3) 7data = model.createData() 8params = odyn.Params() 9solver = odyn.DenseQP() 10solver.solve(model, data, params, print_level=odyn.VerboseLevel.High) 11 12# If we want to run the sparse backend for this random QP problem, we can just convert the model and solver to sparse format. 13sparse_model = model.toSparse() 14sparse_data = sparse_model.createData() 15sparse_solver = odyn.SparseQP() 16sparse_solver.solve( 17 sparse_model, sparse_data, params, print_level=odyn.VerboseLevel.High) 18 19# We can create a dense models from the QP’s matrices and vectors 20other_model = odyn.DenseModel(model.Q, model.c, model.A, model.b, model.G, model.h) Listing 1: Example illustrating Odyn’s API. TABLE I: Odyn’s hyper-parameters for double precision. Parameter Value Description β, θ 0.85,0.950.85,0.95 neighborhood hyper-parameters σmin _min, σmax _max 0.1,0.90.1,0.9 minimum and maximum centering μmin _min 10−1610^-16 minimum barrier value θmin _min 0.10.1 minimum interpolation value θl _l, θu _u 0.3,0.850.3,0.85 min. and max. improvement thresholds θ−θ^-, θ+θ^+ 0.4,0.30.4,0.3 decrease and increase constants ρ 10−910^-9 proximal penalty parameters ρmin _min 10−910^-9 minimum penalty parameters δ 55 regularization decreasing factor η 10−210^-2 line search parameter V Odyn Application in Robotics and AI In this section, we briefly introduce three applications of Odyn: (i) an SQP solver built on top of the Odyn QP backend (OdynSQP), (i) a PyTorch module that embeds Odyn as an implicitly differentiable layer (ODYNLayer), and (i) a contact-physics engine that uses Odyn as its numerical backend (ODYNSim). The following subsections provide essential details on how Odyn is integrated into each of these components using state-of-the-art techniques. Further implementation details of these applications lie beyond the scope of this article and will be released as part of the open-source implementation. V-A OdynSQP: A SQP solver for MPC QP solvers are essential building blocks for solving nonlinear programs through sequential quadratic programming (SQP). The core idea of SQP is to transform a nonlinear program into a sequence of QPs, each solving a local linear-quadratic approximation of the cost and constraints. SQP methods are among the most powerful tools for nonlinear optimal control and estimation in robotics. Specifically, they can be applied to solve optimal control problems defined as follows min, _X,U ℓN(N)+∑k=0N−1ℓk(k,k) _N(x_N)+ _k=0^N-1 _k(x_k,u_k) (56) subject to 0=^0, _0= x_0, k+1=k(k,k), _k+1=f_k(x_k,u_k), k(k,k)=, _k(x_k,u_k)=0, k(k,k)≤, _k(x_k,u_k) 0, where =(,)∈⊆ℝnxx=(q,v) ^n_x denotes the system state, ∈ℝnuu ^n_u the control input, ℓN(N) _N(x_N) the terminal cost, and ℓk(k,k) _k(x_k,u_k) is the stage cost at node k. The system evolution is governed by the dynamics :×ℝnuf:X×R^n_u, while :×h:X×u and :×g:X×u encode equality and inequality path constraints. In robotics, equality constraints commonly enforce end-effector constraints [20, 51] or embed inverse dynamics [43, 18], whereas inequality constraints capture joint limits (e.g., [41]), friction cones, physical realism of dynamics parameters (e.g., [39]), and other safety or feasibility conditions (e.g., [11, 49]). This sequence of QP problems are sparse with block-banded structure. At each iteration, SQP forms a quadratic approximation of the Lagrangian, yielding the following QP: min _w 12δ(i)δ+∇(i)⊺δ 12 \,H(w_i)\, + (w_i) (57) subject to ¯(i)+∇¯(i)δ=, h(w_i)+∇ h(w_i) =0, ¯(i)+∇¯(i)δ≤, g(w_i)+∇ g(w_i) 0, where w stacks the state trajectory and control sequence. The matrices (i)H(w_i) and ∇(i) (w_i) are assembled from the derivatives of the stage and terminal costs at the current iterate iw_i, i.e., ℓ _x, ℓ _u, ℓ _xx, ℓ _xu, and ℓ _uu. Similarly, ∇¯(i)∇ h(w_i) is obtained by linearizing the dynamics (f_x, f_u) and stagewise equality constraints (h_x, h_u), and ∇¯(i)∇ g(w_i) is built from the linearized stagewise inequality constraints (g_x, g_u). OdynSQP is implemented using Crocoddyl [42] and Odyn’s sparse backend. In addition, we incorporate a Levenberg–Marquardt scheme to enhance convergence robustness, together with a non-monotone Armijo-type line search for step acceptance. V-B ODYNLayer: Differentiable optimization Differentiable optimization [1, 4, 37] has gained attention in differentiable physics simulation [47, 28] and policy learning [64]. This is because it enables gradient-based training through implicit optimization problems. Embedding a differentiable QP layer within a neural network allows joint optimization of perception and control parameters while preserving structural priors such as dynamics, constraints, and energy conservation. This approach underpins many model-based reinforcement learning methods [44], where differentiable rollouts integrate control subproblems (e.g., MPC or contact physics), yielding task-consistent gradients through structured decision processes instead of unstructured policy gradients. V-B1 Parametrized QPs A parameterized quadratic program QP() QP( θ) defines a family of QP solutions depending on parameters θ: min, _x,s 12⊺()+()⊺ 12\,x Q( θ)\,x+c( θ) x (58) subject to ()=(), ( θ)\,x=b( θ), ()+=(), ( θ)\,x+s=h( θ), ≥, 0, where its corresponding KKT point ⋆=(⋆,⋆,⋆,⋆)U =(x ,y ,s ,z ) depend on the parameters θ and its computed by satisfying the perturbed KKT residual ~(⋆;)= r(U ; )=0. For a scalar objective ℓ(⋆,) (U , θ), the derivative is computed implicitly through the optimality conditions rather than by unrolling the solver’s iterative steps. This is enabled by the implicit function theorem (IFT) (see [30]). The resulting approach is both numerically stable and computationally efficient, as it avoids the long backpropagation chains that often cause exploding or vanishing gradients when differentiating through iterative optimization procedures [1]. V-B2 Gradients via IFT Linearizing ~(⋆;,) r(U ; , θ) at ⋆U yields ∂~∂|=⋆d⋆+∂~∂d=. . ∂ r∂U |_U=U dU + ∂ r∂ θ\,d θ=0. (59) According to the IFT, iff ∂~/∂ r/ is non-singular, ∂⋆∂|=⋆ . ∂ θ |_U=U =−[∂~∂|=⋆]−1∂~⋆∂, =- [ . ∂ r |_U=U ]^-1 ∂ r ∂ θ, (60) with∂~⋆∂ ∂ r ∂ θ =[dd⋆+dd+(dd)⊺⋆+(dd)⊺⋆dd⋆−dddd⋆−dd]. = bmatrix dQd θ\,x + dcd θ+ ( dAd θ ) y + ( dGd θ ) z \\ dAd θ\,x - dbd θ\\ dGd θ\,x - dhd θ\\ 0 bmatrix. To avoid explicit Jacobian calculations, we introduce the adjoint λ as the solution of the linear system, [∂~∂|=⋆]⊺=−∂ℓ∂|=⋆, . [ ∂ r |_U=U ] λ=- . ∂ |_U=U , (61) while the hyper-gradient is defined using the solution ⋆ λ , ∂ℓ∂=(∗)⊺∂~∂. ∂ ∂ θ=( λ^*) ∂ r∂ θ. (62) Crucially, ∂~/∂ r/ coincides with the KKT matrix used in the forward Newton step (Eqs. 39 and 40), allowing the backward pass to efficiently compute the gradients of feasible QPs. In cases where the KKT matrix becomes ill-conditioned or nearly singular at convergence, a regularized variant of Eq. 62 is solved to ensure numerical stability. ODYNLayer is integrated into PyTorch. It supports both dense and sparse backends, multiple floating-point precisions, as well as condensed and full KKT factorizations. V-C ODYNSim: contact dynamics as QP Among all forces admissible under the unilateral and friction constraints, the realised contact force is the one that maximises energy dissipation. This dissipative nature of frictional contact is stated by maximum dissipation principle [45]. For rigid contacts with relative velocity v, this leads directly to a convex conic program. However, if we approximate the Coulomb friction cone for dry friction with a pyramidal model, the problem can be cast as a QP of the form min _ λ 12⊺()−1+⊺(()−+()) 12 λ (q)^-1 λ+ λ ( σ(q)^-+ (v)) (63) subject to ≥, λ 0, where ∈ℝ3nc λ ^3n_c describes the contact-force vector, ()=c()()−1c()⊺∈ℝ3nc×3nc (q)=J_c(q)M(q)^-1J_c(q) ^3n_c× 3n_c is the Delassus matrix with ()∈ℝnq×nq⪰M(q) ^n_q×n_q 0 and c()∈ℝ3nc×nqJ_c(q) ^3n_c×n_q as the joint-space mass matrix and contact Jacobian, respectively, ()−=c()−∈ℝ3nc σ(q)^-=J_c(q)v^- ^3n_c is the free velocity of the contact, ()∈ℝ3nc (v) ^3n_c is the de-Saxcé correction [13], and C defines a pyramidal approximation of the Coulomb friction cone. VI Results We begin by benchmarking the convergence behavior of our QP solver against state-of-the-art methods developed over the past decade using the Maros–Mészáros test set (Section VI-A). We then compare the warm-starting capabilities of Odyn with those of leading ALM-based solvers (Section VI-B). Next, we examine the behavior of Odyn on a set of representative degenerate problems (Section VI-D). Finally, we report results for OdynSQP, ODYNSim, and ODYNLayer in Sections VI-E, VI-G and VI-H. Together, these experiments demonstrate the suitability of Odyn for robotics and AI applications. All experiments were conducted on a MacBook Pro equipped with an M2 Pro processor. VI-A Hard QP problems: the Maros–Mészáros test set The Maros–Mészáros test set introduced in [38] comprises 138138 QP problems and serves as a standard benchmark for evaluating the performance of QP solvers. They are widely regarded as the gold standard for benchmarking QP solvers. Using the Maros–Mészáros test set, we evaluated Odyn against state-of-the-art sparse QP solvers, including ProxQP, PiQP, OsQP, Mosek, and Gurobi, under both medium- and high-accuracy configurations. We reported statistic of performance such as failure rate and performance profiles with absolute and relative tolerances set to 10−610^-6 or 10−910^-9, i.e., medium and high accuracy, respectively. VI-A1 Failure rates Table I summarizes the failure rates under both medium- and high-accuracy settings. The results indicate that Odyn is among the most competitive solvers and achieves the best performance at medium accuracy. Its strong convergence robustness can be attributed to its close connection with interior-point methods. Notably, IPMs (PiQP, MOSEK, and Gurobi) consistently outperform ALMs (ProxQP and OsQP) in this benchmark. Moreover, the slower convergence observed in Odyn under high-accuracy settings may be attributed to its lack of iterative refinement procedures, which are commonly employed in QP solvers. Note that most failures observed with Odyn were due either to exceeding the time limit (120120 seconds) or, in some cases, to the detection of primal or dual infeasibility. TABLE I: Failure rate (%) of different QP solvers on a subset of the Maros–Mészáros test set at medium and high accuracy. Accuracy Odyn ProxQP PiQP OsQP Mosek Gurobi Medium 9.7% 33.0% 10.7% 44.7% 14.6% 9.7% High 15.5% 35.9% 10.7% 49.5% 14.6% 12.6% VI-A2 Performance profiles We employed performance profiles [14] to evaluate and compare Odyn with state-of-the-art solvers. Profiles based on computation time and iteration count reported in Fig. 3 were obtained by running a set of solvers S on a collection of benchmark problems P. Each performance profile reports the fraction of problems solved by a specific solver, which is expressed as a function of the best observed performance for a given metric. Specifically, we defined the performance ratio of solver s∈s on problem p∈p as: rp,s=tp,smins∈tp,s,r_p,s= t_p,s _s t_p,s, (64) where tp,st_p,s represents computation time or number of iterations (our metrics). When a solver s failed to solve problem p, we set tp,st_p,s to an upper bound like the maximum allowed time to solve a problem or a maximum number of iterations. The performance of a solver s on a test set of problems P was defined as: fs(τ)=1||∑p∈ℐ≤τ(rp,s),f_s(τ)= 1|P| _p I_≤τ(r_p,s), (65) where ℐ≤τ(rp,s)=1I_≤τ(r_p,s)=1 if rp,s≤τr_p,s≤τ or 0 otherwise. Figure 3: (top) Iteration based performance profiles at high accuracy. (bottom) Solve time performance profiles at high accuracy. VI-B Warm-starting performance comparison We evaluated the warm-starting capabilities of Odyn against OsQP and ProxQP on a subset of the Maros–Mészáros test that are feasible and solvable by all the solvers. When doing so, we used the high-accuracy setting. Note that PiQP, Mosek, and Gurobi were excluded as interior-point solvers do not support warm-starting. We assessed the warm-starting capabilities by adopting an approach commonly used in the optimization (see [7]). Concretely, we first solved the original problem from a cold-start initialization using the QP problem =(,,,,,)QP=(Q,c,A,b,G,h). With this primal–dual solution, we recorded the number of iterations required for convergence, denoted by NcoldN_cold. We then constructed a perturbed problem ~=(~,~,~,~,~,~) QP=( Q, c, A, b, G, h) as described in [15]. For instance, the perturbed right-hand side ~ b is computed as ~=+δ(∘) b=b+δ(m η ), where δ>0δ>0 describes the different perturbation levels (0.001,0.01,0.1)(0.001,0.01,0.1), ∼Uniform([−1,1])m η ([-1,1])^m is a uniformly distributed random vector. Moreover, m is a random mask vector defined as: i=1,if i<min(0.1, 20/m)0,otherwise,m_i= cases1,&if ζ_i< (0.1,\;20/m)\\ 0,&otherwise cases, (66) where ∼Uniform([−1,1])m ζ ([-1,1])^m is another uniformly distributed random vector. Again, this strategy was used to perturb the entire QP problem. However, when doing so, we preserved the same sparsity structure and symmetric of the original matrices Q, A, and G. For each perturbation level δ, we recorded the number of iterations required to solve 1010 perturbed models ~i QP_i, with i∈1,…,10i∈\1,…,10\. This allowed us to compute the mean warm-to-cold ratio, defined as WCR=Nwarm/NcoldWCR=N_warm/N_cold, which quantifies the reduction in iteration count achieved through warm-starting relative to a cold start. Fig. 4 reports the warm-starting performance across different perturbation levels, showing Odyn’s superiority compared to OsQP and ProxQP. These results are reported for high-accuracy settings; however, a similar behavior was observed for medium accuracy. Additionally, these plots were generated using the standard benchmarking procedure described in Section VI-A2, which is more rigorous than directly plotting the WCR against the fraction of problems below a fixed threshold, as done in [15]. Figure 4: Warm-to-cold ratio performance profiles under high-accuracy settings, evaluated with perturbed QP data. The top and bottom plots correspond to perturbation levels δ=10−3δ=10^-3 and 10−110^-1, respectively. VI-C Dense backend We also benchmarked the dense backends of Odyn, ProxQP, and PiQP. Fig. 5 shows the average per-iteration computation time versus problem size. For each dimension, we generated random QP instances with m=p=n2m=p= n2 and repeated the solve 1010 times. It can be observed that Odyn is competitive. These results are critical for many real-time applications in robotics, including model predictive control, whole-body control, and contact simulation. They are also relevant for training differentiable QP layers in AI, where many applications are expected to involve small- to medium-scale, dense problem instances. Figure 5: Performance of the dense backend of Odyn, ProxQP, and PiQP. We generate random QP problems with dimensions ranging from 11 to 400400 and report the average computation time per iteration over 1010 trials. As expected, all solvers exhibit cubic complexity, whereas Odyn and PiQP demonstrate more deterministic timing behaviour. VI-D Degenerate QP problems We analyzed the robustness of Odyn on a set of degenerate problems. Iteration counts were benchmarked by treating all problems as sparse. The results, reported in Table IV, show that Odyn is competitive with other QP solvers. It also indicates that OsQP, a widely used solver in robotics due to its warm-starting capabilities, exhibits reduced numerical stability under the tested conditions. We selected simple QP problems so that they also serve an educational purpose. The corresponding test cases are described below. TABLE IV: Number of iterations required by each solver on the three degenerate QP problems. We highlight the three best solvers. Solver Redundant Ineq. LICQ Failure Multi. Solutions Odyn 5 10 6 ProxQP 12 30 10 PiQP 5 8 6 OsQP 75 † 200 Mosek 4 14 5 Gurobi 11 8 9 † declared dual infeasible on the LICQ failure problem. VI-D1 Redundant constraints We propose a simple QP that exhibits ill-conditioning due to redundant inequality constraints as follows minx1,x2 _x_1,x_2 12(x12+x22) 12 (x_1^2+x_2^2 ) (67) subject to x1≥1, x_1≥ 1, x1≥0. x_1≥ 0. In this QP problem, we observe that the second inequality constraint is redundant, yet it may still introduce numerical instability in many QP solvers. VI-D2 Linear dependent constraints We consider a QP problem that fails to meet the linear independence constraint qualification (LICQ) at convergence, while being ill-conditioned as its Q matrix and c vector are filled with small entries, i.e., minx1,x2 _x_1,x_2 12(10−10x12+10−6x22) 12\! (10^-10x_1^2+10^-6x_2^2 ) (68) +10−12x1x2+10−4x1−x2 +0^-12x_1x_2+0^-4x_1-x_2 subject to x1=0, x_1=0, x1≥0, x_1≥ 0, x2≥0. x_2≥ 0. In this QP problem, the optimal solution is =[00]x= bmatrix0&0 bmatrix, requiring the activation of all the constraints. This means that, when the inequality x1≥0x_1≥ 0 becomes active at the solution, there is linearly dependency with the equality constraint, thereby violating LICQ. VI-D3 Multiple solutions We also study a QP problem that admits multiple solutions due to a positive semidefinite Q matrix, i.e., minx1,x2 _x_1,x_2 12x12+x1 12x_1^2+x_1 (69) subject to 0x1+0x2≤0, 0x_1+0x_2≤ 0, 0≤x1≤3, 0≤ x_1≤ 3, 0≤x2≤3. 0≤ x_2≤ 3. In this case, the variable x1x_1 is free to take any value in the interval [1,3] bmatrix1,3 bmatrix, and all such values yield the same optimal cost. Additionally, the first constraint is a redundant inequality. This is because it is always active but has no effect on the feasible set. Moreover, its gradient is zero, which makes the KKT conditions ill-posed, since an active constraint with zero gradient violates the regularity conditions (e.g., LICQ) typically required for well defined Lagrange multipliers. Again, the results of solving these QP problems are summarized in Table IV. VI-E Optimizing agile maneuvers via OdynSQP We evaluated OdynSQP on a set of challenging trajectory optimization problems. Specifically, we considered three representative robotic systems with nonlinear dynamics subject to both inequality constraints (joint and control limits, friction-cone constraints) and equality constraints (end-effector targets and forward-dynamics consistency). In each of the examples, we included state and control regularization terms. Snapshots of the resulting optimal trajectories are reported in Fig. 6 and the feasibility evolution is depicted in Fig. 7. OdynSQP successfully optimized highly dynamic maneuvers, such as backflips on the G1 humanoid robot, while enforcing joint position, velocity, and torque limits together with contact wrench-cone constraints. It also enabled loco-manipulation, including the “chicken-head” task on the Unitree B1 quadruped equipped with the Z1 manipulator. In addition to joint and friction-cone constraints, we imposed bounding-box constraints on the end-effector, as illustrated in Fig. 6. Furthermore, OdynSQP handled aerial manipulation scenarios requiring end-effector target constraints and thruster limits. From Fig. 7, we observe that OdynSQP decreased the constraint feasibility to high accuracy (i.e., 10−9)( i.e., 10^-9) while allowing temporary increases in feasibility through its non-monotone line-search globalization strategy exhibiting strong convergence capabilities. These experiments demonstrated the ability of OdynSQP to efficiently solve highly nonlinear optimal control problems across diverse robotic platforms. Figure 6: Snapshots of optimal trajectories computed with OdynSQP. (Top) G1 humanoid performing a backflip under joint position, velocity, and torque limits, together with contact wrench-cone constraints. (Bottom-left) Borinot aerial robot reaching a target end-effector position (green sphere) while satisfying actuator torque limits. (Bottom-right) B1–Z1 executing a chicken-head stabilization task during lateral locomotion while enforcing end-effector and friction-cone constraints; the robot’s hand is constrained to remain within the green bounding box. Figure 7: Feasibility evolution on challenging constrained optimal control problems computed by OdynSQP. VI-F Model predictive control via OdynSQP We evaluated OdynSQP in a model predictive control setting on two Unitree robotic platforms: the B1 quadruped and the Z1 manipulator. These experiments assess the tracking performance and real-time capabilities of OdynSQP in MPC loops (see Fig. 8). The first experiment considers the B1 quadruped performing a quasi-static walking-in-place motion under contact constraints. Our MPC formulation enforced friction-cone constraints together with state and input constraints, while tracking a desired center-of-mass and angular momentum trajectory, as shown in Fig. 9. The MPC horizon was defined with nine nodes discretized at 33ms33\,ms, corresponding to a prediction horizon of 0.3s0.3\,s, and the control loop runs at 30Hz30\,Hz. Despite frequent changes in the active set induced by contact transitions, OdynSQP maintained stable convergence across successive MPC iterations, resulting in smooth and consistent motion generation. Figure 8: Snapshots of MPC trajectories generated by OdynSQP. (Top) B1 quadruped executing a quasi-static walking-in-place motion under contact and friction-cone constraints, while enforcing state and input constraints. (Bottom) Z1 manipulator tracking a three-dimensional ellipsoidal end-effector trajectory (blue dashed curve) under state and input constraints. Despite changing active sets and high-frequency updates, OdynSQP achieved accurate tracking and stable behavior across successive MPC iterations.To watch the video, click the picture or see https://w.youtube.com/watch?v=fJqjgdT2De0. Figure 9: Model predictive control tracking performance on the Unitree B1 quadruped using OdynSQP. The plots show the tracking of the center-of-mass position (bottom) and angular momentum (top) along the y-axis during a quasi-static walking-in-place motion. Figure 10: Model predictive control tracking performance on the Unitree Z1 manipulator using OdynSQP. The plots show the tracking of the end-effector position along selected axes while following a three-dimensional ellipsoidal trajectory. Solver statistics further highlight the efficiency of OdynSQP. For the Z1 manipulator, the average solve time was 13.5ms13.5\,ms, with a maximum of 19.6ms19.6\,ms, remaining strictly below the 20ms20\,ms control period. For the B1 quadruped, the average solve time was 19.7ms19.7\,ms, with a maximum of 29.4ms29.4\,ms, also below the 33.3ms33.3\,ms control period. In both cases, the solver was limited to a single SQP iteration per MPC step, demonstrating the capability of OdynSQP to operate reliably under real-time constraints. These results demonstrated that OdynSQP operates reliably across different MPC regimes, ranging from locomotion to manipulation, and provides an efficient backend for solving sequences of structured QPs in real-time robotics applications. VI-G Fast resolution of contact physics via ODYNSim We evaluated Odyn, ProxQP, and OsQP in the context of contact-physics simulation. The benchmark consists of a pushed-cube scenario in which contact forces are computed by solving a quadratic program at each simulation step. This setup provides a controlled environment to compare solver robustness, convergence behavior, and computational efficiency. A single cube (edge length 0.20m0.20\,m, mass 15kg15\,kg) rests on a fixed plane, an external horizontal force was applied throughout the simulation, causing it to slide across the surface. The force magnitude was set to 2μmg≈147.15N2μ mg≈ 147.15\,N, and its direction alternated every 0.5s0.5\,s between the world x- and y-axes. Within each interval, the force was ramped linearly over the first 0.3s0.3\,s and then held constant; the sign of the force was deterministically flipped at each interval. At every time step, the simulator enforced non-penetration and Coulomb friction constraints as formulated in Eq. 63. Contacts were modeled using pyramidal friction cones with coefficient μ=0.5μ=0.5. The time step was set to dt=0.004sdt=0.004\,s, and this QP problem was solved at each step with ϵa=10−8 _a=10^-8 and ϵr=0 _r=0. Figure 11: Pushed cube contact simulation resolved via QP-based contact dynamics. (Top) illustrate the cube’s interaction with the plane under unilateral frictional contact constraints at different time instants. (Bottom) compares the number of solver iterations required per time step for Odyn, ProxQP, and OsQP. As shown in Fig. 11, Odyn consistently required fewer iterations per time step than ProxQP and OsQP. In particular, OsQP exhibited significantly higher iteration counts, reflecting the limited efficiency of first-order operator-splitting methods in exploiting warm starts at high accuracy. This behavior is evident even when successive QPs differ only marginally due to incremental state and contact updates. When transitions occur between sticking and sliding contact, the optimal solution changes more substantially. Even in these cases, Odyn converged in fewer iterations than ProxQP. Overall, these results indicate that Odyn is well suited to real-time contact-dynamics simulation, where sequences of related QPs must be solved efficiently. VI-H Learning to play Sudoku via ODYNLayer ODYNLayer remains differentiable, even with degenerate problems, thanks to its smoothed NCP regularization and proximal primal–dual structure. This brings a more principled way to develop differentiable without relying on sub-gradients. To demonstrate ODYNLayer’s capabilities, we trained a QP optimization to learn playing Sudoku as similarly done in [1]. (a) Double Precision (float64) (b) Single Precision (float32) (c) Ground Truth vs. Model Prediction Figure 12: Effect of ODYNLayer precision on convergence and invalid-Sudoku rate. (a) With double precision (float64), it exhibited stable convergence with a rapid decline in invalid Sudoku solutions. (b) With single precision (float32), it showed slightly slower convergence in terms of invalid Sudoku count, a degradation that is attributed to Odyn’s single-precision convergence tolerance and randomized-parameter initialization. (c) Example reconstruction comparing ground truth with the model prediction; digits shown in blue denote the fixed givens supplied to ODYNLayer. These results highlight how numerical precision impacts both convergence speed and true Sudoku constraint satisfaction (invalid Sudoku count). True Sudoku constraint satisfaction implies that the learned linear constraint basis accurately captures the underlying Sudoku structure. Sudoku is a logic-based puzzle played on a N×N× N grid divided into n×n× n subgrids (“blocks”), with n=3n=3 and N=n2=9N=n^2=9 typically. Sudoku is inherently a binary problem as it encodes digit assignments using binary variables xi,j,k∈0,1x_i,j,k∈\0,1\ indicating whether cell (i,j)(i,j) that is located inside a n×n× n sub-block b, contains digit k for i,j,k∈1,…,Ni,j,k∈\1,…,N\. Theses rules are captured by equality constraints enforcing that every row, column, and n×n× n sub-block contains each digit exactly once: ∑i=1Nxi,j,k=1, _i=1^Nx_i,j,k=1, ∀i,k, ∀ i,k, ∑j=1Nxi,j,k=1, _j=1^Nx_i,j,k=1, ∀j,k, ∀ j,k, (70) ∑(i,j)∈bxi,j,k=1, _(i,j)∈ bx_i,j,k=1, ∀b,k, ∀ b,k, with pre-filled cells (“givens”) imposed by fixing xi,j,k0=1x_i,j,k_0=1. VI-H1 Formulation To integrate this structure into a differentiable pipeline, we relaxed the integer nature of the decision variables by considering bounded real numbers (i.e.. 0≤xi,j,k≤10≤ x_i,j,k≤ 1), transforming it into a LP. This relaxation is exact, since the constraint matrix defined by Eq. 70 is totally unimodular. This in turn implies that all extreme points of the feasible polytope are integer-valued [8]. Thus, solving this LP recovers valid Sudoku solutions without enforcing integrality explicitly. Specifically, we proposed to formulate a strictly feasible, regularized relaxation of this LP via ODYNLayer, i.e., min∈ℝn6 _x ^n^6 ϵ2‖22+⊺ ε2\,\|x\|_2^2+c x (71) subject to ()=(), ( θ)\,x=b( θ), ≤, 0 1, where x encodes the relaxed digit-assignment variables of the Sudoku instance and ϵ=0.1ε=0.1 is a small Tikhonov regularization parameter [60] to follow OptNet’s problem setup [1]. Givens are incorporated as soft constraints through a linear penalty vector =−k0c=-x_k_0, which biases the optimizer toward satisfying the given entries while preserving differentiability. To enforce strict feasibility, we set ()b( θ) to be in the range of ()A( θ). This enabled smooth implicit differentiation through the QP, while maintaining the exact combinatorial structure of the Sudoku constraints. VI-H2 Training capabilities We evaluated ODYNLayer on a diverse collection of 4×44× 4 Sudoku instances. For training, we generated 500500 randomly sampled Sudoku instances and reserved a disjoint set of 10001000 additional instances for testing. This setup yields a substantially smaller training set than that used by OptNet (90009000 Sudoku instances). However, we employed additional training epochs to ensure reliable convergence under this more limited data regime. The results are summarized in Fig. 12. Training proceeded by minimizing the discrepancy between the predicted digits ^i,j() x_i,j( θ) and the ground-truth values i,j∗x^*_i,j for each cell (i,j)(i,j) across batches of size B, namely, ℒ()=1BN2∑b=1B∑i=1N∑j=1N‖^i,j(b)()−i,j∗(b)‖22.L( θ)= 1BN^2 _b=1^B _i=1^N _j=1^N \| x^(b)_i,j( θ)-x^*(b)_i,j \|_2^2. (72) This loss function served as a surrogate objective, guiding the identification of a linear-constraint basis that faithfully encodes the Sudoku rules defined in Eq. 70. VI-H3 Effect of multiple setups We analyzed the capabilities of ODYNLayer across different numerical configurations, considering both sparse and dense backends as well as full-KKT and condensed-KKT factorizations. The results of this ablation study are reported in Table V. ODYNLayer produces zero invalid Sudoku predictions in all configurations except for the condensed-KKT factorization in single precision. This behavior is consistent with the well-known numerical instability associated with condensation-based approaches. Furthermore, for the full-KKT factorization, the faster convergence rate observed with single precision is explained by the need to relax the solver stopping tolerances—set to 10−910^-9 in double precision and 10−510^-5 in single precision—when operating at reduced numerical accuracy. TABLE V: Sudoku ablation study with ODYNLayer: We indicate the average solver iterations per epoch (iters) during training and the number of invalid Sudoku predictions (#errors) in the out-of-distribution test-set (10001000 Sudoku instances). The full-KKT factorization setup is indicated by “Full” and the condensed by “Cond.”. Double (float64) Float (float32) #errors iters (±σ±σ) #errors iters (±σ±σ) Full Dense 0 11.38(±0.31)11.38\;(± 0.31) 0 8.91(±0.30)8.91\;(± 0.30) Sparse 0 11.16(±0.29)11.16\;(± 0.29) 0 8.98(±0.30)8.98\;(± 0.30) Cond. Dense 0 11.36(±0.28)11.36\;(± 0.28) - (max iters) - (max iters) Sparse 0 11.14(±0.28)11.14\;(± 0.28) - (max iters) - (max iters) VII Conclusion This paper presented Odyn, a novel all-shifted non-interior-point QP solver based on shifted NCP functions and proximal primal–dual regularization. By relaxing strict interior-feasibility requirements while preserving path-following behavior, Odyn provides a principled alternative to classical interior-point methods. The resulting framework enables robust treatment of degenerate, ill-conditioned, and warm-started problems that frequently arise in robotics and AI workloads. Extensive numerical experiments on standard benchmarks demonstrate that Odyn achieves competitive convergence performance on par with modern interior-point methods while exhibiting superior warm-start capabilities compared to state-of-the-art augmented Lagrangian and operator-splitting approaches. Beyond standalone QP benchmarks, Odyn was validated in representative robotic and AI applications, including sequential quadratic programming, contact-dynamics simulation, and differentiable optimization layers, illustrating its versatility across robotics, physics simulation, and learning pipelines. In summary, Odyn blends the advantages of primal–dual interior-point methods and generalized augmented Lagrangian techniques, offering a powerful and flexible foundation for the development of the next generation of optimization solvers for robotics, AI, and numerical optimization. A. All-shifted complementarity constraints From the KKT conditions of the perturbed NCP Lagrangian (Eq. 32), namely ∇~= _ ξ P=0 and ∇~= _w P=0, we obtain −μ−1++ρn2(−E)+ρn2(−E)=.-μ ξ^-1+w+ _n2(s-s_E)+ _n2(w-w_E)=0. (73) The condition ϕ(,;μ)= φ(s,w;μ)=0 implies =s= ξ. Substituting into Eq. 73 yields the all-shifted complementary constraints: ∘=μ+ρn2∘(E−+E−).s = 1+ _n2\,s (s_E-s+w_E-w). (74) Finally, non-negativity of s and w follows directly from ϕ(,;μ)= φ(s,w;μ)=0. References [1] B. Amos and J. Z. Kolter (2017) OptNet: Differentiable Optimization as a Layer in Neural Networks. In Int. Conf. Mach. Learn., Cited by: §V-B1, §V-B, §VI-H1, §VI-H. [2] R. Andreani, E. G. Birgin, J. M. Martínez, and M. L. Schuverdt (2008) On Augmented Lagrangian Methods with General Lower-Level Constraints. SIAM J. Optim. 18. Cited by: §I. [3] A. Bambade, S. El-Kazdadi, A. Taylor, and J. Carpentier (2022) ProxQP: Yet another quadratic programming solver for robotics and beyond. In Rob.: Sci. Sys. (RSS), Cited by: §I, §I-A, TABLE I, §I-B, §IV-F, §IV-G. [4] A. Bambade, F. Schramm, A. Taylor, and J. Carpentier (2024) Leveraging augmented-Lagrangian techniques for differentiating over infeasible quadratic programs in machine learning. In Int. Conf. on Learn. Repr., Cited by: §I-A, §V-B. [5] G. Banjac, B. Stellato, N. Moehle, P. Goulart, A. Bemporad, and S. Boyd (2017) Embedded code generation using the OSQP solver. In IEEE Int. Conf. on Dec. Cntr. (CDC), Cited by: §I-A. [6] G. Banjac, P. Goulart, B. Stellato, and S. Boyd (2019) Infeasibility Detection in the Alternating Direction Method of Multipliers for Convex Optimization. J. Optim. Theory Appl 183. Cited by: §IV-I. [7] H. Y. Benson and D. F. Shanno (2007) An exact primal–dual penalty method approach to warmstarting interior-point methods for linear programming. Computational Optimization and Applications 38. Cited by: §I-B, §VI-B. [8] C. Blair (1988) Theory of Linear and Integer Programming (Alexander Schrijver). SIAM Review 30. Cited by: §VI-H1. [9] J. V. Burke and S. Xu (1998) A non-interior predictor-corrector path-following method for LCP. In Reformulation: Nonsmooth, Piecewise Smooth, Semismooth and Smoothing Methods, Cited by: §I-C. [10] J. V. Burke and S. Xu (1998) The global linear convergence of a noninterior path-following algorithm for linear complementarity problems. Mathematics of Operations Research 23. Cited by: §I-C. [11] S. Caron, Q. Pham, and Y. Nakamura (2015) Stability of surface contacts for humanoid robots: Closed-form formulae of the Contact Wrench Cone for rectangular support areas. In IEEE Int. Conf. Rob. Autom. (ICRA), Cited by: §V-A. [12] A. R. Conn, N. I. M. Gould, and P. L. Toint (2010) Lancelot: A Fortran Package for Large-Scale Nonlinear Optimization (Release A). Springer Publishing Company, Incorporated. Cited by: §I. [13] G. De Saxcé and Z.-Q. Feng (1998) The bipotential method: A constructive approach to design the complete contact law with friction and improved numerical algorithms. Mathematical and Computer Modelling 28. Cited by: §V-C. [14] E. D. Dolan and J. J. Moré (2002) Benchmarking optimization software with performance profiles. Mathematical programming 91. Cited by: §VI-A2. [15] A. Engau, M. F. Anjos, and A. Vannelli (2009) A primal-dual slack approach to warmstarting interior-point methods for linear programming. In Operations Research and Cyber-Infrastructure, Cited by: §I-B, §VI-B, §VI-B. [16] S. Engelke and C. Kanzow (2002) Predictor-corrector smoothing methods for linear programs with a more flexible update of the smoothing parameter. Computational Optimization and Applications 23. Cited by: §I-C. [17] H. J. Ferreau, C. Kirches, A. Potschka, H. G. Bock, and M. Diehl (2014) qpOASES: a parametric active-set algorithm for quadratic programming. Mathematical Programming Computation 6. Cited by: §I-A, TABLE I. [18] H. Ferrolho, V. Ivan, W. Merkt, I. Havoutis, and S. Vijayakumar (2021) Inverse Dynamics vs. Forward Dynamics in Direct Transcription Formulations for Trajectory Optimization. In IEEE Int. Conf. Rob. Autom. (ICRA), Cited by: §V-A. [19] G. Frison and M. Diehl (2020) HPIPM: a high-performance quadratic programming framework for model predictive control. In Int. Fed. of Aut. Cntrl. (IFAC), Cited by: §I-A. [20] M. Giftthaler and J. Buchli (2017) A projection approach to equality constrained iterative linear quadratic optimal control. In IEEE Int. Conf. Hum. Rob. (ICHR), Cited by: §V-A. [21] P. E. Gill and M. Zhang (2024) A projected-search interior-point method for nonlinearly constrained optimization. Computational Optimization and Applications 88. Cited by: §I-B. [22] D. Goldfarb and A. Idnani (1983) A numerically stable dual method for solving strictly convex quadratic programs. Mathematical Programming 27 (1), p. 1–33. Cited by: §I-A. [23] G. H. Golub and C. F. V. Loan (2013) Matrix computations. fourth edition, The Johns Hopkins University Press. Cited by: §I, §IV-C. [24] J. Gondzio and A. Grothey (2008) A New Unblocking Technique to Warmstart Interior Point Methods Based on Sensitivity Analysis. SIAM J. Optim. 19. Cited by: §I-B. [25] G. Guennebaud, B. Jacob, et al. (2010) Eigen. Note: https://libeigen.gitlab.io Cited by: §IV-J. [26] Gurobi Optimization, LLC (2024) Gurobi Optimizer Reference Manual. External Links: Link Cited by: §I, TABLE I. [27] Z. Huang, J. Han, and Z. Chen (2003) Predictor-corrector smoothing Newton method, based on a new smoothing function, for solving the nonlinear complementarity problem with a P 0 function. J. Optim. Theory Appl 117. Cited by: §I-B. [28] K. M. Jatavallabhula, M. Macklin, F. Golemo, V. Voleti, L. Petrini, M. Weiss, B. Considine, J. Parent-Levesque, K. Xie, K. Erleben, L. Paull, F. Shkurti, D. Nowrouzezahrai, and S. Fidler (2021) gradSim: Differentiable simulation for system identification and visuomotor control. In Int. Conf. on Learn. Repr., Cited by: §V-B. [29] C. Kanzow (1996) Some noninterior continuation methods for linear complementarity problems. SIAM J. Mat. Analy. App. 17. Cited by: §I-C, §I-B, §IV-E. [30] S. G. Krantz and H. R. Parks (2002) The Implicit Function Theorem: History, Theory, and Applications. Birkhäuser. Cited by: §V-B1. [31] J. Lee and S. Leyffer (2011) Mixed integer nonlinear programming. Vol. 154, Springer Science & Business Media. Cited by: §I. [32] D. Liao-McPherson and I. Kolmanovsky (2019) FBstab: A stabilized semismooth quadratic programming algorithm with applications in model predictive control. arXiv preprint arXiv:1901.04046. Cited by: §IV-I. [33] D. Liao-McPherson and I. Kolmanovsky (2020) FBstab: a proximally stabilized semismooth algorithm for convex quadratic programming. Automatica 113. Cited by: §I-C, TABLE I. [34] K. Lin and T. Ohtsuka (2025) A non-interior-point continuation method for the optimal control problem with equilibrium constraints. Automatica 171. Cited by: §I-C. [35] P. Ma, J. Chen, C. Huang, and C. Ko (2018) Discovery of new complementarity functions for NCP and SOCCP. Computational and Applied Mathematics 37. Cited by: §I-C. [36] M. Macklin, K. Erleben, M. Müller, N. Chentanez, S. Jeschke, and V. Makoviychuk (2019) Non-smooth newton methods for deformable multi-body dynamics. ACM Trans. Graph. (TOG) 38. Cited by: §I-C. [37] C. W. Magoon, F. Yang, N. Aigerman, and S. Z. Kovalsky (2025) Differentiation Through Black-Box Quadratic Programming Solvers. External Links: 2410.06324, Link Cited by: §V-B. [38] I. Maros and C. Mészáros (1999) A repository of convex quadratic programming problems. Optimization Methods and Software 11. Cited by: §I-A, §VI-A. [39] S. Martinez, R. Griffin, and C. Mastalli (2025) Multi-Contact Inertial Parameters Estimation and Localization in Legged Robots. IEEE Robot. Automat. Lett. (RA-L). Cited by: §I, §V-A. [40] S. Martinez, S. Tonneau, and C. Mastalli (2026) System Identification under Constraints and Disturbance: A Bayesian Estimation Approach. Cited by: §I. [41] C. Mastalli, J. Marti-Saumell, W. Merkt, J. Sola, N. Mansard, and S. Vijayakumar (2022) A Feasibility-Driven Approach to Control-Limited DDP. Autom. Robots. 46. Cited by: §V-A. [42] C. Mastalli, R. Budhiraja, W. Merkt, G. Saurel, B. Hammoud, M. Naveau, and N. Mansard (2020) Crocoddyl: An efficient and versatile framework for multi-contact optimal control. In IEEE Int. Conf. Rob. Autom. (ICRA), Cited by: §V-A. [43] C. Mastalli, S. P. Chhatoi, T. Corbéres, S. Tonneau, and S. Vijayakumar (2023) Inverse-Dynamics MPC via Nullspace Resolution. IEEE Trans. Robot. (T-RO) 39. Cited by: §V-A. [44] T. M. Moerland, J. Broekens, A. Plaat, and C. M. Jonker (2023) Model-based Reinforcement Learning: A Survey. Found. Trends Mach. Learn. 16. Cited by: §V-B. [45] J. J. Moreau and P. P.D. (1988) Nonsmooth mechanics and applications. Springer-Verlag. Cited by: §V-C. [46] MOSEK ApS (2024) MOSEK Optimizer API for Python, Version 10.0. External Links: Link Cited by: §I, TABLE I. [47] Newton: GPU-accelerated physics simulation for robotics, and simulation research Newton a Series of LF Projects, LLC. External Links: Link Cited by: §V-B. [48] J. Nocedal and S.J. Wright (2006) Numerical Optimization. second edition, Springer, New York, USA. Cited by: §I, §IV-E. [49] D. Ordoñez-Apraez, G. Turrisi, V. Kostic, M. Martin, A. Agudo, F. Moreno-Noguer, M. Pontil, C. Semini, and C. Mastalli (2025) Morphological symmetries in robotics. The Int. J. of Rob. Res. (IJRR) 44. Cited by: §V-A. [50] A. G. Pandala, Y. Ding, and H. Park (2019) qpSWIFT: A Real-Time Sparse Quadratic Program Solver for Robotic Applications. IEEE Robot. Automat. Lett. (RA-L). Cited by: §I-A, TABLE I. [51] M. Parilli, S. Martinez, and C. Mastalli (2025) Endpoint-Explicit Differential Dynamic Programming via Exact Resolution. In IEEE Int. Conf. Rob. Autom. (ICRA), Cited by: §V-A. [52] L. Qi and D. Sun (1993) A nonsmooth version of newton’s method. Mathematical Programming 58 (1), p. 353–367. Cited by: §I-C. [53] R. T. Rockafellar (1976) Monotone operators and the proximal point algorithm. SIAM Journal on Control and Optimization 14 (5), p. 877–898. Cited by: §I-C. [54] M. Schubiger, G. Banjac, and J. Lygeros (2020) GPU acceleration of ADMM for large-scale quadratic programming. Journal of Parallel and Distributed Computing 144. Cited by: §I-A. [55] R. Schwan, Y. Jiang, D. Kuhn, and C. N. Jones (2023) PIQP: A proximal interior-point quadratic programming solver. In IEEE Int. Conf. on Dec. Cntr. (CDC), Cited by: §I, §I-A, TABLE I, §IV-E, §IV-F, §IV-H. [56] A. Skajaa, E. D. Andersen, and Y. Ye (2013) Warmstarting the homogeneous and self-dual interior point method for linear and conic quadratic problems. Mathematical Programming Computation 5. Cited by: §I-B. [57] Stack-of-Tasks Eiquadprog: a c++ implementation of the goldfarb–idnani quadratic programming algorithm. Note: Online. Available: https://github.com/stack-of-tasks/eiquadprog Cited by: §I-A, TABLE I. [58] B. Stellato, G. Banjac, P. Goulart, A. Bemporad, and S. Boyd (2020) OSQP: an operator splitting solver for quadratic programs. Mathematical Programming Computation 12. Cited by: §I, §I-A, TABLE I. [59] D. Sun (1999) A regularization Newton method for solving nonlinear complementarity problems. Applied Mathematics and Optimization 40. Cited by: §I-C. [60] A. N. Tikhonov and V. Y. Arsenin (1977) Solutions of Ill-Posed Problems. Winston, Washington, DC. Cited by: §VI-H1. [61] E. Todorov (2010) Implicit nonlinear complementarity: A new approach to contact dynamics. In IEEE Int. Conf. Rob. Autom. (ICRA), Cited by: §I-C. [62] V. Tsounis, R. Grandia, and M. Bächer (2025) https://arxiv.org/abs/2504.19771On Solving the Dynamics of Constrained Rigid Multi-Body Systems with Kinematic Loops. External Links: 2504.19771 Cited by: §I. [63] A. Wächter and L. T. Biegler (2006) On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. Mathematical Programming 106, p. 25–57. Cited by: §I. [64] J. Xu, V. Makoviychuk, Y. Narang, F. Ramos, W. Matusik, A. Garg, and M. Macklin (2021) Accelerated Policy Learning with Parallel Differentiable Simulation. In Int. Conf. on Learn. Repr., Cited by: §V-B. [65] H. Zhang and W. W. Hager (2004) A Nonmonotone Line Search Technique and Its Application to Unconstrained Optimization. SIAM J. Optim. 14. Cited by: §IV-D. [66] R. Zhang, X. Liu, and Y. Dai (2023) IPRQP: a primal-dual interior-point relaxation algorithm for convex quadratic programming. Journal of Global Optimization 87. Cited by: §I-C. [67] R. Zhang, X. Liu, and Y. Dai (2024) IPRSDP: a primal-dual interior-point relaxation algorithm for semidefinite programming. Computational Optimization and Applications 88. Cited by: §I-C.