Differentiation Through Black-Box Quadratic Programming Solvers

## Abstract

Abstract In recent years, many deep learning approaches have incorporated layers that solve optimization problems ( e.g. , linear, quadratic, and semidefinite programs). Integrating these optimization problems as differentiable layers requires computing the derivatives of the optimization problem’s solution with respect to its objective and constraints. This has so far prevented the use of state-of-the-art black-box numerical solvers within neural networks, as they lack a differentiable interface. To address this issue for one of the most common convex optimization problems – quadratic programming (QP) – we introduce dQP , a modular framework that enables plug-and-play differentiation for any QP solver, allowing seamless integration into neural networks and bi-level optimization tasks. Our solution is based on the core theoretical insight that knowledge of the active constraint set at the QP optimum allows for explicit differentiation. This insight reveals a unique relationship between the computation of the solution and its derivative, enabling efficient differentiation of any solver, that only requires the primal solution. Our implementation, which will be made publicly available, interfaces with an existing framework that supports over 15 state-of-the-art QP solvers, providing each with a fully differentiable backbone for immediate use as a differentiable layer in learning setups. To demonstrate the scalability and effectiveness of dQP, we evaluate it on a large benchmark dataset of QPs with varying structures. We compare dQP with existing differentiable QP methods, demonstrating its advantages across a range of problems, from challenging small and dense problems to large-scale sparse ones, including a novel bi-level geometry optimization problem.

## 1 Introduction

Computational methods rely heavily on solving optimization problems, i.e., finding an optimum of a given function, under some given constraints on the solution.
Optimization is arguably the most popular method to approach computational problems that do not admit a closed-form solution.
This in turn has led to the development of both open-source and commercial numerical solvers specialized for different classes of optimization, especially constrained convex optimization (Wright, ; Boyd & Vandenberghe, ).
It is, thus, quite enticing to incorporate optimization as a “layer” within machine learning architectures, e.g., where a neural network’s intermediate output defines the optimization problem, and the solution of that optimization problem is taken as the final output of the neural network (Amos & Kolter, ; Agrawal et al., ; Blondel & Roulet, ). This approach has proven successful on many practical tasks including image classification (Amos et al., ), optimal transport (Rezende & Racanière, ; Richter-Powell et al., ), zero-sum games (Ling et al., ), tessellation (Chen et al., ), control (Amos et al., ; de Avila Belbute-Peres et al., ; Ding et al., ), decision-making (Tan et al., ), robotics (Holmes et al., ), biology (Zhang et al., ), and natural language processing (Thayaparan et al., ).

In general, training a neural network requires the ability to *backpropagate* gradients to optimize the network’s weights and biases. Hence, in case the network includes an optimization layer as described above, one needs to have a way to *differentiate* that layer, i.e., compute the gradients of the solution of the optimization problem with respect to the parameters of the optimization problem itself. Gradients can be obtained through optimality conditions, which provide a characterization that allows for the application of the implicit function theorem (Krantz & Parks, ) and implicit differentiation. However, this approach requires the dual solution, which many optimizers do not provide, and additionally yields a linear system for the gradients that can be costly to invert for large problems. As a result, many existing differentiable methods tightly couple the differentiation to a custom tailor-made solver that outputs a dual solution, allows information from the solution algorithm to be re-used, or uses GPU acceleration.

Figure: Figure 1: Comparison over QP problems of increasing size. Ours (using Gurobi) outperforms OptNet and QPLayer in terms of time and accuracy as problem size increases. OptNet and QPLayer become intractable for larger problems.
Refer to caption: x1.png

The existing tight coupling between neural architecture and optimizer severely limits the applicability of neural optimization methods: in general, solving optimization problems is a hard, challenging task, requiring state-of-the-art solvers such as Gurobi (Gurobi Optimization, LLC, ) and MOSEK (Andersen & Andersen, ) which have been developed through years of commercial and academic research. These solvers provide the capability to efficiently and reliably handle problems at a scale that non-optimized implementations cannot achieve. More so, even having one of these solvers at hand is not enough, since none of them is a “catch-all” solution.
Instead, it is necessary to choose and swap between specific solvers for specific structures of optimization problems that may emerge in different learning tasks. To solve this issue and obtain a general, efficient way to interface between general solvers and neural networks via backpropgation, one needs to devise a “bridge” to differentiate through these blackbox solvers.

In this paper, we focus on completely removing the above limitation to one of the most canonical and important convex optimization problems, quadratic programing (QP), which minimize a quadratic objective under linear inequality constraints.
Our framework, which we dub dQP (as in differential notation), augments any QP solver into a differentiable one and seamlessly integrates it as a differentiable layer.
dQP stems from our main novel theoretical observation: the gradients of the optimal point of a QP can be obtained explicitly from a primal solution provided by the optimizer, and the active set of constraints, which can be deduced from the solution.

Figure: Figure 2: dQP solves a bi-level optimization problem to compute the bijective planar embedding of a large-scale ant mesh (15k vertices).
Refer to caption: x2.png

We draw this conclusion by leveraging classic sensitivity theory, and by clarifying the role of the active set which has otherwise appeared in recent work, but without emphasis and not in this form. Notably, our explicit perspective recasts the traditional implicit differentiation approach into an *explicit* method, which provides a straightforward pathway to complete solver modularity. Namely, we avoid the costly linear solve for the necessary gradients by showing the full gradients can be expressed solely through the (much smaller) active set. Furthermore, we show that this reduced system can also be used to solve for the active dual variables, if not provided by a solver. This enables us to implement dQP on top of a minimal open-source interface (Caron et al., ), which provides direct access to over 15 free and commercial QP solvers and easily supports the integration of additional solvers.

Using the modularity of our method, we conducted a comprehensive evaluation on a diverse benchmark dataset comprising over 2,000 QPs, comparing dQP’s performance against existing differentiable QP methods. As highlighted in Figure 1, dQP demonstrates a significant advantage in structured QPs when paired with state-of-the-art sparse QP solvers. We show the superiority of our method on tasks such as sparse projections, as well as in a novel geometric bi-level optimization experiment that was intractable for previous methods, whereas our method excels.

To summarize, our contributions are:

- 1.
We prove that QPs can be explicitly differentiated using only the primal solution via a locally-equivalent linear system.
- 2.
Building on this, we devise and implement a fully modular differentiable layer compatible with any QP solver, allowing for plug-and-play flexibility where users can easily select the best solver for their specific task. Our open-source implementation will be made publicly available.
- 3.
We demonstrate state-of-the-art performance in solving and differentiating large-scale, sparse QPs using various solvers across a series of extensive experiments.

## 2 Related Works

Implicit Layers.
Optimization layers are an example of recently introduced implicit layers that leverage implicit differentiation to compute gradients of solution mappings without requiring closed-form expressions (Duvenaud et al., ). This category also includes deep equilibrium models, which represent fixed-point mappings and can be viewed as networks of infinite depth (Bai et al., ; Kawaguchi, ; El Ghaoui et al., ; Gurumurthy et al., ; Winston & Kolter, ; Bai et al., ). Extending beyond algebraic equations, similar techniques are applied in neural ordinary differential equations using the adjoint state method from parametric control for partial differential equations (Lions, ; Xue et al., ; Beatson et al., ; Chen et al., ). Implicit differentiation is also used in bi-level programming where optimization problems are nested in one-another (Colson et al., ; Kunisch & Pock, ; Gould et al., ; Alesiani, ) and meta-learning where the outer learning process is optimized (Finn, ; Andrychowicz et al., ; Hochreiter et al., ; Hospedales et al., ; Rajeswaran et al., ; Sambharya et al., ).
Alternative approaches avoid implicit differentiation by using approximation techniques. For instance, they apply automatic differentiation directly to iterative algorithms through loop unrolling (Belanger & McCallum, ; Belanger et al., ; Metz et al., ; Scieur et al., ) or, in the case of fixed-point mappings, by differentiating a single iteration or employing a Jacobian-free method (Geng et al., ; Fung et al., ; Bolte et al., ).

Sensitivity Analysis and Parametric Programming.
There is extensive mathematical theory on the local behavior of optimization problems under perturbations (Rockafellar, ; Rockafellar & Wets, ), particularly in assessing the sensitivity and stability of solutions (Fiacco, ; ; Fiacco & McCormick, ; Lee et al., ; Bonnans & Shapiro, ). For this, the implicit function theorem is a central analytical tool, but unlike fixed point mappings, an intermediate step is required before it can be applied. Particularly, one must first pass from the optimization problem itself to its optimality conditions, often requiring a number of assumptions in order for them to be completely equivalent. The theoretical results of sensitivity theory underpin applications in multi-parametric programming (Pistikopoulos et al., ), like model predictive control, where the problem is solved for various input parameters, leading to intense computations. To address this, Bemporad et al. () observed that QP systems have a closed-form solution if the active or binding set is known beforehand, allowing it to be pre-computed offline. This approach requires partitioning the parameter space into regions of fixed active set (Spjøtvold et al., ), inside which the active set is stable to perturbations. Methods based on this idea continue to be developed for solving parametric QPs (Ferreau et al., ; Narciso et al., ; Arnström & Axehill, ).

Differentiable Programming.
Our work follows OptNet (Amos & Kolter, ; Amos, ) which differentiates QPs through their optimality conditions, and focuses on small dense problems for GPU batching. They solve the full Jacobian system efficiently by reusing the factorization employed in their custom interior-point method. However, as noted in (Bambade et al., ), this comes at the cost of ill-conditioning due to symmetrization. More recent differentiable QP methods include Alt-Diff and SCQPTH (Sun et al., ; Butler & Kwon, ; Butler, ) which use first-order ADMM and approximately differentiate the fixed point map, and QPLayer (Bambade et al., ) focusing on accommodating infeasibility via extended conservative Jacobians.
Several works have highlighted the importance of the active constraint set in differentiating constrained optimization problems (Amos et al., ; Gould et al., ; Paulus et al., ), as well as in the context of quadratic programming (Amos et al., ; Bambade et al., ; Pan et al., ; Niculae et al., ). A common observation is that the algebraic system obtained through implicit differentiation can be simplified by removing rows corresponding to inactive constraints. Amos et al. (); Pan et al. () have additionally observed that backpropagation can be cast as an equality-constrained QP parameterized by incoming gradients. However, existing approaches do not directly utilize the formation of a significantly dimension-reduced symmetric linear system to efficiently differentiate arbitrary black-box QP solvers, thus missing the opportunity to effectively decouple optimization and differentiation.
Other classes of optimization problems, such as convex cone programs (Agrawal et al., ) and mixed-integer programs (Paulus et al., ), have also been differentiated. Frameworks proposed in (Agrawal et al., ; Blondel et al., ; Pogančić et al., ; Besançon et al., ) provide a differentiable interface to broader classes of optimization problems. However, these frameworks often offer limited support for specialized efficient solvers, may require access to a dual solution, or invoke costly symbolic compilation, making them most suitable for prototyping.

## 3 Approach

We now detail the theoretical underpinning of our method that can transform any black-box QP solver into a differentiable layer. We begin by formulating the problem concretely, move on to establishing basic theory of differentiation of QP’s (via sensitivity analysis and KKT conditions), and finally connect those with our main novel theoretical observation, leading to a straightforward algorithm. We note that various subparts of the theoretical background discussed in the following have been used in recent years to develop differentiable QP layers, see Section 2. However, no single work has fully leveraged this theory to completely decouple optimization and differentiation in a manner that supports arbitrary state-of-the-art QP solvers while also ensuring efficient, robust differentiation.

### 3.1 Problem Setup: Differentiating Quadratic Programs

We consider a quadratic program which is feasible and strictly convex (i.e., $P\succ 0$) in standard form,

$$ $\displaystyle z^{*}(\theta)=\operatorname*{arg\,min}_{z}$ $\displaystyle\frac{1}{2}z^{T}P(\theta)z+q(\theta)^{T}z$ (1) subject to $\displaystyle A(\theta)z=b(\theta)$ $\displaystyle C(\theta)z\leq d(\theta),$ $$

where $P\in\mathbb{R}^{n\times n},q\in\mathbb{R}^{n},A\in\mathbb{R}^{p\times n},b\in
\mathbb{R}^{p},C\in\mathbb{R}^{m\times n}$ and $d\in\mathbb{R}^{m}$ are smoothly parameterized by some $\theta\in\mathbb{R}^{s}$. This $\theta$ can either be the output of a previous layer in a neural network, or otherwise learnable parameters. To simplify notation, in the following we omit $\theta$.

To motivate our work, consider the case in which a QP of the form Equation (1) is incorporated as the $\ell$-th layer of a neural network, i.e., the QP layer receives an input vector $x_{\ell}$ and outputs a vector $x_{\ell+1}$ satisfying the relation
$x_{\ell+1}=z^{*}\left(x_{\ell}\right)$.
In other words, the QP layer’s input, $x_{\ell}$, serves as the parameters $\theta$ that control the objective and constraints of the QP, and the optimal point $z^{*}(x_{\ell})$ is the layer’s output. Training a neural network requires backpropagating gradients through it, which involves computing the Jacobian of the layer’s output with respect to its input, $\frac{\partial x_{\ell+1}}{\partial x_{\ell}}$. In the case of a QP layer, these gradients are exactly $\frac{\partial z^{*}}{\partial{\theta}}$, i.e., the derivative of the optimal point $z^{*}$ with respect to the parameters $\theta$. The same derivative is also essential when using descent methods to solve certain bi-level optimization problems (Colson et al., ).

In this work, we focus on computing $\partial_{\theta}z^{*}(\theta)=\frac{\partial}{\partial\theta}z^{*}(\theta)$, the derivative of the optimal point of the QP Equation (1) with respect to the parameters $\theta$. Intuitively, this derivative quantifies the change in the optimal point of the QP in response to a perturbation of its parameters $\theta$. Our goal is to efficiently compute $\partial_{\theta}z^{*}(\theta)$ independently of the method used to approximate the optimal point $z^{*}(\theta)$.

### 3.2 Theoretical Differentiation of QPs via KKT Conditions and Sensitivity

Our goal is to devise a method for differentiating QPs based solely on the solution provided by a black-box numerical solver. First, we need to establish the theoretical foundations necessary for the desired derivatives. These derivatives, as is common in optimization, are obtained through sensitivity analysis applied to the KKT conditions. In this section, we elaborate on these concepts, synthesizing key theoretical insights from optimization, sensitivity analysis, parametric programming, and differentiable programming techniques, distilling them in the context of QPs to lay the groundwork for our results and the development of dQP.

##### Optimality Conditions.

The first-order Karush–Kuhn–Tucker (KKT) conditions (Karush, ; Kuhn & Tucker, ; Boyd & Vandenberghe, ; Wright, ) provide a useful algebraic characterization of the optimal points of constrained optimization problems. In essence, they are an extension of the method of Lagrange multipliers for problems that include inequalities. For the QP Equation (1), the KKT conditions take the form,

$$ $\displaystyle Pz^{*}+q+A^{T}\lambda^{*}+C^{T}\mu^{*}$ $\displaystyle=0$ (2) $\displaystyle Az^{*}-b$ $\displaystyle=0$ $\displaystyle Cz^{*}-d$ $\displaystyle\leq 0$ $\displaystyle\mu^{*}$ $\displaystyle\geq 0$ $\displaystyle D(\mu^{*})(Cz^{*}-d)$ $\displaystyle=0,$ $$

where $D(\mu^{*})=\textrm{diag}(\mu^{*})$, and the additional added variables $\lambda^{*}\in\mathbb{R}^{p}$ and $\mu^{*}\in\mathbb{R}^{m}$ are called the optimal dual variables of the linear equalities and inequalities, respectively. With these dual variables, one considers the extended primal-dual solution $\zeta^{*}(\theta)=(z^{*}(\theta),\lambda^{*}(\theta),\mu^{*}(\theta))$ of the QP Equation (1). Crucially, under strict convexity and feasibility, the QP Equation (1) has a unique solution $\zeta^{*}(\theta)$, and the KKT conditions Equation (2) are necessary and sufficient for its optimality.

##### Active Set and Complementary Slackness.

A main point of interest in this work lies in the last equation of Equation (2), which is the nonlinear *complementary slackness* condition. Intuitively, it encodes the two situations in which each original inequality constraint from Equation (1), $(Cz^{*}-d)_{j}\leq 0$, may be. Either (1) the constraint is *active*, i.e., it is satisfied as an equality $(Cz^{*}-d)_{j}=0$, in which case $\mu^{*}_{j}\geq 0$; or (2) the constraint is *inactive*, i.e., it is satisfied with a strict inequality, in which case $\mu^{*}_{j}=0$. Importantly, an inactive constraint implies that the same optimal solution $z^{*}$ would be obtained even if that specific constraint were removed from the QP. We denote by $J(\theta)=\{j:\left(C(\theta)z^{*}(\theta)-d(\theta)\right)_{j}=0\}$ the set of active constraints.

##### Derivatives via Sensitivity Analysis.

To define derivatives of QPs, we turn to the Basic Sensitivity Theorem (Theorem 2.1 in Fiacco ()), which provides the foundation for differentiating the KKT conditions with respect to $\theta$. To differentiate at $\theta$, the theorem requires the additional condition of strict complementary slackness; this prohibits the degenerate case where both $(Cz^{*}-d)_{j}=0$ and $\mu^{*}_{j}=0$, ensuring that a small perturbation of the solution does not alter the active set. Under strict complementary slackness, it establishes that *in a neighborhood of $\theta$*, the primal-dual point $\zeta^{*}(\theta)=(z^{*}(\theta),\lambda^{*}(\theta),\mu^{*}(\theta))$ is a differentiable function of $\theta$, optimal for the QP Equation (1), uniquely satisfies the KKT conditions Equation (2), and maintains strict complementary slackness. Crucially, the active set $J(\theta)$ remains fixed within this neighborhood.

Since the active set is stable, the equality conditions in Equation (2) suffice to provide a local characterization of $\zeta^{*}(\theta)$. Implicit differentiation of these yields the Jacobians of the solution $\partial_{\theta}\zeta^{*}$ in terms of the following linear system,

$$ $\displaystyle\begin{bmatrix}P&A^{T}&C^{T}\\ A&0&0\\ D(\mu^{*})C&0&D(Cz^{*}-d)\\ \end{bmatrix}\begin{bmatrix}\partial_{\theta}z^{*}\\ \partial_{\theta}\lambda^{*}\\ \partial_{\theta}\mu^{*}\end{bmatrix}=-\begin{bmatrix}\partial_{\theta}Pz^{*}+ \partial_{\theta}q+\partial_{\theta}A^{T}\lambda^{*}+\partial_{\theta}C^{T}\mu ^{*}\\ \partial_{\theta}Az^{*}-\partial_{\theta}b\\ D(\mu^{*})(\partial_{\theta}Cz^{*}-\partial_{\theta}d)\end{bmatrix}.$ (3) $$

Under the conditions for the Basic Sensitivity Theorem, the linear system Equation (3) is invertible. It degenerates exactly in the presence of weakly active constraints $\mu_{j}^{*}=(Cz^{*}-d)_{j}=0$, for which the QP is non-differentiable (see, e.g., Amos & Kolter ()). For any inactive constraint $j\notin J$, the dual variable $\mu_{j}^{*}$ vanishes, and thus the corresponding rows and columns of Equation (3) can be removed, simplifying it into the reduced form

$$ $\displaystyle\begin{bmatrix}P&A^{T}&C_{J}^{T}\\ A&0&0\\ C_{J}&0&0\\ \end{bmatrix}\begin{bmatrix}\partial_{\theta}z^{*}\\ \partial_{\theta}\lambda^{*}\\ \partial_{\theta}\mu_{J}^{*}\end{bmatrix}=-\begin{bmatrix}\partial_{\theta}Pz^ {*}+\partial_{\theta}q+\partial_{\theta}A^{T}\lambda^{*}+\partial_{\theta}C_{J }^{T}\mu_{J}^{*}\\ \partial_{\theta}Az^{*}-\partial_{\theta}b\\ \partial_{\theta}C_{J}z^{*}-\partial_{\theta}d_{J}\end{bmatrix},$ (4) $$

where $\mu^{*}_{J}$, $C_{J}$ and $d_{J}$ denote restriction to rows corresponding to active inequality constraints $j\in J$.

### 3.3 Extracting derivatives from a QP solver’s solution

Through the above theory, we can obtain our main theoretical results and introduce dQP, a straightforward algorithm for efficient and robust differentiation of any black-box QP solver.

Our approach stems from two straightforward yet powerful insights: (1) given the primal solution of a QP, it is easy to identify the active set of a QP; (2) once the active set is known, both the primal-dual optimal point *and* its derivatives can be explicitly derived in closed-form. Furthermore, the computation of these quantities can then be achieved efficiently using a single matrix factorization of a reduced-dimension symmetric matrix.

These observations in turn lead to a simple algorithm that is easy to implement: first, solve the optimization problem using *any* QP solver; then, use the solution to identify the active set and solve a linear system to compute the derivatives. Consequently, we can define a “backward pass” for any layer that uses a QP solver, allowing for the seamless integration of any solver best suited to the problem, thus leveraging years of research and development invested in state-of-the-art QP solvers.

Figure: Figure 3: Schematic active set differentiation. Left: a QP is shown by its quadratic level sets and polyhedral feasible set; the solution lies on a facet of the boundary; perturbations of the constraints lead to perturbations in the solution. Right: the perturbation of the solution remains the same when inactive constraints are eliminated.
Refer to caption: x3.png

##### Explicit Active Set Differentiation.

Consider a QP and its optimal point $\zeta^{*}(\theta)$, along with the set $J(\theta)$ of active constraints (see Section 3.2). We define the reduced equality-constrained quadratic program, obtained by removing inactive inequalities and converting active inequality constraints into equality constraints,

$$ $\displaystyle z^{*}(\theta)=\operatorname*{arg\,min}_{z}$ $\displaystyle\frac{1}{2}z^{T}P(\theta)z+q(\theta)^{T}z$ (5) subject to $\displaystyle\begin{bmatrix}A(\theta)\\ C(\theta)_{J(\theta)}\end{bmatrix}z=\begin{bmatrix}b(\theta)\\ d(\theta)_{J(\theta)}\end{bmatrix}.$ $$

Under the assumptions of Section 3.2, this simpler QP is, in fact, *locally* equivalent to the QP Equation (1). Moreover, it provides an explicit expression for both the primal-dual optimal point and its derivatives:

###### Theorem 1 .

The QP Equation (5) is locally equivalent to the reduced equality-constrained QP Equation (1) and its solution $\zeta^{*}(\theta)=(z^{*}(\theta),\lambda^{*}(\theta),\mu^{*}(\theta))$ admits the explicit form

$$ $\begin{bmatrix}z^{*}\\ \lambda^{*}\\ \mu_{J}^{*}\end{bmatrix}=\begin{bmatrix}P&A^{T}&C_{J}^{T}\\ A&0&0\\ C_{J}&0&0\\ \end{bmatrix}^{-1}\begin{bmatrix}-q\\ b\\ d_{J}\end{bmatrix}.$ (6) $$

Furthermore, the optimal point can be explicitly differentiated to obtain

$$ $\begin{bmatrix}\partial_{\theta}z^{*}\\ \partial_{\theta}\lambda^{*}\\ \partial_{\theta}\mu_{J}^{*}\end{bmatrix}=-\begin{bmatrix}P&A^{T}&C_{J}^{T}\\ A&0&0\\ C_{J}&0&0\\ \end{bmatrix}^{-1}\left(\begin{bmatrix}\partial_{\theta}P&\partial_{\theta}A^{ T}&\partial_{\theta}C_{J}^{T}\\ \partial_{\theta}A&0&0\\ \partial_{\theta}C_{J}&0&0\\ \end{bmatrix}\begin{bmatrix}z^{*}\\ \lambda^{*}\\ \mu_{J}^{*}\end{bmatrix}-\begin{bmatrix}-\partial_{\theta}q\\ \partial_{\theta}b\\ \partial_{\theta}d_{J}\end{bmatrix}\right).$ (7) $$

A proof of this Theorem, based on the Basic Sensitivity Theorem (Fiacco, ), is provided in Appendix A, along with a calculation of the derivatives using differential matrix calculus (Petersen & Pedersen, ; Magnus & Neudecker, ). We note that this result is closely related to analyses studied in multi-parametric programming (Bemporad et al., ; Pistikopoulos et al., ; Spjøtvold et al., ; Arnström & Axehill, ; Narciso et al., ).

Notably, in the case of quadratic programming, the Basic Sensitivity Theorem allows one to bypass the need for implicit differentiation techniques. Theorem 7 underscores a critical computational insight: once a black-box solver provides the primal solution to the QP, the active set can be determined, and additionally the derivatives can be computed via Equation (7). Furthermore, if the solver provides only the primal solution and not the primal-dual pair, the dual can be completed through Equation (6), with negligible extra cost due to prefactorization. All these insights lead up to the key algorithm of dQP, summarized in Algorithm 1.

Figure: Algorithm 1 – dQP: Differentiation through Black-box Quadratic Programming Solvers

###### Theorem 1 .

##### Numerical Computation.

Our approach leads to a compact and efficient computation of gradients.
Indeed, the linear system Equation (7) that we factorize to compute the derivatives and dual solution is symmetric and reduced in size. In contrast, implicit differentiation of the full KKT conditions Equation (1) leads to a significantly larger, asymmetric system Equation (3). Furthermore, our approach not only simplifies computation but also allows for the use of fast, specialized linear solvers that exploit the problem’s structure (e.g., (Schenk & Gärtner, ; Amestoy et al., )).

[Uncaptioned image]: x4.png

Empirically, we observe that the reduced linear system Equation (7) is often significantly better conditioned than its full counterpart. The inset figure illustrates this with an example of a QP governed by two parameters $\theta=(\theta_{1},\theta_{2})$ from Spjøtvold et al. (), computed using DAQP (Arnström et al., ). The figure visualizes (a) regions in which the active set is constant, and (b) the conditioning of the full and reduced linear systems along a cross-section in parameter space (right), demonstrating that eliminating inactive constraints improves conditioning.

This figure also highlights the challenge of calculating derivatives near singularities, where the active set changes and some inequalities turn into weakly active, leading to ill-defined derivatives. Near such singularities, implicit differentiation suffers from severe ill-conditioning. This affects our approach as well, manifesting in the challenge of determining the active set at an approximate solution. Various methods have been proposed to address this issue, such as specialized algorithms for active set identification (Cartis & Yan, ; Oberlin & Wright, ; Burke & Moré, ). Our implementation includes an optional heuristic for active set refinement to address this instability, described in Appendix B. However, we found that simple hard thresholding of the primal residual $r_{j}=(Cz^{*}-d)_{j}\geq-\epsilon_{J}$ was sufficiently robust in all of our experiments, as shown in Section 4.

##### Implementation.

Our open-source implementation will be made publicly available. We implement dQP, Algorithm 1, as a fully differentiable module in PyTorch (Paszke et al., ), providing a simple-to-use interface for easily integrating differentiable QPs into machine learning algorithms or bi-level programming. As a PyTorch module, it is necessary to render Equation (7) as a backpropogation step, which we describe in Appendix C. To ensure modularity, we offer complete flexibility in selecting a QP solver for the forward pass by interfacing with the open-source qpsolvers library (Caron et al., ). Their library provides a minimal-overhead interface supporting over 15 free and commercial QP solvers, and easily supports the integration of additional solvers. We provide flexibility in choosing the linear solver that will be used for differentiation: our code supports several popular direct linear solvers, as well as solvers designed for symmetric indefinite KKT systems, such as QDLDL (Stellato et al., ; Davis, ). Additionally, our implementation offers full end-to-end support for both dense and sparse problems with appropriate QP and linear solvers. For users who wish to determine the “best” QP solver for their problem, dQP includes a simple profiling tool (see Appendix D). More details are given in Appendix E including constraint normalization, handling non-differentiable points, and options like warm-starting for bi-level optimization.

## 4 Experimental Results

We have extensively tested dQP to verify its robustness, evaluate its performance against other methods for differentiable quadratic programming, and demonstrate its applicability and advantages in large-scale structured problems. Our evaluation includes a large benchmark consisting of over 2,000 dense and sparse challenging QPs taken from public datasets as well as randomly generated problems, designed to test dQPs robustness and performance. We also present two prototype applications, demonstrating the applicability of dQP in a learning experiment and in bi-level optimization. See Appendix G for the full details on each experiment’s configuration.

Figure: Figure 4: Accuracy versus total forward/backward solve for the (a) MPC and (b) Maros-Meszaros datasets. Each point represents a solved problem; point size illustrates dimension; problems solved solely by dQP circled. Legend shows percentages of success rates, and counts of forward solvers used by dQP for each problem.
Refer to caption: x5.png

##### Modularity and Performance.

We tested dQP on nearly 200 QPs from the QP benchmark (Caron et al., ), including 65 small Model Predictive Control (MPC) problems and 126 challenging, moderately-sized sparse problems from the standard Maros-Meszaros (MM) dataset (Maros & Mészáros, ), designed as a stress test for QP solvers. We compared dQP with other differentiable QP methods: OptNet (Amos & Kolter, ), QPLayer (Bambade et al., ), and SCQPTH (Butler, ), each intergrated with its specialized QP solver. For its forward pass, dQP was paired with the leading QP solver for each problem as reported by QP Benchmark (Caron et al., ).
The total runtime (forward and backward passes), accuracy (duality gap), and dimension (illustrated by point size) are reported in the scatter plots in Figure 4 for each problem and each differentiable QP solver. For small, dense MPC problems, dQP is typically faster than or comparable to QPLayer while achieving much higher accuracy. MM problems, being significantly more challenging, often cause competing methods to fail. OptNet and SCQPTH solved less than 50% of the problems, while dQP successfully solved all MM problems and was the *only* differentiable solver to succeed in 38 of them (circled in the figure). Moreover, dQP was the fastest and most accurate in 81% and 83% of all problems, respectively. It particularly excelled in larger problems (dimension over 1000), where it was the fastest and most accurate in 98% and 95% of cases, respectively. Further technical details, along with additional experiments on 450 random dense QPs and 625 sparse QPs with dimensions ranging from $10$ to $10^{4}$, are provided in Appendix G.1.1, showing similar results.

##### Scalability.

We evaluated dQP on large-scale sparse problems, a regime where state-of-the-art QP solvers hold a significant advantage over less optimized solvers. We tested dQP and other available differentiable QP solvers on two prototypical projection layers expressed as constrained QPs:

$$ $P_{1}(x)=\operatorname*{arg\,min}_{z}\left\|x-z\right\|_{2}^{2}\ {\textrm{ subject to}}\ 0\leq z\leq 1,{\textstyle\sum}z_{i}=1\textrm{, and}$ $$

$$ $P_{2}(x_{1},\ldots,x_{n})=\operatorname*{arg\,min}_{z_{1},\ldots,z_{n}}{ \textstyle\sum}\left\|x_{j}-z_{j}\right\|_{2}^{2}\ {\textrm{subject to}}\ \left\|z_{j}-z_{j+1}\right\|_{\infty}\leq 1.$ $$

Results for $P_{1}$ are shown in Figure 1, demonstrating dQP’s scalability compared to OptNet and QPLayer. Other methods fail on all but small problems (see Appendix G.1.2). In dimensions greater than 2000, dQP outperforms competing methods by 2-3 orders of magnitude in both speed and accuracy. Competing methods are limited to dense calculations and fail in dimensions beyond $10^{4}$. It’s worth noting that $P_{1}$ is the projection onto the probability simplex, also known as SparseMAX, for which more efficient, non-QP-based methods exist (Martins & Astudillo, ). Results for $P_{2}$, representing projection onto ”chains” with bounded links, exhibit similar scalability and are detailed in Appendix G.1.3.

[Uncaptioned image]: x6.png

##### Learning Sudoku.

We evaluated dQP in a popular learning setting, first introduced in OptNet (Amos & Kolter, ). In this experiment, linear constraints model the rules of the Sudoku game, which are then learned from examples of solved Sudoku boards via differentiable QPs. We reproduced the experiment from (Amos & Kolter, ) by integrating dQP into their code. As shown in the inset, the training and testing losses achieved by OptNet and dQP paired with PIQP (Schwan et al., ) are comparable. This experiment similar to the MPC results shown in Figure 4(a), validates dQP’s ability to a perform on par with leading differentiable QP packages. However, we note that, compared to tightly integrated forward-backward implementations, dQP has some disadvantages, e.g., QPLayer supports differentiation of infeasible QPs, and OptNet natively supports GPU batching, which is not available for black-box state-of-the-art sparse QP solvers, and thus cannot be easily integrated with dQP.

Figure: Figure 5: Solver speed under mesh refinement for mapping into a non-convex perturbed square, visualized at different resolutions.
Refer to caption: x7.png

##### Bi-Level Geometry Optimization.

We further test dQP in a non-learning, optimization-based setting inspired by the geometric problem of intersection-free straight-edge planar graph drawing: embed a planar graph representing a triangular mesh into a non-convex domain, such edges are drawn as straight non-overlapping lines. Kovalsky et al. () formulate a linear inequality constraint-satisfaction problem for which they show exists an (unknown) Laplacian $M$ that defines a quadratic energy and thus a QP which, when solved, yields exactly such a straight-edge drawing. However, their conditions are nonconstructive and have remained theoretical. We cast this problem as the following bi-level optimization problem:

$$ $\displaystyle M^{*}=\operatorname*{arg\,min}_{M}\$ $\displaystyle\left\|\mu^{*}(M)\right\|_{2}^{2}$ subject to $\displaystyle(v^{*}(M),\lambda^{*}(M),\mu^{*}(M))=\operatorname*{arg\,min}_{v} \left\{\mathrm{tr}\left(v^{T}Mv\right)\ {\textrm{subject to}}\ Bv=u,\ CMv \succeq 0\right\}$ $$

where $v\in\mathbb{R}^{n\times 2}$ represents the $n$ coordinates of the mesh vertices, $M\in\mathbb{R}^{n\times n}$ is a parameterized Laplacian, and $B$,$u$ and $C$ encode the boundary conditions of Kovalsky et al. (). The results of Kovalsky et al. () then imply that $v^{*}(\mathcal{L^{*}})$ represents a straight line drawing if the dual variable $\mu^{*}(M)$, corresponding to linear inequalities of the nested optimization, vanishes.

[Uncaptioned image]: x8.png

In our experiments, we solve this bi-level problem using dQP paired with PIQP.
The inset shows an example of this experiment: (a) the triangulated unit square is the chosen graph; (b) an invalid embedding produced by choosing an arbitrary Laplacian; (c) a valid embedding which minimizes the above bi-level problem; (d) with additional regularization on the shape of the triangles. Figure 5 shows a refinement experiment showing that dQP scales favorably as mesh size increases compared to OptNet, QPLayer and SCQPTH; in particular, only dQP scales up to problems with over $10^{4}$ vertices. We only report forward (QP) time for OptNet, QPLayer and SCQPTH because OptNet and SCQPTH do not output, nor differentiate the duals, and while QPLayer does, it suffers poor scaling from dense operations as the others. Lastly, Figure 2 presents the large-scale bijective embedding of an ant mesh.

## 5 Conclusion

dQP is shown to provide a differentiable interface to any QP solver, and yield an extremely efficient QP-based layer which can be used in, e.g., neural architectures. We believe this work is the first step in providing similar differentiable layers for other popular optimization problems (e.g., semidefinite programming), which we plan to tackle next. Additionally, we note that our current method does not enable neither full parallelization nor GPU support, and we mark these as important challenges to tackle.

## Appendix A Proof of Theorem 1

In this section we provide a proof of Theorem 1, which we restate below:
See 1

###### Proof.

We begin by establishing that the QP Equation (1) and the equality-constrained reduced QP Equation (5) are equivalent. For any $\theta$ satisfying the assumptions of the theorem, the QP Equation (1) has a unique solution characterized by the KKT system

$$ $\displaystyle P(\theta)z^{*}(\theta)+q(\theta)+A(\theta)^{T}\lambda^{*}(\theta )+C(\theta)^{T}\mu^{*}(\theta)$ $\displaystyle=0$ (8) $\displaystyle A(\theta)z^{*}(\theta)-b(\theta)$ $\displaystyle=0$ $\displaystyle C(\theta)z^{*}(\theta)-d(\theta)$ $\displaystyle\leq 0$ $\displaystyle\mu^{*}(\theta)$ $\displaystyle\geq 0$ $\displaystyle D(\mu^{*}(\theta))(C(\theta)z^{*}(\theta)-d(\theta))$ $\displaystyle=0.$ $$

Complementarity implies that active constraints $j\in J(\theta)$ have $\mu^{*}(\theta)_{j}>0$ and therefore must be satisfied with an equality $\left(C(\theta)z^{*}(\theta)-d(\theta)\right)_{j}=0$, while inactive constraints $j\notin J(\theta)$ have $\mu^{*}(\theta)_{j}=0$ and thus can be eliminated, without altering the solution. Therefore, the unique solution $\zeta^{*}(\theta)=(z^{*}(\theta),\lambda^{*}(\theta),\mu^{*}(\theta))$ of Equation (8) is also the unique solution of the reduced system

$$ $\displaystyle P(\theta)z^{*}(\theta)+q(\theta)+A(\theta)^{T}\lambda^{*}(\theta )+C(\theta)_{J(\theta)}^{T}\mu^{*}(\theta)_{J(\theta)}$ $\displaystyle=0$ (9) $\displaystyle A(\theta)z^{*}(\theta)-b(\theta)$ $\displaystyle=0$ $\displaystyle C(\theta)_{J(\theta)}z^{*}(\theta)-d(\theta)_{J(\theta)}$ $\displaystyle=0,$ $$

which are exactly the KKT conditions of the equality-constrained reduced QP Equation (5). Uniqueness of solution then implies that Equation (1) and Equation (5) are pointwise equivalent at $\theta$. Moreover, since $P,q,A,b,C,d$ are smoothly parameterized by $\theta$, the Basic Sensitivity Theorem (Fiacco, ) asserts that the primal-dual solution $\zeta^{*}(\theta)$ for Equation (1) is a differentiable function of $\theta$ in a neighborhood of $\theta$, defined implicitly through the KKT’s equality conditions. Furthermore, the active set $J(\theta)$ is fixed in this neighborhood, therefore Equation (1) and Equation (5) are locally equivalent.

Equation (9) implies that the reduced primal-dual solution $\zeta_{J}^{*}(\theta)=(z^{*}(\theta),\lambda^{*}(\theta),\mu^{*}_{J}(\theta))$ satisfies $K_{J}(\theta)\zeta_{J}^{*}(\theta)=v_{J}(\theta)$, where

$$ $K_{J}(\theta)=\begin{bmatrix}P(\theta)&A(\theta)^{T}&C(\theta)^{T}_{J(\theta)} \\ A(\theta)&0&0\\ C(\theta)_{J(\theta)}&0&0\\ \end{bmatrix},\quad v_{J}(\theta)=\begin{bmatrix}-q(\theta)\\ b(\theta)\\ d(\theta)_{J(\theta)}\end{bmatrix}.$ (10) $$

Under the assumptions of the theorem, the reduced KKT matrix $K_{J}(\theta)$ is invertible and

$$ $\zeta_{J}^{*}=K_{J}^{-1}v_{J},$ (11) $$

yielding Equation (6). Moreover, since $J(\theta)$ is locally constant, the Basic Sensitivity Theorem establishes that $\zeta_{J}^{*}(\theta)$ is differentiable. Using the formal derivative of the matrix inverse (Magnus & Neudecker, ; Petersen & Pedersen, ) we explicitly differentiate Equation (11) to obtain

$$ $\partial_{\theta}\zeta_{J}^{*}=(-K_{J}^{-1}(\partial_{\theta}K)K_{J}^{-1})v_{J }+K_{J}^{-1}(\partial_{\theta}v_{J})=-K_{J}^{-1}(\partial_{\theta}K_{J})\zeta_ {J}^{*}+K_{J}^{-1}(\partial_{\theta}v_{J}),$ (12) $$

yielding Equation (7).

∎

###### Proof.

## Appendix B Active Set Refinement

Figure: Figure 6: The set-up in figure 3.3 with looser solver tolerance $\epsilon_{\textrm{abs}}=10^{-4}$, active tolerance $\epsilon_{J}=10^{-7}$, and solver PIQP. (a) The computed active set is degraded due to the inaccurate solution. (b) Our heuristic active set refinement algorithm recovers the ground truth active sets.
Refer to caption: x9.png

Inaccuracy in a solution may lead to instability in the active set near weakly active constraints, degrading the gradient quality. To show this, we repeat the experiment in Figure 3.3 which has a simple polyhedral active set parameter space. One setup where instability appears is illustrated in Figure 6 where we use absolute solver tolerance $\epsilon_{\textrm{abs}}=10^{-4}$ and active tolerance $\epsilon_{J}=10^{-7}$. Qualitatively, the active set at each solution is severely degraded, even for points away from the boundaries where the set changes. We provide a *optional* heuristic algorithm to address this, which recovers the desired set in this problem. First, we order the constraints by increasing residual and select an initial active set from the tolerance $\epsilon_{J}$. Then, we progressively add constraints by checking if the residual of the system 6 for $\zeta_{J}^{*}$ decreases, and greedily accepting until adding constraints no longer improves the residual. At each step, we keep the primal solution from the forward fixed, and solve for the new active dual variables. While this algorithm works well on simple examples, more sophisticated and efficient techniques may be desired for harder problems. We did not use this refinement algorithm in any of our experiments.

## Appendix C Backpropagation

Like other differentiable QP layers implemented within automatic differentiation frameworks such as PyTorch (Paszke et al., ), we do not directly compute the derivative $\partial_{\theta}\zeta^{*}$. Specifically, dQP directly receives the QP parameters $P,q,A,b,C,d$ and not $\theta$, and so in backpropogation we are not concerned with $\theta$. This is rather accounted for in the next step outside dQP, usually by automatic differentiation. Instead, backpropogation requires that we compute a so-called Jacobian-vector product which are products of the Jacobians with an “incoming” gradient of a quantity or loss $\ell$ that depends on $\zeta^{*}$. This requires less computation and does not require the formation of a 3-tensor. Since $\zeta_{J}^{*}=K_{J}^{-1}v_{J}$ is a formal matrix-vector multiplication, the Jacobian-vector product is well-known,

$$ $\nabla_{v_{J}}\ell=(K_{J}^{-1})^{T}\,\,\nabla_{\zeta_{J}^{*}}\ell,$ (13) $$

and

$$ $\nabla_{K_{J}}\ell=-\nabla_{v_{J}}\ell\,\,{\zeta_{J}^{*}}^{T},$ (14) $$

with respect to $K_{J},v_{J}$, respectively. Although backpropogation introduces a transposition, the re-use of a factorization from solving for the active duals is unaffected. This follows from the symmetry of the reduced KKT matrix which simplifies Equation (13) into $\nabla_{v_{J}}\ell=K_{J}^{-1}\,\,\nabla_{\zeta_{J}^{*}}\ell$. Next, we extract the gradients with respect to the parameters by the chain rule. This amounts to tracking their position in the blocks and accounting for symmetry constraints. It is helpful to write $(d_{z},d_{\lambda},d_{\mu_{J}})=-\nabla_{v_{J}}\ell$ so that we express

$$ $\displaystyle\nabla_{P}\ell$ $\displaystyle=\frac{1}{2}\left(d_{z}\,z^{*T}+z^{*}\,d_{z}^{T}\right)$ $\displaystyle\nabla_{q}\ell=d_{z}$ $\displaystyle\nabla_{A}\ell$ $\displaystyle=d_{\lambda}z^{*T}+\lambda^{*}d_{z}^{T}$ $\displaystyle\nabla_{b}\ell=-d_{\lambda}$ (15) $\displaystyle\nabla_{C_{J}}\ell$ $\displaystyle=d_{\mu_{J}}z^{*T}+\mu_{J}^{*}d_{z}^{T}$ $\displaystyle\nabla_{d_{J}}\ell=-d_{\mu_{J}},$ $$

similar to OptNet (Amos & Kolter, ). We note that the gradient with respect to $P$ is constrained to lie within the subspace of symmetric matrices. Similarly, if the matrices $P,A,C$ are sparse, then we project the gradient to lie within the non-zero entries, which can be implemented efficiently in Equations 15. Although the above argument is for a scalar loss $\ell$, the same approach is naturally adapted if $\zeta^{*}$ is mapped to a vector in the immediate next layer.

## Appendix D Choosing a Solver

Figure: Figure 7: Evaluating the best QP solver for the cross geometry problem using our diagnostic tool. The solution tolerance regimes are varied between $\epsilon_{\textrm{abs}}=10^{-8},10^{-5},10^{-2}$.
Refer to caption: x10.png

Since our work enables users to choose any QP solver as the front-end for their differentiable QP applications, we include a simple diagnostic tool for quantitatively measuring solver performance. We present an example result in Figure 7 for the cross geometry experiment in section 4, finding PIQP, OSQP, and QPALM to be the most efficient. For this reason, we choose PIQP in the geometry experiments. We also include tools for checking the solution and gradient accuracy.

## Appendix E Implementation Details

Tolerances In addition to the active set tolerance $\epsilon_{J}$, QP solvers often support additional user-provided tolerances. These include the primal residual which measures violations of feasibility, the dual residual which measures violations of stationary, and for some solvers also the duality gap, which provides a direct handle on solution accuracy. We inherit the structure of qpsolvers for setting custom tolerances on different QP solvers, though we set a heuristic default which is sufficient for many of the experiments in this work.

Convexity and Feasibility Two key assumptions of our method are strict convexity and feasibility. However, these are often violated in practice. We include optional checks that $P$ is symmetric positive definite. On the other hand, we do not perform any special handling for infeasibility – a limitation of our method compared to, for example, QPLayer (Bambade et al., ).

Non-differentiable Points For non-differentiable problems, we solve for the derivatives in the least-squares sense, plugging the system into qpsolvers which can handle least-squares, or a standard least-squares solver. We attempt to anticipate weakly active constraints which cause non-differentiability by measuring the norms of the primal residual and the dual. The reduced KKT is also non-invertible if the active dual solution is not unique. To detect this, we check a necessary condition: the total number of active constraints plus the number of equalities must be less than the dimension. Otherwise, if these necessary checks are passed, we attempt the standard linear solve and pass to least-squares if it fails.

Normalization Some problems have large variations in scale between different rows within the constraints. This influences the primal residual and thus the active set, which is determined by comparing with an absolute threshold tolerance. To address this issue for these problems, we include an *optional* differentiable normalization step on the constraints before Algorithm 1 is carried out. Under this choice, the resulting relative primal residual becomes the scale-invariant distance to the constraint.

Equality Constraints While we include equality constraints in our general formulation, they are not required.

Warm-Start Since qpsolvers supports warm-starting, we inherit it as an option and store data in the PyTorch module from previous outer iterations, which can be used as initialization. This is useful for bi-level optimization problems where the input $\theta$ changes little between outer iterations.

Fixed Parameters For fixed parameters, we do not compute the corresponding derivative. This saves the cost of unwrapping the linear solve as in Equation (15) and saves the memory to form the loss gradients, which are matrices for $P,A,C$.

Active Set Refinement See the discussion in Appendix B.

Solvers Throughout this work, we use a number of solvers available in qpsolvers including Clarabel (Goulart & Chen, ),
DAQP (Arnström et al., ),
Gurobi (Gurobi Optimization, LLC, ),
HiGHS (Huangfu & Hall, ),
HPIPM (Frison & Diehl, ),
MOSEK (Andersen & Andersen, ),
OSQP (Stellato et al., ),
PIQP (Schwan et al., ),
ProxQP (Bambade et al., ),
QPALM (Hermans et al., ),
qpSWIFT (Pandala et al., ),
quadprog (Goldfarb & Idnani, ), and
SCS (O’Donoghue et al., ).

## Appendix F Validating Algorithm 1

Figure: Figure 8: Forward and backward evaluation on the MPC dataset with the dual output by the QP solver (red) and the one obtained by solving 6.
Refer to caption: x11.png

The solvers currently available in qpsolvers provide dual solutions. Thus, to validate our modular algorithm which does not require them, we repeat the experiment in Figure 4(a). We ignore the dual solution received from the forward solver, and instead perform the optional step of computing them from the reduced KKT system 6. The results are shown in Figure 8. The additional computation of the duals has a small effect on the total backward time, as we prefactorize $K_{J}$ and use it for the derivatives as well. Using the reduced KKT to solve for the duals also impacts the duality gap, which can be seen for the larger problems in the MPC dataset, but still respect the absolute tolerance set on the duality gap.

## Appendix G Experimental Details

For completeness and reproducibility, we include additional details on the experiments. We run all experiments and methods on CPU, including methods that support GPU such as OptNet.

### G.1 Performance Evaluation

All experiments in this section were run on a Macbook Air with Apple M2 chips, 8 cores, and 16GB RAM.

In our QP benchmark experiments, we evaluate the solution accuracy using the primal residual $r_{p}$ (the maximum error on equality and inequality constraints), dual residual $r_{d}$ (the maximum error on the dual feasibility condition), and duality gap $r_{g}$ (the difference between primal and dual optimal values).

$$ $\displaystyle r_{p}$ $\displaystyle=\max\left(\left\|Az-b\right\|_{\infty},\left[Cz-d\right]_{+}\right)$ $\displaystyle r_{d}$ $\displaystyle=\left\|Pz+q+A^{T}\lambda+C^{T}\mu\right\|_{\infty}$ $\displaystyle r_{g}$ $\displaystyle=|z^{T}Pz+q^{T}z+b^{T}\lambda+d^{T}\mu|$ $$

Throughout our experiments, we present results for the duality gap to indicate the solution accuracy since, for a strongly convex QP, a zero duality gap $r_{g}=0$ is a necessary and sufficient condition for optimality.

For the forward, we set the absolute residual tolerance to $10^{-6}$.
We set the active constraint tolerance to $\epsilon_{J}=10^{-5}$. We run each problem separately with batch size 1.

In our benchmark, we regard a problem as successfully solved if it meets the following criteria:

- 1.
The solve time is less than a practical 800s time limit.
- 2.
The primal residual, dual residual, and duality gap are less than 1.0. This is a coarse check, less stringent than the imposed tolerances.
- 3.
The differentiation is executed, and does not lead to a fatal error (e.g. due to non-invertibility of a linear system).

Experimental results are averaged over 5 independent samples.

Since SCQPTH does not support equality constraints, we convert them into an corresponding set of inequality constraints.

#### G.1.1 Random dense/sparse problems

We generated two sets of random QPs: dense and sparse. For the dense set, the data is generated as $P=Q^{T}Q+10^{-4}I$, where $Q\in\mathbb{R}^{n\times n}$ with $Q_{ij}\sim\mathcal{U}(0,1)$, $C\in\mathbb{R}^{m\times n}$ with $C_{ij}\sim\mathcal{U}(0,1)$, $d=C\mathds{1}+\mathds{1}$, and $A\in\mathbb{R}^{p\times n}$ with $A_{ij}\sim\mathcal{U}(0,1)$, $b=A\mathds{1}$. The set, with 450 problems, has dimensions $n\in$ {10, 20, 50, 100, 220, 450, 1000, 2100, 4600}, $m=n$, and $p=n/2$. Each dimension contains 50 problems. We use DAQP and ProxQP as the forward solvers. Figure 9 shows that our method is comparable to OptNet and QPLayer in both time and accuracy. For smaller dimensions ($n\leq 1000$), DAQP provides higher accuracy, while for larger problems, ProxQP is more efficient.

For the sparse set, $P$ is generated as $P=L^{T}L$, where $L$ is the standard Laplacian matrix of $k$-nearest graph ($k=3$). Entries of $C$ and $A$ are filled by $\mathcal{N}(0,1)$ random numbers with density of $5\times 10^{-4}$ and zero row is avoided. The vectors $d$ and $b$ are generated similarly to the dense set. The set, with 625 problems, has dimensions $n\in$ {100, 220, 450, 1000, 2100, 4600, 10000}, with $m=n$ and $p=n/2$. For $n\leq 4600$, each dimension contains 100 problems and 25 problems for $n>4600$. KKT systems in these problems tend to be ill-conditioned. We use Gurobi as the forward solver and employ least squares solver for backward. The numerical results indicate that OptNet fails on all problems, and SCQPTH is the slowest and fails when $n\geq 4600$. Thus we excluded them in Figure 9 , which demonstrates our superior accuracy and efficiency over QPLayer.

Figure: Figure 9: Time and accuracy performance of random dense/sparse problems
Refer to caption: x12.png

#### G.1.2 Projection onto the probability simplex

Figure: Figure 10: Projection onto the probability simplex: Forward/ Backward time, referred to Figure 1
Refer to caption: x13.png

Figure: Figure 11: Projection onto chain
Refer to caption: x14.png

This formulation projects a vector onto the probability simplex, as formulated in $P_{1}$. We set $x\in\mathbb{R}^{n}$ with $x_{i}\sim\mathcal{N}(0,1)$. The set, with 500 problems, has dimensions $n\in$ {10, 20, 50, 100, 220, 450, 1000, 2100, 4600, 10000, 100000}. For $n\leq 4600$, each dimension contains 50 problems and 25 problems for $n>4600$. Gurobi serves as our forward sparse solver. Fig 1 shows the median performance within the $1/4$ and $3/4$ quantiles for each dimension. As for SCQPTH, excluded from the figure, failed at all problems when $n>50$. Figure 10 shows that we outperform in both forward and backward times.

#### G.1.3 Projection onto chain

As formulated in $P_{2}$, this experiment projects the input point cloud $x_{1},...,x_{m}\in\mathbb{R}^{d}$ onto a chain with link of length bounded by 1 in $\infty$-norm. We set $x_{i}\sim\mathcal{N}(0,100I_{d})$, with the number of points $m=100$. By varying the dimension of the vector, $d$, we generated 300 problems in this set with dimensions $n\in$ {200,500,1000,2000,4000,10000,100000}. For $n\leq 4000$, each dimension contains 50 problems and 25 problems for $n>4000$. Gurobi was used as our forward solver. Figure 11 shows performance similar to that shown in Figure 1 in terms of efficiency. In addition, dQP successfully solves large-scale problems other solves fail to solve.

### G.2 Sudoku

The Sudoku experiment was run on an Intel(R) Core(TM) i7-8850H CPU @ 2.60GHz with 6 cores and 16GB RAM.

The set-up of the Sudoku problem is a perturbed linear program

$$ $\displaystyle z^{*}(q;A,b)=\operatorname*{arg\,min}_{z}$ $\displaystyle\alpha z^{T}z+q^{T}z$ (16) subject to $\displaystyle Az=b$ $\displaystyle z\geq 0,$ $$

where $q$ encodes the input unsolved board and $z^{*}(q)$ encodes the solved board. We distinguish $q$ from the other input data, the constraints $A,b$, which model the Sudoku rules are learnable parameters that are optimized by minimizing the mean squared error to the ground truth solution for training boards. Instead of treating $A,b$ as completely independent, they are parameterized to ensure feasibility. The perturbation $\alpha=0.1$ makes the problem amenable to differentiable quadratic programming.

Our reproduction of the 2x2 Sudoku experiment from OptNet follows closely with their original settings (Amos & Kolter, ), except that we use batch size one, run exclusively on CPU, and modify the solution tolerances. For OptNet and dQP, we use the same solution tolerance $\epsilon_{\textrm{abs}}=10^{-6}$, and for dQP, we use the active tolerance $\epsilon_{J}=10^{-5}$. We use the optimizer Adam with learning rate $10^{-3}$ for both methods for the training over 10000 samples, split into 9000 training and 1000 test samples (Kingma & Ba, ).

### G.3 Bi-Level Geometry Optimization

Figure: Figure 12: The evolution of the loss for the mappings of the square into the cross. Iterations for which the active set change are denoted with vertical red lines. (a) Without regularization, the loss is driven monotonically to the tolerance. (b) With a competing regularizing loss term (dashed) convergence to the tolerance is slowed but not prevented.
Refer to caption: x15.png

The geometry experiments were run on an Intel(R) Core(TM) i7-8850H CPU @ 2.60GHz with 6 cores.

The cross and ant (Figure 2) meshes and boundary constraints are obtained from the datasets in (Du et al., ). We create the mesh refinement example in Figure 5 by perturbing the corner of a square mesh. Importantly, all of the boundary maps selected in our experiments are free of self-intersections, so that preventing triangle inversions implies the global bijectivity of the maps.

To optimize over Laplacians $M$, we directly parameterize the space of Laplacians; we impose that the diagonals are the absolute row sums during optimization and that the off-diagonals are negative. We also constrain $M$ to have the same sparsity pattern as the combinatorial Laplacian $M_{c}$. We note that the original conditions of (Kovalsky et al., ) were formulated in terms of negative semi-definite Laplacians, and so we must transform the problem into the standard form 1. Since the Laplacian $M$ which takes the place of the quadratic term in 1 has a trivial eigenvalue, the resulting QP does not have strict convexity. To address this, we perturb $M$ by a small scaling of the identity $10^{-4}I$.

Throughout the geometry experiments, we use the same solution tolerance $\epsilon_{\textrm{abs}}=10^{-5}$ and active tolerance $\epsilon_{J}=10^{-4}$ with the forward solver PIQP as determined in Appendix D. For the outer optimization, we use the optimizer Adam with learning rate $10^{-2}$ (Kingma & Ba, ). We initialize the bi-level optimization with $M_{c}$. The optimization for the cross experiment is shown in Figure 12(a) where the unregularized loss is driven to the desired tolerance, accompanied by sudden changes in the active set as the dual variables are driven to zero. We terminate the optimization at convergence, once all constraints are inactive to guarantee a bijective map. For the regularized optimization (Figure 12(b)), we penalize deviations from the initial combinatorial Laplacian up to a rescaling using the regularization $\lambda\|\frac{M}{\|M\|_{F}}-\frac{M_{c}}{\|M_{c}\|_{F}}\|_{\infty}$. In the cross shown in section 4, we choose the regularization hyper-parameter to be $\lambda=10$ after sample testing. This regularization is initially weak and so the duals are driven down, eventually crossing the regularization loss as it increases. Yet, while this slows convergence, it does not prevent it – crucial to reach a bijective map because the conditions in (Kovalsky et al., ) require all of the inequality constraints to be inactive.

The backward timing that we report in Figure 5 is for the backpropogation through only dQP, as described in Appendix 15. Thus, we remove any contribution coming from the set-up of the parameterized Laplacian and directly report the time to solve the reduced KKT and extract the gradients with respect to $M$.

## Peer Reviews

confidence: 4
contribution: 2
initial_rating: 5
presentation: 3
questions: See weakness.
soundness: 2
strengths: 1. The paper achieves a commendable integration and decoupling of algorithms. Previous works in differentiable quadratic optimization often entangled forward and backward algorithms, or required backward algorithms to utilize dual variables, which were only provided by a few solvers (like gurobi and cplex). By leveraging the KKT conditions to recover dual variables and offering a corresponding backpropagation method, this paper resolves the coupling issue present in earlier approaches. 

2. By addressing the coupling of forward and backward processes, the paper enables the selection of optimal solvers to enhance forward performance in the optimization layer. The extensive experiments presented demonstrate improved success rates and times in the forward process when using existing differentiable optimizers such as QPTH and SCQPTH. Thus, the findings of this paper can be applied to larger problem sizes.
summary: This paper addresses a common issue in Quadratic Programming (QP) by proposing a decoupled algorithm for computing gradients. The core principle is to utilize the KKT (Karush-Kuhn-Tucker) conditions to recover the dual variables $\lambda$, $\mu$, given only the primal variables. Additionally, the paper employs the active set method, allowing for the reconstruction of an equivalent new problem based on the provided primal solution $z^*$, which includes only the original equality constraints and the active inequality constraints. Consequently, this work presents a generalized scheme for backpropagation in computational optimization layers, expanding the options for solvers.
weaknesses: 1. While this work builds upon established methods effectively, there may be room for further originality in certain aspects. Though it successfully resolves the coupling issue in related works, it is limited to the specific case of QP. The recovery of KKT conditions in QP simplifies to solving a linear equation; however, in more general cases, altering the function (f) could lead to a non-linear problem, necessitating solvers that handle such complexities, which may incur additional computational costs. 

2. There are concerns regarding the fairness of the experimental comparisons. Although the paper demonstrates faster solver speeds due to the ability to select optimal solvers, it should be noted that these optimal solvers (e.g., Gurobi) are implemented in C/C++, while QPTH and SCQPTH are implemented in Python. Given the inherent speed differences between these languages, the observed performance advantage may stem from the implementation language rather than the algorithms themselves. It would be beneficial for the authors to present additional experimental data showcasing the time taken by QPTH, SCQPTH, and dQP during both the forward and backward processes. A separate demonstration using QPTH/SCQPTH for the forward process combined with dQP for the backward process would further ensure a fair comparison. 

3. The paper mentions a lack of GPU support. From the perspective of the backward process, the proposed backpropagation algorithm appears straightforward to apply to batches of data, similar to the capabilities of OptNet in this context. However, the absence of GPU support seems to arise from the external solvers utilized, which do not support GPU implementation. Given the significance of batch processing on GPU for deep learning training, it may be advantageous to develop Python versions of some solvers to enable batch mode support.

---

confidence: 4
contribution: 2
initial_rating: 5
presentation: 3
questions: - I am confused regarding the authors interpretation of explicit vs implicit, which is highlighted many times throughout the manuscript. If I understand correctly, the authors use the term explicit to refer to the fact that by using a factorization of the KKT matrix inverse (e.g. eq (6) and (7)) no optimization procedure has to be used to obtain the result on the left hand side. But why is this in any way fundamentally different to the (according to the authors implicit) differentiation in eq. (3) by multiplying with the inverse of the KKT matrix and solving using a factorization? After all, eq. (7) is exactly the result obtained from the *implicit* function theorem applied to an equality-constrained QP?
- The factorization used to solve eq (6) and (7) depends on the active set $J$, which means we cannot re-use it if we encounter a different active set by changing only the vectors $(q,b,d)$ while keeping $(P, A, C)$ constant. This is in contrast to computing a factorization of the whole KKT matrix as in OptNet, which can be re-used for any value of $(q,b,d)$. Have the authors thought about this potential drawback? I assume it could be especially relevant in some applications where $(P, A, C)$ is considered constant throughtout the training procedure and only $(q,b,d)$ is learned. Is it e.g. possible to re-use the factorization of the KKT matrix with $(P, A, C_J)$ to more efficiently compute the factorization of a KKT matrix for a different active set $(P, A, C_{J'})$?
- The authors write: "Furthermore, our approach not only simplifies computation but also allows for the use of fast, specialized linear solvers that exploit the problem’s structure" (line 322) Could the authors elaborate on this? Is this already used in the software package?
soundness: 3
strengths: - dQP as a programming package seems very useful, the ability to choose between any state-of-the-art solver is a strong selling point.
- Re-using the matrix factorization that from the derivative computation for calculating the optimal dual variables from the primal variables is useful insight.
- The bi-level geometry optimization problem is an interesting new application of differentiable optimization methods.
summary: This paper proposes dQP, a framework for differentiation of blackbox quadratic programming solvers. The method is based on the observation that an inequality-constrained QP can be re-phrased to a locally equivalent equality-constraint QP given the active set at the solution. This eliminates the inactive inequality constraints from the differential KKT matrix, which is substantially reduced in dimension and becomes symmetric, allowing for faster derivative computation. The required factorization can also be used to compute dual variables from the primal variables, enabling support for solvers that only produce primal variables as outputs.
Experimentally, the authors test their method on a diverse benchmark of QP problems. Especially when the QPs are structured, dQP inherits the benefits of the superiority of the enabled blackbox solvers compared to the simpler dense solvers in existing differentiable QP layers, which allows dQP to be applied to much larger sparse QPs.
weaknesses: - The theoretical contribution of the paper is minor. The only insight is that the optimal dual variables can be obtained by re-using the KKT matrix factorization.
- The experiments do not sufficiently disentangle forward and backward computation time. Whenever the forward+backward time is reported, the authors should also report the the forward and backward time individually. It is clear that the solvers used in dQP can solve larger problems faster on the forward pass, but it is also important to see how the required compute for the derivative computation compare to exisitng differentiable QP layers. 
- Currently, the "modularity and performance" and "scalability" experiments feel like simply benchmarking the different solvers, the gradient computation time is reported but no actual learning takes place. The Sudoku experiment does not really show any advantage over the previous OptNet method (also, here the training times are missing to jusge how big the disadvantage of not allowing GPU parallelization in dQP is). The paper would be much stronger if these results would be moved to the appendix, and more results of the kind of "Bi-Level Geometry Optimization" would be included instead.
- dQP does not allow GPU parallelization, as opposed to previous methods.

---

confidence: 4
contribution: 2
initial_rating: 5
presentation: 3
questions: - It is mentioned that gradients are obtained 'explicitly' from the primal solution. How?
- It is frequently promoted that something is 'explicit' ("Notably, our explicit perspective recasts the traditional implicit differentiation approach into an explicit method" l.87, "derivatives can be explicitly derived" l.272, "optimal point can be explicitly differentiated" l.301, "it provides an explicit expression for both the primal-dual optimal point" l.293). I don't see how the formulas are more explicit than those in the implicit function theorem -- It requires solving a linear system. I recommend avoiding this promotion.
- The paper claims that it is the first method of this kind. For instance, "LPGD: A General Framework for Backpropagation through Embedded Optimization Layers" unifies several methods and is applicable also to QP (they also do the Sudoku experiment)


- l.15 “This has so far prevented the use of state-of-the-art black-box numerical solvers within neural networks, as they lack a differentiable interface.” There are several solutions to this problem, mainly of LP, see the LPGD paper and the references therein.
- l. 244 “ensuring that a small perturbation of the _solution_ does not alter the active set” shouldn't it be 'parameters'?
- l. 250 "Implicit differentiation of these yields..." What is implicit here? It's just taking the derivatives of the equations and the chain rule.
- l.312"...the basic Sensitivity Theorem allows one to bypass the need for implicit differentiation techniques." If I understand correctly, the Theorem is based on the implicit function theorem. Hence, there is no bypass, in fact.
- l.322 I don't understand the use of 'linear solvers to exploit the problem's structure'.
- Fig 3 is not referred to in the text, and I don't understand its purpose.
- I think the overall idea of dQP is repeated too many times, and then the place for more valuable things is missing. It also feels like marketing in some places.
- Section "Approach" may be called standardly "Method."
- Related work might also mention "Unrolling methods."
soundness: 2
strengths: - The material is presented with care, and it is easy to read and follow. The problem is clearly described in great detail, and I believe the text is accessible to a broad audience.
- The idea of an independent and universal backward pass for any QP is interesting both theoretically and practically, and this paper contains a valid contribution to this topic.
- Although the idea of solving the differentiated KKT system with removed inactive constraints already appeared in prior works, this paper systematically describes and formalizes it.
- Theorem 1 allows us to compute the backward pass efficiently solely from the forward primal solution. This seems to be a new contribution.
summary: The paper belongs to the area that focuses on embedding optimization problems into neural networks. Here, optimization algorithms are integrated directly within neural networks, allowing the optimization problem to be treated as a layer with gradients that can propagate through it. This enables end-to-end learning, where neural networks learn not only representations but also parameters of embedded optimization problems.
The paper aims to integrate any quadratic programming solver in a black-box manner without specific requirements (like the availability of a dual solution), i.e., it implements the backward pass for any optimizer, using only its primal solution. The method is based on differentiating the KKT optimality conditions and solving the reduced system, where the unnecessary inactive constraints are removed.
weaknesses: - Although I personally accept the motivation that drives this work, the presented arguments should be backed by some evidence:
	- At the top of page 2, it is claimed that "many optimizers" do not provide dual solutions without mentioning any. The listed ones (Gurobi and MOSEK) provide them for QP problems.
	- There, it is also claimed that many existing methods are coupled to a tailor-made solver. Please name them.
	- At the end of related work/differentiated programming, it is claimed that the mentioned frameworks often offer limited support for solvers and may require access to a dual solution. I'm unfamiliar with all the details of these methods, and it would be beneficial if the paper could list the limitations explicitly.
- Some fragments of the theory and algorithm are missing. The overall idea is (in my opinion) clear and known. The strengths could have been in filling in the important details.
	- Computing the dual solution is said to be optional (step 4 in Algo. 1). I don't see how the derivatives (step 5) are computed if step 4 is skipped -- Eq. 7 contains both primal and dual solution.
	- Factorization (l. 320) and prefactorization (l.317, l.330) are mentioned, but never explained
- Unfortunately, any quantitative evaluation of the method is missing. All three presented experiments demonstrate and compare the performance of a selected solver. The proposed backward pass (Alg. 1) is not evaluated and I cannot confirm the proposed computational benefits of dQP.
	1. (Modularity and Performance) The "best-performing solver for each problem" with dQP backward is shown to be faster and more accurate than certain methods with their fixed solvers dedicated to certain problems. The backward pass times are not shown in separation. It it just a solver benchmark (and moreover not a fair one; OptNet was not run on GPU, which is its strength.)
	2. (Scalability) Also, OptNet was designed for lot of small dense problems and has a parallel implementation on GPU. The comparison on large sparse problems on CPU is not much valuable.
	3. (Sudoku) Loss over epoch is not informative since both methods should do the same thing (GD using analytic gradients). Training over time can tell more if dQP would prove to have more efficient backward. However (as noted) dQP suffers from infeasible solutions and OptNet is faster due to batched GPU implementation.
	4. (Bi-Level Geometry Optimization) It is just a solver comparison of forward times, not related to dQP.

Overall, the problem is well-motivated, and the approach is theoretically and practically promising. While the novelty of deriving the backward pass from the forward primal solution alone is valuable, gaps in the theoretical explanation, missing evidence to support the claims about limitations in existing methods, as well as the lack of a quantitative evaluation directly comparing the proposed backward pass with other methods, limit the paper's impact.

---

confidence: 5
contribution: 2
initial_rating: 3
presentation: 2
questions: 1. What are the benefits of dQP over OptNet? Is it more efficient? If so, why? Intuitively, it seems that compared to existing methods, the number of computation steps hasn’t decreased—rather, additional steps related to tolerance testing have been introduced. For example, line 5 in Algorithm 1 appears more complicated.
   
2. How the tolerance $\epsilon_J$ is set? Couldn’t this lead to infeasibility? I believe the authors should analyze the impact of this parameter in the experiments.
   
3. Regarding the duality gap, why does OptNet show a larger duality gap in Figure 4, even though the forward method is the same for both OptNet and dQP?
   
4. The Bi-Level Geometry Optimization in the last part of the experiment is not a common setup. I’m curious about how it can be formulated into KKT conditions. Could you provide more details on this?
   
5. How is Sparse QP defined in the paper, and what type of problems qualify as sparse? OptNet has a "lsqr" mode that significantly accelerates performance when handling sparse optimization problems. It seems the authors did not correctly utilize the different modes of the existing methods.

6. Besides, I think the authors should elaborate more on why dQP is particularly effective for large-scale problems. What mechanism specifically helps address this challenge?

If my question are well-addressed, I'll be happy to raise my score.
soundness: 2
strengths: The authors present a solid idea aimed at improving the efficiency of the backward process, making it especially effective for large-scale problems. 

They conducted numerous experiments and provided a theoretical proof to support their methods.
summary: This paper presents dQP, a framework for differentiating QP solvers, and conducted the experiments through blackbox solvers for comparing with many existing solvers.The authors demonstrate that dQP outperforms existing models across a range of benchmark problems.
weaknesses: I believe the authors do not clearly address their contribution. If the goal of this paper is to propose a generic framework, then it should be clarified how dQP incorporates other solvers. On the other hand, if dQP is intended as a specific method of differentiate opt layer, it should be compared to various state-of-the-art approaches to demonstrate its advantages.


From my understanding, dQP seems to be a novel method for differentiating optimization problems. In Eq (10), it seems all the inequality constraints are tightened into equality constraints. How does this approach ensure the feasibility of the problem? Wouldn’t this compromise the feasibility, potentially resulting in the optimal solution of the modified optimization problem being different from the original one? If that's the case, this framework may not be applicable to predict-then-optimize problems because the optimal solution of the optimization problem has been altered. 

However, if this is not the case and a black-box solver is used to solve the problem, then essentially a new forward process is designed for backward optimization but is not actually used, there will be two forward processes, which makes this step seem redundant. Additionally, existing solvers such as OptNet don’t seem to encounter this issue. If this is the main contribution of the paper, I think this it may not be sufficient.


Besides, I’m concerned about the practical performance of dQP (maybe a inexact gradient approach?) I would like to see experimental results comparing dQP with many other existing methods, especially from the aspect of computational efficiency. In the current manuscript, only graphical results are provided, such as in Figure 5, where for around 100 nodes, the results all fall within 1 second, showing minimal differences. Including more numerical results, such as the similarity of gradients computed by different solvers, running times, and differences in the optimal solutions of the optimization problems, and so on, could make the results more convincing.


One more issue is that, in the abstract (and line 355), the authors wrote, "Our implementation, which will be made publicly available, interfaces with an existing framework that supports over 15 state-of-the-art QP solvers." Maybe I think the implementation hasn’t been fully prepared before submission, so the paper may not be fully ready for publication, as it seems there is still some unfinished work.

Meta-review: Summary: The paper presents dQP, a framework for differentiating blackbox quadratic programming solvers, which allows for the integration of various solvers into neural networks without specific requirements like providing dual solutions. The method leverages KKT conditions to recover dual variables and presents a generalized scheme for backpropagation in computational optimization layers. The experiments indicate that dQP can handle larger sparse problems more efficiently than existing differentiable QP layers.

Strengths:

dQP offers a novel approach to integrating optimization problems into neural networks by differentiating blackbox solvers, which expands the range of available solvers for such tasks.

The framework's ability to work with only primal solutions is a significant advantage, as it broadens the applicability of differentiable optimization layers to solvers that do not provide dual solutions.

Drawbacks:

The paper lacks a clear comparison to state-of-the-art methods, particularly cvxpy, which is known for its efficiency in differentiating sparse problems.

The experiments do not sufficiently demonstrate the computational efficiency of dQP, with some results showing minimal differences or even underperformance compared to existing solvers.

There is a concern about the practical performance of dQP, especially regarding the potential impact of the tolerance parameter on feasibility and the lack of GPU support, which is crucial for deep learning applications.

Based on the above points , I must reject this work due to the insufficient demonstration of its advantages over existing methods, the lack of clarity on its contributions, and the concerns regarding its practical performance and applicability.

Decision: Reject