Differentiable Optimization of Generalized Nondecomposable Functions using Linear Programs
Authors: Zihang Meng
Lopamudra Mukherjee
Vikas Singh
Sathya N. Ravi
Venue:  (2021)
Keywords: linear programming
nondecomposable functions
differentiable
AUC
Fscore

Abstract: We propose a framework which makes it feasible to directly train deep neural networks with respect to popular families of task-specific non-decomposable per- formance measures such as AUC, multi-class AUC, F -measure and others, as well as models such as non-negative matrix factorization. A common feature of the optimization model that emerges from these tasks is that it involves solving a Linear Programs (LP) during training where representations learned by upstream layers influence the constraints. The constraint matrix is not only large but the constraints are also modified at each iteration. We show how adopting a set of influential ideas proposed by Mangasarian for 1-norm SVMs – which advocates for solving LPs with a generalized Newton method – provides a simple and effective solution. In particular, this strategy needs little unrolling, which makes it more efficient during backward pass. While a number of specialized algorithms have been proposed for the models that we de- scribe here, our module turns out to be applicable without any specific adjustments or relaxations. We describe each use case, study its properties and demonstrate the efficacy of the approach over alternatives which use surrogate lower bounds and often, specialized optimization schemes. Frequently, we achieve superior computational behavior and performance improvements on common datasets used in the literature.


Under review as a conference paper at ICLR 2021
DIFFERENTIABLE OPTIMIZATION OF GENERALIZED NONDECOMPOSABLE FUNCTIONS USING LINEAR PRO- GRAMS
Anonymous authors Paper under double-blind review
ABSTRACT
We propose framework which makes it feasible to directly train deep neural networks with respect to popular families of task-specific non-decomposable per- formance measures such as AUC, multi-class AUC, F-measure and others. A common feature of the optimization model that emerges from these tasks is that it involves solving Linear Programs (LP) during training where representations learned by upstream layers influence the constraints. The constraint matrix is not only large but the constraints are also modified at each iteration. We show how adopting a set of influential ideas proposed by Mangasarian for 1-norm SVMs which advocates for solving LPs with generalized Newton method - provides a simple and effective solution. In particular, this strategy needs little unrolling, which makes it more efficient during backward pass. While a number of specialized algorithms have been proposed for the models that describe here, our module turns out to be applicable without any specific adjustments or relaxations. We descrihe each use case, study its properties and demonstrate the efficacy of the approach over alternatives which use surrogate lower bounds and often, specialized optimization schemes. Frequently, achieve superior computational behavior and performance improvements on common datasets used in the literature.
INTRODUCTION
Commonly used losses such as cross-entropy used in deep neural network (DNN) models can be expressed as sum over the per-sample losses incurred by the current estimate of the model. This allows the direct use of mature optimization routines, and is sufficient for a majority of use cases. But in various applications ranging from ranking/retrieval systems to class imbalanced learning, the most suitable losses for the task do not admit a "decompose over samples" form. Examples include Area under the ROC curve (AUC), multi-class variants of AUC, F-score, Precision at fixed recall (P@R R) and others. Optimizing such measures in a scalable manner can pose challenges even in the shallow setting. Since the compromise involved in falling back decomposable loss when non-decomposable objective may be more appropriate for the task at hand can range from negligible to concerning depending the application, the last few years have approaches proposed to efficiently deal with such structured and non-decomposable losses. To this end, recent algorithms for AUC maximization have been developed based on convex surrogate losses Liu et al. (2018); Natole et al. (2018) in a linear model or in conjuction with a deep neural network Liu al. (2019) as well stochastic and online variations (Ataman et al. (2006); Cortes & Mohri (2004); Gao et al. (2013); Liu et al. (2018; 2019)) are available. Methods for measures other than the AUC have also been studied exact algorithms for optimizing F-score Nan et al. (2012); Dembezynski et al. (2011) 1), optimizing average precision through direct optimization Song et al. (2016), scalable methods for non-decomposable objectives Eban al. (2017); Venkatesh et al. (2019) and using a structured hinge-loss upper bound to optimize average precision and NDCG Mohapatra et al. (2018). Recently, the AP-Perf method Fathony & Kolter (2020) showed how custom non-decomposable performance measures be conveniently incorporated into differentiable pipelines.
number of interesting
is known that a number of these -on-decomposable objectives can be expressed in the form of an integer program that can be relaxed to linear program (LP). Earlier approaches adapted methods for structured SVMs Joachims et al. (2009) or cutting plane techniques Yue et al. (2007) and were interesting but had difficulty scaling to larger datasets. More recently, strategies have instead focused

---

Under review as a conference paper at ICLR 2021
on the stochastic setting where we operate mini-batch of samples. A common strategy, which is generally efficient, is to study the combinatorial form of one or more losses of interest and derive surrogate lower bounds. These are then tackled via specialized optimization routines. The arguably more direct alternative of optimizating as a module embedded in a DNN architecture remained difficult until recently but OptNet Amos & Kolter (2017) and CVXPY now offer support for solving certain objectives as differentiable layers within a network. Further, methods based on implicit differentiation have been applied to various problems, showing impressive performance.
Our approach is based on the premise that tackling the LP form of the non-decomposable objective as a module within the DNN, one which permits forward and reverse mode differentiation and can utilize in- buili support for specialized GPU hardware in modern libraries such as PyTorch, is desirable. First, as long as a suitable I.P formulation for an objective available, the module may be directly used. Second, based on which scheme is used to solve the LP, one may be able to provide guarantees for the non-decomposabll objective based on simple calculations (e.g.. number of constraints, primal-dual gap). The current tools, however, do not entirely address all these requirements, we briefly describe next. A specific characteristic of the LPs that arise from the losses mentioned above is that the constraints modified at each iteration data in the upstream layers. Further, the mini-batch of samples changes at each iteration. Solvers within CVXPY, effective but due their OptNet is quite efficient but designed for quadratic programs (QP): the theoretical results and its efficiency depends on factorizing a matrix in the quadratic term in the objective (which is zero/non- invertible for LPs). The primal-dual properties and implicit differentiation for QPs do not easily translate LPs, and efficient ways of dealing with constraints that are iteratively updated are not widely available at this time. In principle, of course, backpropagating through a convex optimization model (and in particular, LPs) not an unsolved problem. For LPs, we can take derivatives of the optimal value (or the optimal solution) of the model with respect to the LP parameters, and this can be accomplished by calling powerful external solver. Often, this would involve running the solver the CPU, which introduces overhead. The ideas in Meng et al. (2020) are relevant in that the optimization steps for the LP are unrolled and only involve simple linear algebra operations but the formulation is applicable when the number of constraints are about the same as the number of variables- an assumption that does not hold for the models we will study.
function of the updates to the representations of the
nature, rely on interior point methods.
In $3, we show that the modified Newton's 's algorithm in Mangasarian (2004) can be used for deep neural network (DNN) training in an end-to-end manner without requiring an external solvers where support for GPUs remains limited. Specifically, by exploiting self-conccrdance of the objective, we show that the algorithm can converge globally without line search strategies. On the practical side, we analyze the gradient properties, and some modifications improve stability during backpropagation. We show that this scheme based on Mangasarian's 's parametric exterior penalty formulation of the primal LP can be a computationally effective and scalable strategy to solve LPs with a large number of constraints.
2 NONDECOMPOSABLE FUNCTIONS AND CORRESPONDING LP MODELS
We first present a standard LP form and then reparameterize several generalized nondecomposable objectives in this way, summarized in Table in the appendix. We start with the binary AUC, extend it to multi-class AUC, and then later, show a ratio objective, F-score. Some other objectives are described in the appendix.
2.1 NOTATIONS AND GENERALIZED LP FORMULATION
Notations. We use the following notations:
(i) n: number of samples used in training.
(ii) X E Raxdd the explanatory features fed classifier (e.g., parameterized by w);
(iii) (xx) (or f(i)): a score function for the classifier where I; E X such that f(X)= = wX:
(iv) Y E {0, 1]: target label and Y E {0, 1): predicted label for binary classification, both in RR:;
(v) (()): non linear function applied on f(X):
(vi) A @ B: Kronecker product of matrices A and B.
(vii) I,: Identity matrix of size r and 1 is the indicator function.
(viii) Bg.- (and B),a) gives the k-th row (and k th) column of B.
2

---

Under review as a conference paper at ICLR 2021
LP formulation. We consider a general linear program (LP) that contains nonnegative variables as well as inequality and equality constraints. The form the LP is given as
max g" + h s.t Eu + Fu Sp + Gu =q M,80
U, 2 u
(I)
We can write it more succinetly as variable æ = [u u]:
coefficient c= l-g h]: constraints A
E F
B
- B T G
constants = [pq - 9
The corresponding primal LP becomes ming of x s. Ar < b, x 20
2.2 MAXIMIZING AUC
The Area under the ROC Curve (AUC) calculates the probability that a classifier f(-) will rank randomly chosen positive sample higher than randomly chosen negative sample. Since AUC varies between 0 and where represents all positives being ranked above the negatives, AUC may be estimated using the Vilcoxon-MMan-Ve Whitney (WMW) Statistic Hanley & McNeil (1982), as
Definition 2.1 (AUC). Let n be the number of samples. Let X+ (and X resp.) be the set of positive (and negative resp.) samples such that |X+1 + |X-| where | is the cardinality of the set. Then,
AUC is given as Σ Σ |X- 1 (())>>((,)))
X
[|XX||X_]] for 4; : E (1,-..,7
Here, we follow Ataman et al. (2006) to calculate the AUC based on the WMW statistic as follows.
X+I|X_
min Σ Adj
s.t.
f(x,) - f(zz) - 2i where æ; E X4, E E-=; Biy 2 0, (2)
where is given constant. The model in (2), maximizes AUC indirectly, by minimizing the number of pairs (one each from the positive and negative classes) where the positive sample not ranked higher than the negative sample: so, the number of zero entries in equals the number of pairs where this condition is not true. To compute the AUC, we do:
AUC (n ||=110) (IX- X-I) (n ΣΣ E relul( 22 + (XX (3)
If zi is 0, then e- rrelul(0 - -ig + equals 1. Otherwise 2gy > 0, follows from the first constraint in (2), that 2ij 2 e, so r0lu -zy + <) equals Observe that in (2), the number of constraints is X-| which is quadratic in
2.3 MAXIMIZING MULTI-CLAS AUC
An extension of AUC to the multi-clss. setting, AUC is defined in Kleiman & Page (2019). The , objective optimizes an indicator matrix computed on the orientation function, O,,5 defined as,
Definition 2.2 (Orientation Function; Kleiman & Page (2019)). Assume we have K classes {y1, } J. Let f(x;, -) E R indicate the model's prediction on æ; for each of the K classes (class-specific probability) Let ot provide the index of ır; 's true class label. Let P E Rkxk x be partition matrix where P & is the cost of classifying a sample as class k when it should be k. Define Vkk = P PP. and = Vztzz E R Then, Oi, = (Va; Vaz (VV f(r,, -)) (V. f -))).
Here, V&E ranks the instances by their cost difference between assignments to class k and k.. For two classes, say indexed by and 2, V1,2 ff is the decision boundary between the classes. The ranking is correct for 2; and xy (with correct labels 1 and resp.) if the orientation of the points w.r.t. the decision boundary is the same as those of the class labels (converted to one-hot vectors).
Then AUC,, Eneew Ok.k' EiED 10, where D₂ is the index set of samples having same class label and is a constant derived from class sizes. In formulating the linear program, how P is defined. One way is to set P(k, k) 0 Yk with arbitary entries or formulate AUC a
observe that the model is dictated and 1 for all other entries. We can also define one-vs-all setting. Here, we show the model in the first case (the other two cases discussed in
extenssion of f(.) multi-class setting
3

---

Under review as a conference paper at ICLR 2021
Appendix A.1). Similar to previous formulation, the goal is to minimize the sum of negative values of the orientation function, rather than maximize the sum of positive values. Let f(i, f(r;, x;) - f(Ij. ot) + f(z;, - f(I,,27) Then the LP formulation is.
A
bin
min
s.t. d,f((,)) 2 2 - Sij. Zij
<2
Zij 2 0
(4)
==1 1= <aj
Here d.j = V2; Vo. This can be seen as an extension of our binary AUC model, where the 2ij is the ranking between a pair of points defined for multiple classes.
2.4 MAXIMIZING F-SCORE
The F-score (or F-measure) representative of objectives expressed as ratios of some combination of True positives (TP). False positives (FP), True negatives (TN) and False negatives (FN). The general form of the ratio functions and formulations for other objectives is in the Appendix A.2. Specifically, F-score is defined as follows:
Definition 2.3 (F-score). F-score
2(Precision Recall) Procision++
21P FPP-
2(Y
1TY+1TÝ
The second part of the equality comes from simplying the precision TP4FF and recall TPTFN based formula Dembezynski et al. (2011). The last part is obtained by replacing TP with (Y" Y) FP with (1 Y) Y and FN with (Y)T x (1 Y) as functions of Y and Y. This leads to the following integer fractional optimization model,
TP
TP
CTY Y IT +
F-score = max
s.t. i (0, 1], i= I...., n where 2Y and b Σ Y.
(5)
To solve this, we first relax the constraint on y and reformulate the model as the following LP, by introducing two variables 2 E R" and & E R where
b8 t b and E (1,..,, .n):
14
J2 1T + bt b:
max 2, b
o(f(x))t S (1 ((f(x:))) 12 2,
(b)
(c)
Remark 1. (a) ensures the appropriate relation benween 2, t and Ý and is essentially reformulaton of the ratio objective as a linear function with fixed denominator. Y is recovered from the solution to the linear program by computing Y, = z: when t O and Y; 24 otherwise. (b) sets an upper bound for V, < 1. (c) ties the output of the previous layer (x:():) classifier score for Y, see the definition Section 2.1) as a input this layer. Assume ((.) E 11,, (ensured if d is sigmoid or tanh) is the indicator of the class label (based on sign). We want Y; 2 ((((z)) and Y, < 1+ ((f(())) If is {0,1}, 1}, these two constraints ensures that Y; = 0 when (xf(x.) S O and Y, = 1 when of o(f(x)) > 0. When Y4 is relaxed and replaced with and t, get the equivalent form (c).
i
This model imposes 4n constraints for samples. Since this is maximization, a solution to the LP, Oy, is an upper bound on the integer objective opt", and serves as the loss.
3 BACKPROPAGATION VIA FAST EXTERIOR PENALTY OPTIMIZATION
Unlike traditional feedforward networks, where the output of each layer is relatively simple (though non-linear) function of the previous layer, a LP layer must solve constrained optimization problem, therefore implementing scalable and efficient backpropagation schemes that minimizes computational overhead requires more care and as reviewed in $1, is an active topic of research. This problem is, of course, not unique to LPs and manifests in differentiable sorting Mena et al. (2018) and formulating quadratic or cone programs Amos et al. (2017). One may unroll gradient descent steps Amos et al. (2017); Goodfellow et al. (2013); Metz et al. (2016) use projections Zeng al. (2019). Recently
4

---

Under review as a conference paper at ICLR 2021
Agrawal et al. (2019) introduced a package for differentiable constrained convex programming, which includes LPs as a special case. For LPs, Meng et al. (2020) presents an unrolled scheme and Blondel et al. (2020) shows that we can differentiate through LP formulations of sorting/ranking exactly by using smooth approximations of projection steps. Berthet et al. (2020) describes an interesting approach where one computes approximate gradients through ranking/shortest path problems by stochastic perturbation techniques.
Remark 2. Some previous works Zeng et al. (2019) have considered LPs where the constraints are deterministic (for fixed input dimension), i.e., do not depend on the data X, which is different from the LPs in 222--2.4.
Note that such perturbation techniques in Berthet al. (2020) are applicable to our LPs as well. For example, the Fenchel Young losses as defined in Berthet et al. (2020) is attractive because there is no need to compute the Jacobian. From the implementation standpoint, one could simply think of the backward pass as a function given the input and output of the forward pass. But the gradient expressions of the losses is an expectation and hence may require multiple calls to LP solver in order to approximate the expectation. although parallelization and warm starts were shown to alleviate this dependency by sampling in parallel.
Consider LP with large m number of constraints in fixed dimensions (n & m). This assumption holds in all formulations in $2. This is because we assume that the architecture is fixed whereas minibatch size depends on the complexity of the task (stable gradient when noise in gradient is high). Hence, solving such LPs using off-the-shelf solvers as in Berthet et al. (2020) may slow down the training process. The strategy in Agrawal al. (2019) does offer benefits over Amos & Kolter (2017) for sparse QPs. Our strategy is to run Mangasarian's 's Newton's 's method on an exterior penalty function of the LP. There are two advantages: (i) during forward pass, quadratic local convergence of Newton's 's method indicates that unrolling the method may be reasonable choice; and further (ii) based on the relationship between dual and primal variables, and the exactness of the exterior penalty we can show that backward pass is independent of m. We will discuss both these results and some modifications to deal with discontinuous Hessian (and its inverse) that is required for Newton's 's method. A similar approach is proposed in Amos & Kolter (2017) in which Primal- Dual Interior Point methods with implicit differentiation is used for differentiation purposes. But the exterior penalty in (6) satisfies an interesting property: primal and dual solutions are related by a closed form expression which can be exploited for efficient backpropagation.
Rationale.
3.1 FORWARD PASS USING NEWTON'S S ALGORITHM
key requirement for fast automatic (forward or reverse mode) differentiation is that we can perform the forward pass efficiently. In our setting, seek to solve and backpropagate through an LP. We will focus on reverse mode differentiation since it is the most suitable for DNN training.
Given primal LP, for fixed accuracy 0, Mangasarian (2004) solves unconstrained problem,
min (v)) = 2 |Io(Ay - b)]l + ec y,
(6)
where a(.) max(: 0) represents the elementwise relu function. A modified Newton's method can be use to solve equation 6 that performs the following iterations:
+ Xd where d = ()) = (V2 g(v) + p/)-iVg(v))
In large scale settings of A, b, such Newton methods are known to perform empirically better than gradient descent Mangasarian (2006); Keerthi al. (2007). We will discuss if this holds for our purposes shortly.
Why is Newton's method applicable for minibatches? In general, the convergence of Newton's method depends strongly on initialization, i.e., we can only provide local convergence results. However, this is not the case for our problems since in our examples, either the level sets are bounded from below; or the feasible set is compact, as noted in Mangasarian (2004). There are two reasons why the result, by itself, is insufficient for our purposes: (i) it assumes that we can perform line search to satisfy Armijo condition; (ii) even with line search, the result does not provide rate of convergence. In DÑN training, such line search based convergence results can be prohibitively expensive. The main difficulty is handling the discontinuity in the Hessian. As remedy, we use self concordance of
5

---

Under review as a conference paper at ICLR 2021
(6) to guarantee global convergence of (7) iterations for the exterior penalty formulation in (6). To do SO, we first show a result (proof in Appendix) that characterizes the discrepancy between the actual Hessian of (6) and the modified one in (7) when A is randomly distributed.
Lemma 1. Assume that A is a random malrix, and fix some y E RR Then with probubility one, g in equation 6 is constant (given by H, V (yl over sufficiently small neighborhood of
Intuitively, Lemma I states that with probability one, each y has neighborhood in which the Hessian is constant. In addition, the modified Hessian is nonsingular at all points (in particular the optimal y²), and so we can then show the following global convergence result.
Theorem 2. Newton's 's method converges globally at linear rate with local quadratic convergence.
Proof: (Sketch) First, we use the fact that the objective in equation 6 is piecewise quadratic, and hence self concordant. Second, observe that the possible choice of Hessian are finite, and so we can choose p > 0 so that there is descent direction without line search, that is, there exist step size > 0 such that \Vg(z)* < O. Finally, we use Theorem 4.1.12 in Nesterov (2013) to claim the result.
D
Please see Appendix for the proof. The loss function is almost surely quadratic in the neighborhood of any y, thus intuitively suggests local quadratic convergence independent of the starting point.
Remark 3. Convergence in Thm. 2 is guaranteed under standard constraint qualification assump- tions. Linear Independent Constraint Qualification (LICQ) is satisfied for AUC, and Multi-class AUC formulations in $2. But the F-score formulation does not satisfy LICQ, hence we need safeguarding principles in the initial iterations (until iterates get close to the optimal solution).
3.2 BACKWARD PASS USING OPTIMAL DUAL VARIABLES AIDED BY UNROLLING
The advantage of optimizing the exterior penalty in (6) is that given an iterate yt. accuracy c, we can extract the optimal dual solution e by simple thresholding, that is, 1/E(A" (Ay b)). By complementarity slackness, the nonzero coordinates of us specify the set of active constraints in
Ax I b. So, given an approximate solution y such that g (y)) (Vgg)) I c, to compute the primal solution x", we solve the "active" linear system given by Ab, where denotes the active rows of A and the corresponding subvector U. Hence, backpropagation through the layer reduces to computing derivatives of Á b which is simple via automatic differentiation.
How to choose e? If we can successfully retrieve the active constraints at the optimal solution, we do not need to store the intermediate iterates at all during the forward pass (memory efficient). However, setting correctly can be tricky for arbitrary polyhedra since depends on the geometric properties such as facets and vertices that may be difficult to enumerate. One possible way to get around this is use a "burn-in" period in which we increase c slightly in each iteration (of deep network training) and backpropagate through the unrolled Newton's 's iterations during this period. Once we see that the convergence profile has stabilized, we can fix e at that value and start using the complementarity conditions and derive the active linear system À - as discussed above.
How to backpropagate through unrolled iterations? We assume that the chain rule is applicable up to this LP layer and is aL d2 (for one of the parameters A), and note that is possible to find ar. de 84
da 8n which involves the LP (either directly or using chain rule). Therefore we focus our attention on 3 layer. Indeed, unrolling each iteration in (7) is equivalent a "sublayer". So order to backpropagate we have to show the partial derivatives of each operation or step wrt to the LP parameters A, b, and c.
Our goal is to ad ø4 where d = 9 Q H and u = Vg(y). We can use the product rule to arrive at: ad = (u Q- @ Q-1) -Q----Q-Q--- -'u. To see this, note that we have used the chain rule to differentiate through the inverse in the first term. The second term is easy to compute similar to the computation of Hessian. For each of these terms we eventually have to compute 2Q On where E {c, A, b} which can also be done by another application of chain rule. Please see Appendix A.4 for empirical verification of unrolled gradient and the one provided by A -b
or dz dz
Before proceeding, we should note an issue that comes up when differentiating each step of the unrolled algorithm due to the fact that the Hessian is piecewise linear (constant) as a function of

---

Under review as a conference paper at ICLR 2021
the input to that particular layer. Here, some possible numerical approximations needed, describe below.
aQ is nondifferentiable due to the presence of Remark 4. Note that the diagonal matrix term in 2A the step function. However, the step function is a piecewise constant function, and hence has zero derivative almost surely, that is, in any bounded set S, I E S. if a ball (of radius r 0,) (x) S S, then the Lebesgue measure of the set of nondifferentiable points on zero. Please see Appendix B.3 for a formal justification where show this by approximating the step function using sequence of logistic functions with increasing slope parameter at the origin.
Therefore, in this setting, Remark 4 provides way to compute an approximate sub-gradient when using Newton's method based LP layers. The function is a piece-wise quadratic function and differentiable everywhere, and the inverse of the Hessian acts preconditioner.
Summary, Our forward pass involved three steps: 1. finite steps of Newton's method using which we 2. computed the dual variable by thresholding operation, and 3. finally, get the primal solution, these dual variables are first used to identify the active constraints followed by solving a linear system. In order to backpropagate through these three steps, must differentiate through each layer of our procedure including 4-1b, independent of whether we use unrolling Danskin's 's theorem. Using Danskin's 's theorem in this setting would involve differentiating through the fixed point of the Newton's 's iterations similar to (regularized) gradient descent iterations considered in the iMAML work Rajeswaran et al. (2019).
EXPERIMENTS
In this section, we conduct experiments on commonly used benchmarks to show that our framework can be used to optimize multiple different objectives within deep neural networks and lead to performance gain. We start with binary AUC optimization, and then extend to multi class AUC optimization and F-score optimization. We also show that nonnegative matrix factorization can be optimized in linear programming form in our framework.
Optimizing Binary AUC
We follow the current state-of-the-ar work on AUC optimization Liu et al. (2019) to conduct experiments on optimizing AUC score directly with deep neural networks. The baseline algorithms we compare with for binary AUC are cross-entropy loss and two algorithm (PPD-SG and PPD- AdaGrad) from Liu al. (2019).
Datasets: Cot& CIFAR1O, CIFAR100, and STL10. Cat&Dog is dataset from Kaggle which contains 25000 images of cats and dogs. 80% of the dataset used as training set and the rest 20% as test set. STL10 is inspired by the CIFAR-I0 dataset but with some modifications. Each class in STL10 has fewer labeled training examples than in CIFAR-10. We follow Liu et al. (2019) to use 19k/1k, 45k/5k, 45k/5k, 4k/lk training/validation split Cat&Dog, CIFAR10. CIFAR100, STL10 respectively.
Table 1: Binary AUC optimization results on four benchmark datasets.
Construction of imbal- anced datasets: We construct imbalanced binary classification task by using half classes as positive class and another half as negative class, and dropping samples from negative class by a certain ratio, which is reflected by the positive ratio (the ratio of the majority class to the minority class) in Table 1.
AUC(%)
Cat&Dog
CIFARIO
Positive Ratio
91%
83%
71%
50%
91%
83%
71%
50%
Cross-Entropy
67.6
74.6
85.1
87.4
65.2
73.3
78.1
83.7
PPD-SG
79.1
81.5
85.5
87.1
69.8
73.9
79.1
82.6
PPD-AdaGrad
77.3
80.6
83.7
86.3
69.7
74.1
78.4
83.1
Ours
78.6
81.3
85.6
87.8
72.5
74.4
78.3
82.7
AUC(%)
CIFARIO
STL10
Positive Ratio
91%
83%
71%
50%
91%
83%
71%
50%
Cross-Entropy
57.8
58.4
62.2
66.3
63.3
67.1
72.7
80.8
PPD--G
56.5
58.9
61.6
65.2
70.7
71.6
75.1
77.4
PPD-AdaGrad
56.2
59.0
62.6
67.6
68.5
72.4
76.7
78.5
Ours
58.2
60.5
64.5
69.0
68.4
71.1
76.7
81.6
Experimental setting. We use Resnet-18 18 He et al. (2016) as the deep neural network for all algorithms. During optimization, the batch size is set to 64. The initial learning rate is tuned in {0.1, 0.01, 0.001) and decays 2/3 at 2k, 10k, 25k-th iteration. We train 40k iterations in total. The in
7

---

Under review as a conference paper at ICLR 2021
Newton's method is 0.001. We the same random seed, learning rate and total number of iterations in all of our experiments including multi class AUC and F-score experiments.
Results. The results shown in Table 1. We can sec that our method slightly outperforms Liu et al. (2019) and outperforms cross-entropy loss by a large margin, especially on imbalanced datasets, where AUC objective shows superior- ity over cross-entropy loss.
Table Ablation study of e on Cat&Dog dataset.
Positive Ratio
91%
83%
71%
50%
Ours(c 0.1)
71.3
77.0
84.4
87.3
Ours(e = 0.01)
78.6
81.3
85.6
87.8
Ours(e 0.001)
65.9
71.3
71.8
76.1
Influence of c. We report the influence of our algorithm in Table 2, where choose Cat& Dog as an example to test different e. We can see that e = 0.1 gets slightly worse performance than <= 0.01, while e 0.001 performs much worse. To choose e, we follow the approach proposed by the value of the perturbed quadratic
Mangasarian (2004). If for two successive values of <1 function is the same, then it the least 2-norm solution the dual. Therefore, we simply choose an that satisfies this property, which is chosen be 0.01 in our experiments.
Optimizing AUC
Multiclass
Table Multiclass AUC optimization results on STL1O and CIFAR Drop rate is the proportion used when dropping samples from two of three classes.
We further demonstrate our method for optimizing mul- ticlass AUC. Similar to previous section binary AUC, we construct imbal- anced multiclass datasets by dividing datasets into 3 classes and drop samples from 2 of them and report the one- mer-ussus-all AUC (de- noted as AUCuaa) and AUChin score For STLIO, we group class - 2, 3 - 5, 6 - 9 into the three classes, and drop samples from the first two classes. For CIFR 100, we group class 0 32, 33 65, 66 - 99 into three classes, and also drop samples from the first two classes.
(U(
CIFAR R100
STL.10
Drop rate
90%
80%
60%
0%
90%
80%
60%
Cross-Entropy
54.3
59.4
62.7
63.5
66.9
68.0
74.8
81.0
Ours
58.4
59.2
64.1
65.7
72.9
72.5
75.7
82.7
AUC)
CIFAR100
STL10
Drop rate
90%
80%
60%
0%
90%
80%
60%
0%
Cross-Hntropy
55.1
60.6
65.0
64.0
68.9
69.6
75.8
82.2
Ours
61.2
66.0
67.2
76.1
74.4
77.7
84.5
Results. Results are in Table 3. addition one-versus-al all AUC metric, we also report the performance in terms of AUC,, Kleiman & Page (2019) which is specifically designed for measuring multiclass AUC and keeps nice properties of binary AUC such as being insensitive to class skew. We see that our method outperforms cross-entropy loss all four datasets and under all different skewed ratios. Specifically, the performance gain tends to be larger when the dataset becomes more imbalanced.
Class as positive clase
Class an positive class
Cless as positive class
3.0
Entrody
1.0
Cross-B Ents-Entropp
3.0
Entrady
Thu Prsstvve Rate
Stue Poaltive Rute
Thue Pesitive Rute
0.7
80
False Postttiv Rate
Talse Posttive Rate
Tatoe Posttive Ratt
(a) Positive class:
(b) Positive class:
(c) Positive class:
Figure 1: ROC curve of multiclass AUC optimization STL10 with 90% drop rate. We divide STLI0 into 3 classes and use one as positive class and other two as negative class to plot the ROC.
Optimizing F-score
Table 4: F-score on four datasets.
We show that by directly op- timizing F-score, we can achieve better perfor- mance on this than when us- ing cross entropy loss. In addition to cross entropy loss, perform evalua-
F-score(f6 Cross-Entrope CVXPY-SCS
Cat&
CIFAR10
CIFAR100
STL10
76.0
70.3
60.4
71.8
70.1
66.6
66.7
66.6
AP-Perf
65.3
66.7
66.4
67.2
Ours
77.8
72.6
63.4
72.7
8

---

Under review as a conference paper at ICLR 2021
tions with two other methods that also directly optimize the Fi-score. First, replace solver with CVXPY-SCS Agrawal et al. (2019), which is a differentiable general purpose linear programming solver, second, we perform comparisons with AP-Perf Fathony & Kolter (2020) which offers differentiable optimization of F-score using an adversarial prediction framework. The datasets and our setup to group them into classes remain the same in binary AUC section. The results in Table shows that our method generally yields improvement over cross-entropy loss in terms of F-score. Note that different from optimizing cross entropy loss, when we optimize F'scoore directly, there exists local optimal point where assigning all examples to the positive class leads to F-score of 66. 7%. We see this behavior on CIFARIO, CIFARIO0, and STL10.
CONCLUSIONS
We demonstrated that various non-decomposable objectives be optimized within deep neural networks in a differentiable way under the same general framework of LPs using a modified Newton's s algorithm proposed by Mangasarian. A number of recent papers have studied the general problem of backpropagating through convex optimization modules, and this literature provides several effective approaches although scalability remains a topic of active research. Our work complements these results and shows that the operations needed can be implemented to utilize the capabilities of modern deep learning libraries. While our experimental results suggest that promising results on binary AUC, multi-class AUC and F-score optimization within DNNs is achievable, we believe that the module may have other applications where the number of constraints are large and data data-dependent.
REFERENCES
Akshay Agrawal, Brandon Amos, Shane Barratt, Stephen Boyd, Steven Diamond, and J Zico Kolter. Differentiable convex optimization layers. In Advances in neural information processing systems, pp. 9562-9574, 2019. 5,
Brandon Amos and J Zico Kolter. Optnet: Differentiable optimization layer in neural networks arXiv preprint
1et.:0003304333 2017. 2,5
Brandon Amos, Lei Xu, and J Zico Kolter. Input convex neural networks. International Conference on Machine Learning, pp. 146-155, 2017. 4
Sanjeev Arora, Rong Ge, and Ankur Moitra. Learning topic models-going beyond svd. In 2012 IEEE 53rd anmual symposium foundations of computer science, pp. 1-10. IEEE, 2012. 14
Kaan Ataman, W Nick Street, and Yi Zhang. Learning to rank by maximizing auc with linear programming. In The 2006 IEEE International Joint Conference on Neural Network Proceedings, pp. 123-129. IEEE, 2006. 1,3
Quentin Berthet, Mathieu Blondel, Olivier Teboul, Marco Cuturi, Jean-Philippe Vert, and Francis Bach. Learning with differentiable perturbed optimizers. arXiv preprint urXiv:2002.08677 2020.
Mathieu Blondel, Olivier Tehoul, Quentin Berthet, and Josip Djolonga. Fast differentiable sorting and ranking. In International Conference on Machine Learning, 2020.
Mike Brookes. The matrix reference manual. Imperial College London, 3, 2005. 16
Edo Collins, Radhakrishna Achanta, and Sabine Susstrunk. Deep feature factorization for concept
discovery. In Proceedings of the European Conference on Computer Vision (ECCV), pp. 336-352, 2018. 13, 14
Corinna Cortes and Mehryar Mohri, Aue optimization vs. error rate minimization. In Advances in neural information processing systems, pp. 313-320, 2004.
Krzysztof J. Dembezynski, Willem Waegeman, Weiwei Cheng, and Eyke Hüllermeier. An ex- act algorithm for f-measure maximization. In J. Shawe- Taylor, R. S. Zemel, P. L. Bartlett, F. Pereira, and K. Q. Weinberger (eds.), Advances in Neural Information Processing Systems 24, pp. 1412 Curran Associates, Inc., 2011. URL http: htpppaappssss/ee/ nips. ccc cc/paper/ 4.89 1:39-annexacc---- exact-al gori hmm- for- f mea su re wwwrre-mxxinne ma..mmiiimmmmmmmmmmmmmmmmmmmmm zat.on..oioiiio pdf. 1, 1,4, 12

---

Under review as a conference paper at ICLR 2021
Elad Eban, Mariano Schain, Alan Mackey, Ariel Gordon, Ryan Rifkin, and Gal Elidan. Scalable learning of non-decomposable objectives. In Artificial Intelligence and Statistics, pp. 832-840, 2017. 1
Rizal Fathony and Zico Kolter. Ap-perf: Incorporating generie performance metrics in differentiable learning. In International Conference Artificial Intelligence and Statistics, pp. 4130-4140. PMLR, 2020. 1, 9
Wei Gao, Rong Jin, Shenghuo Zhu, and Zhi-Hua Zhou. One-pass auc optimization. In International conference on machine learning, pp. 906-914, 2013.
Alexander I Golikov and Igor E Kaporin. Inexact newton method for minimization of convex piecewise quadratic functions. In Numerical Geometry. Grid Generation and Scientific Computing.
pp. 139-155. Springer, 2019. 15
Ian Goodfellow, Mehdi Mirza, Aaron Courville, and Yoshua Bengio. Multi-prediction deep boltzmann machines. In Advances in Neural Information Processing Systems, pp. 548-556, 2013. 4
James A Hanley and Barbara J McNeil. The meaning and use of the area under receiver operating characteristic (roc) curve. Radiology. 143(():29-36) 1982. 3
Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recognition. In Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 770-778, 2016.7
Thorsten Joachims, Thomas Finley, and Chun-Nam John Yu. Cutting- -lanne training of structural svms. Machine learning, 77(1):27-59, 2009.
S Sathiya Keerthi, Vikas Sindhwani, and Olivier Chapelle. An efficient method for gradient-based adaptation of hyperparameters in svm models. In Advances in neural information processing systems, pp. 673-680, 2007. 5
Ross Kleiman and David Page. Aucu: A performance metric for multi-class machine learning models. In International Conference on Machine Learning, pp. 3439-3447, 2019.
12
Mingrui Liu, Xiaoxuan Zhang, Zaiyi Chen, Xiaoyu Wang, and Tianbao Yang. Fast stochastic auc maximization with o(1 convergeeee rate. In International Conference Machine Learning. pp. 3189-3197, 2018.
Mingrui Liu, Zhuoning Yuan, Yiming Ying, and Tianbao Yang. Stochastic auc maximization with deep neural networks. arXiv preprint 1exiiv:908.. 10831, 2019. 1,7,8
OL Mangasarian. A newton method for linear programming. Journal of Optimization Theory and Applications, 121(1):1-18, 2004. 5,8
Olvi L. Mangasarian. Exact 1 1-norm support vector machines via unconstrained convex differentiable minimization. Journal of Machine Learning Research, 70u)):1517-1530, 2006.
Gonzalo Mena, David Belanger, Scott Linderman, and Jasper Snock. Learning latent permutations with gumbel-sinkhom networks. arXiv preprint arXiv: arXiv: 180200665, 2018. 4
Zihang Meng, Sathya Ravi, and Vikas Singh. Physarum powered differentiable linear programming layers and applications. arXiv preprint arXiv:2004.14539, 2020. 2.5
Luke Metz, Ben Poole, David Pfau, and Jascha Sohl- fohl-Dickstein. Unrolled generative adversarial networks. arXiv preprint arXiv: /6/ 063 2016. 4
Pritish Mohapatra, Michal Rolinek, CV Jawahar, Vladimir Kolmogorov, and M Pawan Kumar. Efficient optimization for based loss functions. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 3693-3701, 2018.
Ye Nan, Kian Ming Chai, Wee Sun Lee, and Hai Leong Chieu. Optimizing f-measure: A tale two approaches. arXiv preprint arXiv: /206.6625, 2012.
10

---

Under review as a conference paper at ICLR 2021
Michael Natole, Yiming Ying, and Siwei Lyu. Stochastic proximal algorithms for maximization. In Interaational Conference on Machine Learning, pp. 3710-3719, 2018.
Yurii Nesterov, Introductory lectures on convex optimization: A basic course, volume 87. Springer Science & Business Media, 2013. 15
Aravind Rajeswaran, Chelsea Finn, Sham M Kakade, and Sergey Levine. Meta-learning with implicit gradients. In Advances in Neural Information Processing Systems, pp. 113-124, 2019. 7
Ben Recht, Christopher Re, Joel Tropp, and Victor Bittorf. Factoring nonnegative matrices with linear programs. In Advances in Neural Information Processing Systems, pp. 1214-1222, 2012. 14
Yang Song, Alexander Schwing, Raquel Urtasun, et al. Training deep neural networks via direct loss minimization. In International Conference Machine Learning, pp. 2169-2177, 2016.
George Trigeorgis. Konstantinos Bousmalis, Stefanos Zafeiriou, and Bjoern Schuller. A deep semi- nmf model for learning hidden representations. In International Conference on Machine Learning. pp. 1692-1700, 2014.
Sathya Ravi Abhay Venkatesh, Glenn M Fung, and Vikas Singh. Optimizing nondecomposable data dependent regularizers via lagrangian reparameterization offers significant performance and efficiency gains. arXiv preprint arXiv:1909. /2398, 2019.
Yisong Yue, Thomas Finley, Filip Radlinski, and Thorsten Joachims. A support vector method for optimizing average precision. In Proceedings of the 30th annual international ACM SIGIR conference on Research and development in information retrieval, pp. 271-278 2007.
Xiaohui Zeng, Renjie Liao, Li Gu, Yuwen Xiong, Sanja Fidler, and Raquel Urtasun. Dmm-net: Differentiable matc-ng network for video object segmentation. In Proceedings of the IEEE International Conference Computer Vision, pp. 3929-3938, 2019. 4,

---

Under review as a conference paper at ICLR 2021
A APPENDIX
A. LP FORMULATION FOR MULTI-CLASS AUC
One way to extend binary AUC multi-class is by considering multiple one-versus-all pairs. This leads us to the following formulation:
Σ AUC :mi
Zij
(8)
1= == Ang
S.t. (f(Z, - f(zz, 2 Zij
Vi 2 + xJ,,
24j 2 0
In multi-class AUC experiment, performance both one-versus-all AUC and A In addition, we can also consider the setting of AUC,, where P is set arbitarily. In this case, the exact terms in orientation function O proposed by Kleiman & Page (2019) can be written as follows:
this one- ver-versus-all AUC training loss and report
AUCarbit min
Zij
(9)
1: <aag
v())(f(z;) k(z,,k)) ((),,, - Zij
Vi.j :2 < zzj, Zij 2 0
Note that AUC has the same number of constraints and variables as AUC Once the LPs are solved, the loss function is calculated the same way as binary AUC.
A.2 FORMULATING RATIO OBJECTIVES
In this section, we study a subset of non-decomposable metrics, which are typically expressed as ratios of some combination of True Positive(TPP), False Positives(FP), True Negatives(TN) and False Negatives(FN). These can be expressed in general form 021 AurPPaa2 constants/cofficients which if set to 0, means the term is absent and not equal to zero in other cases. This formulation can used to define Fscore, F3. Jaccard, IOU and Precision at fixed recall. In the following section, describe the formulation of Fscore as representative of this approach, other metrics can be formulated similarly.
FN+ where apy
Given Y the groud truth, our goal is to compute Ý both of length n, which aligns with Y based on the specific metric. We first show how write TP, FP, TN and FN wrt to these vectors.
TP Y x & TN = (1 (1-1)xxxx
FP (1- y)" x FN (Y) x (1-
A.2.1 FORMULATING F-SCORE
The F-score or F- --mmasure is routinely used as performance metric for different types of prediction problems, including binary classification and, multi-label classification. Compared measures like error rate in binary classification and Hamming loss, it enforces better balance between performance on the minority and the majority classes, and, therefore, it is more suitable in the case of imbalanced data Dembczynski et al. (2011). F-score is defined as follows:
F-score 4) 2P(Y, &) x R(Y.
(10)
P(Y, Y) R(Y,))
where P is the measure of precision defined as
TP
P(Y, TP + FP
12

---

Under review as a conference paper at ICLR 2021
and R stands for the measure of recall, given
TP
ROY, Y)
TP + FN
Plugging this Eq(I), and replacing the formulations for TP, FP and FN(from Eq A.2) we get
2TP
2(Y x x)
2(Y x Y)
F-score(YYY), 2
(1I)
2i-1 > Zi=1 Vi 1TY + ITV
where y; refers to the ith element of Y (same for ú.). 1 represents an all one vector R Note that assume 1"Y 3 which is constant (the number of
in training. since Y is generally provided, examples in the positive class in the ground truth). We can also represent the the values of 2Y as coefficient matrix c, then the optimization problem for finding F-score be written
maximize Y
(12)
12 +
subject to
Y; E [0, 1}, i 1....., n.
A.2.2 MAXIMIZING JACCARD COEFFICIENT AND F3
The Jaccard Coefficient and Dice Index lead similar formulation F-score. The Jaccard coefficentt can be expressed as:
TP Jacc Y, Y) +
( Ei= V E 0:
(YT x Y)
y6 14
1TY + (1 Y)TY
(13)
This can be equivalently written as a linear factional program as shown in Model where = Y. d d= (1 Y) and b- y. The rest of the construction similar F-score.
Note that Fa which is defined as
P(Y,)) R(Y,Y) (Y (1 32)
(14)
PPP + R(Y,Y)
where 9 is user specified parameter (balaneing the importance of precision and recall) also permits a similar formulation. Here we simply set = (1 B)) d = =1 and b = 1²1TY.
A.2.3 MAXIMIZING PåR
We begin by defining the maximum precision at fixed minimum recall problem as
P Ra maximize
P S.R.Z a
(15)
(YTY)
maximize
s.t. YT 2 21
(16)
This is again a linear fractional objective with linear constraint. So can write it an equivalent Linear program using the same transformation where c Y, d = and b= 0. R@P P on the other hand directly leads to linear program.
A.3 OPTIMIZING NON-NEGATIVE MATRIX FACTORIZATION (NMF)
Nonnegative matrix factorization is different from the other objectives we presented in that it is primarily used in unsupervised learning and does not satify the criteria for metrie loss function. can still be formulated in generalized non-decomposable form because (i) cannot be written sum over individual samples, (ii) leads to a model where the constraints depend on learned features.
Note that purpose of discussing NMF in this context, is not to provide general purpose solver for the problem, and instead to assess whether NMF layers can serve as regularizer or a clustering module, e.g., learning more interpretable attributes, co-segmentation and a substitute for clustering, see Trigcorgis et al. (2014) and Collins et al. (2018). The following description in this section and the experimental validation in the following section is proof of principle instantiation of this idea.
13

---

Under review as a conference paper at ICLR 2021
We briefly review from Arora et al. (2012) and Recht et al. (2012), how NMF is written as LP.
We know from Arora et al. (2012) that a NMF decomposition V = FW where F is s x and W is s' X and V and W have row sum of 1, is 'separable' if the rows of W are simplicial and there is permutation matrix R E R*xs such that RF anchor words. Recht et al. (2012) proposed data-driven model where the most salient features in the data are used to express the remaining features, given as CV, where C' is of size 8 8.
I M 1" The top rows of F contains so-called
RV CV. To show that thsis factorization is possible, we need to first make F square, which is why it is zero-padded to make it size X matrix. Let p be any vector which is used the coefficient in the objective in the following model. According to Recht et al. (2012), any value for the entries of p should suffice as long as they are distinct. with distinct values. Then the LP formulation is follows:
Assuming V admits rank-r separable factorization, then V R.
1. M
min p" diag(C) s.t.
CV tr(C) = CJs S 1 Vj, Cip S CD Vij, C20
(17)
With C' in hand, W is constructed by extracting rows of V for those indices where CkA 1. Fis constructed by extracting rows of C. which correspond to k where Cyk = 1.
A.3.1 EXPERIMENTAL RESULTS ON NONNEGATIVE MATRIX FACTORIZATION
We demonstrate applicability of our strategy to nonnegative ma- trix factorization (NMF) by performing a rank & factorization on Convolutional Neural Network (CNN) activations as an example, following Collins al. (2018). Recall that the activation tensor of an image at some layer in CNN has the shape V Rexhxw where h, w are the spatial sizes and is the number of channels. We can reshape it into V (ex(h-w) and calculate rank NMF for V: V = FW. Each row W, of the resultant W Rkx(h-w) can be reshaped into a heat map of dimension h X w which highlights regions in the image that correspond the factor W,, We show an example for k 1, 2 in Fig. 2. We can see that heatmap consistently captures a meaningful part/concept the examples. Currently, our memory consumption increases quickly with here since the constraint matrix in our LP formulation of size (ce2) x O(c²). This makes our method only work for small a GPU with 11GB memory (here, we use c = 20). This scaling issue can be possibly solved by utilizing sparsity in the constraint
Figure 2: NMF example. Three rows correspond to original im- ages, k and k respec- tively.
matrix, but the sparse matrix operations are currently not well supported on mainstream deep learning platforms like PyTorch and Tensorflow. Since our method provides backward gradients for the NMF operation, the heatmap generated here can, in fact, be used to construct a loss function during training in order learn interpretable models.
Objective AUC
h
E
F
p
B
-1
E
Pij
IIINN
(f(x;)
[ c)
AU n
11,
-1f
Pij
(r) -e)
F-score
0
-1
(((()))
0
d
b
1
oxf(x,)))
Table 5: Table showing the general I. coefficients for each model. t; length based on problem
setting; 1: ft = di (z,, (C(r)) - fIz MC(=;)) + f(z, Bcc(z)) - f(x, Mc(())));; 8; one block for each E [1,..)) We do not include NMF in this table, as its formulation as a general I.P is more verbose including vectorization of matrices and kronecker product calculations.
14

---

Under review as a conference paper at ICLR 2021
A.4 VERIFICATION OF UNROLLING GRADIENT AND THE ONE PROVIDED BY À
We use Fscore formulation as an example. For input sample x, the neural network predicts a score [ar), and then the scores of batch of samples will be used in solving the linear programming form of Fscore and be used to construct the loss function. We compute the gradient from the final loss function back to the predicted scores from the neural network and compare two approaches: one is that we use z = Á-16 1 the solution (the one used in our experiment) where we can compute gradient by only one step, another one is that we directly use y: resulting from the Newton iterations as the solution and compute gradients by unrolling those iterations. We then compute the cosine value between these two gradient vectors. By experiments on 100 randomly sampled batches, the average cosine value is 0.9991, which means the two gradients are highly consistent.
B PROOFS AND DETAILS OF RESULTS IN SECTION
In this section, we will provide the missing proofs and additional calculations in Section 3.
B.1 PROOF OF LEMMA 1.
Lemma I is restated here for convenience.
Lemma 3. Assume that Rmxx a random matrix, and fix some y E RD. Then with probability one, g in equation is quadratic (given by H, Vg(y)) over a sufficiently small neighborhood of y.
Proof. Using the integral form of second order Taylor's expansion of o²()) = (max((,, y))² we can show that,
g + h) g(v) - h22 Volv) 1 h² A diag (d) Ah
(18)
2
where
d -- (/. (a (Ay-b)), b)), ds
2dt.
(19)
See Remark I in Golikov & Kaporin (2019) for details. Without loss of generality, we can assume b = 0 by simply translating the origin. Following the same remark, the diagonal matrix coincides with the step function based diagonal in f under the following condition on h:
(20) Since is fixed, assuming that the entries of A are chosen from a continuous distribution such that ef A uniformly distributed over the sphere, then (eî Ah)² follows Beta (L. ²=²) when h is drawn uniformly at random from the unit sphere, independent of A. This means that no matter what is, there exists sufficiently small h such that the left hand side of equation 20 false with probability one, and in that neighborhood diag(D) H.
An e (Ay) <0 = le Ahl S lej (Ay)l-
B.2 PROOF OF THEOREM 2
Theorem 4. Assume that the primal LP has a unique optimal solution, and that the level set {x : Ar < b, c² S a} is bounded for all a (for dual feasibility). Then short step (no line search) Newton's 's method converges globally at linear rate with local quadratic convergence.
Proof. First, since the objective function is piecewise quadratic since it is a sum of piecewise quadratic functions. In particular, it is self concordant since its third derivative zero everywhere. Now setting p < e, we see that an approximate solution of the problem with the modified Hessian is also an approximate solution to equation 6. Moreover, since the possible values of ft is finite, the local norm (also known as Newton's is decrement) Vy(y)" -VV(() is finite. Hence, we can choose / > so that there is a descent direction d, that is, there exist step size À > 0 such that A \Vg(x)² d 0. Finally, we use Theorem 4. in Nesterov (2013) to claim the desired result.
The assumptions in Theorem 4 are standard: 1. uniqueness easily be satisfied by randomly perturbing the cost vector; 2. in most of our formulations, we explicitly have bound constraints on the decision variables, hence level sets are bounded.
15

---

Under review as a conference paper at ICLR 2021
B.3 DIFFERENTIATING THE STEP FUNCTION IN REMARK 4
We will use slightly modified "suffix" notation as in Brookes (2005) in our calculations, That is, for a matrix A, Á is the same as vec(A), vectorization of A obtained by concatenating all the columns. The following three properties relating the Kronecker product, and differentials will be used often:
1. Fact For two vectors a, b, a ba
matrix, then 0B (0B/0A) ŌÁ where 2. Fact 2: If A is p matrix, and B is
0B/0A is the (mn) X (pq) Jacobian matrix of B with respect A. If A or column vector or scalar, then Thas no effect.
3. Fact 3: 3(AXB) = (BT & A) ax.
Using the above two facts, we can compute the gradients needed backpropagate through the unrolled iterations. We will show the computation for the gradient of Q u with respect to E [Rmxx for a fixed u Ra. We can apply chain rule to the following composition:
A
f1ofa
(A"MA+
f2
AT +
pl
to get, Joon = f,, o Jn. manipulation, we obtain,
Now using Fact 2 on 8 - - -xx and some algebraic
Joon - (n" (1"HAA+) ('TAAA* (( J
(21)
We will now compute Jn. Note that H is also function of A, using product rule, J/ as a sum of three derivatives with respect to each of A" H. The derivatives with respect to A and AT are fairly straightforward to compute, so will focus on computing the derivative with respect to H. To that end, we will use Fact 3, and show to compute the derivative of the step function by approximating it using the logistic function.
write
i a diag ((Ay - b).) 0A
diag (1 Ø (1 exp (N (- b)))) , <> 0.
(22)
Note that these derivatives are used in computing derivatives of upstream network, so using distribu- tional derivatives, and another application of chain rule to the left hand side of equation 22 results in the dirac delta function which is atomic, that is, has all its mass in measure zero set. Hence this calculation provides an mathematical justification that the set of nondifferentiable points has measure zero for our training purposes. It is easy formally verify this argument using differentiable tent functions as approximations to the heaviside step function.
16


Reviews:
[{"review_id": "sYMmZtYtOh", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "Summary\n-------\n\nThe paper makes the observation that various non-decomposable losses in machine learning can be rewritten as linear programs, whose constraints depends on the model output. This is the case for AUC, multi-class AUC, F-score, and to some extend NMF.\n\nThe authors review these losses, and recall how they may be rewritten as LPs. The LP formulation as known for AUC and NMF, but as far as the reviewer understand, they are new for multi-class AUC and F-score.\n\nThen, the authors propose to directly backpropagate through the LP resolution to minimize non-decomposable losses, applied on top of deep architectures. For this, they propose to solve an approximate solution to the LP problem (a quadratic penalization of the constraint violations) using a modified Newton method. They propose either to backpropagate by unrolling the Newton steps, or by using the computed minimizer directly.\n\nReview\n------\n\nThe endeavor of writing non-decomposable losses as LPs, to see these losses as pluggable LP-layer in deep architecture is interesting, albeit not original.\n\nUsing a penalized approximation of the LP to be able to solve them efficiently using a Newton method is also interesting.\n\nThe experiment section shows that it is indeed beneficial to directly optimize over a certain decomposable loss when we measure performance in term of this loss: in particular, it outperform using a simple logistic loss. This was completely expected, but it is good to verify it experimentally.\n\nOn the other hand, the manuscript suffer from many unclear parts, and from a theoretical analysis that is not polished enough. In particular:\n\n    - Phi is not introduced beforehand p. 4, and the F-score part is very hard to understand.\n\n    - the NMF section is very unclear, in particular as the authors use vague terms in their construction, such as \"zero padding ensures a sxs matrix\". I do not understand the role of tilde p in (6).\n\n    - Lemma 1 is not stated properly, as there is no f in equation (7). The authors state that \"each y has a neighborhood in which the Hessian is quadratic\", which does not mean anything. The proof sketch of Theorem 2 is very vague, in particular when the authors state that \"the possible choice of Hessian is finite\".\n\n    - I do not understand whether rho is chosen at every iteration, and what is its importance.\n\nI have trouble understanding why the authors went to such lengths in their\ntheoretical analysis. They modify a LP by making it a \"smooth almost everywhere\"\nproblem, which can then be solved using any methods, and backpropagated through\nusing either unrolling, or the computed minimizer (by virtue of Danskin\ntheorem), or the implicit function theorem. There is therefore not need to backpropagate throught tilde A^{-1} b.\n\nThe fact the the problem is only smooth almost everywhere may be a problem,\nwhich is not addressed by using a Newton method. It implies that the gradient\nbecomes a subgradient, and may hinder optimization performance. Remark 4\ndismisses this problem as unimportant, yet it is, as local convergence rates for\nnon-convex gradient descent requires smoothness.\n\nRelating to experiments:\n\n    - The reported performance does not show std errors across splits, which makes it impossible to compare in between similar methods (PPD-SG, PPD-AdaGrad and Ours). It appears that all three methods are within statistical variations.\n\n    - NMF is a long studied problem, with many powerful methods to handle large inputs. I do not understand the choice of using the input of a deep learning network for the experiment. As it it, the experiment proposed in this manuscript is not polished enough to be valuable.", "extra_scores": {}, "extra_text": {}}, {"review_id": "qzWgcjJWU35", "reviewer": "AnonReviewer4", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This paper addresses the classical topic of directly optimizing non-decomposable loss functions. Since these metrics can be computed via linear programs, it is sufficient to compute gradient through the LP solver. To that end, the authors propose to use a particular method for solving linear programs. In experiments, the resulting implementation outperforms the cross-entropy loss and mildly outperforms one recent baseline.\n\nThe topic of direct optimization of AUC, F1, and all the other related metrics is still of high practical interest and meaningful progress would be of high significance. However, the main reason why this topic is still open is that none (or a handful) of the dozens of top-tier publications addressing it, delivered a method that is reliable, stable, and scalable to practical setups (e.g. in computer vision). For this reason, I think, any new paper in this area needs to undergo a high level of scrutiny.\n\nI believe this paper has merit and constitutes a perfectly reasonable submission. At the same time, in its current form, it does not indicate sufficiently why it should have a different fate than all the other papers that also outperformed cross-entropy and subsequently disappeared.\n\nMy objections to the paper fall into three categories\n\n### Conceptual confusion \n\nIn this section, I want to clarify that gradients to LPs are a mostly solved problem and there is no dire need to resort to a new type of an LP solver -- as the authors do. I will also contradict the claim \"a simple closed-form gradient is not available for backpropagation\" made on page 4\n\nFor clarity let's focus on three scenarios (in the notation of the paper):\n\n1) Taking derivatives of the **optimal value of the objective function** w.r.t the LP parameters $c$, $b$, $A$\n\nIn all these cases, there is a simple closed-form gradient. E.g. if $x^*$ is optimal, the corresponding objective is $z = z(c) = c^T x^*(c)$ where the dependence of $x^*$ on $c$ is highlighted. This is a piecewise linear function ($x^*(c)$ is clearly piece-wise constant) with $dz/dc = x^*(c)$. This gradient $x^*(c)$ is simply the output of the LP and can be computed by any tailored method (not even necessarily by an LP e.g. by quicksort, Dijkstra, etc.) without requiring *any additional computation* on the backward pass.\n\nThe situation with $dz/db$ and $dz/dA$ is similar (via duality). This is all known for over 40 years ([M1, M2, M3, M4])\n\nIn the AUC LP formulation (eq 2), the resulting AUC score basically the negative optimal value of the objective, so **these classical results apply**.\n\n2) Taking derivatives of the **optimal solution**  $x^*$ w.r.t the LP parameters $c$, $b$, $A$\n\nIn case of $dx^*/db$ and $dx^*/dA$ there are again simple closed-form gradients. One can either refer to eq (6) in (Amos, 2017) \n where closed-form gradients of more general QPs computed, or to an argument the authors make themselves in Sec 3.2. $x^*$ is a solution to the linear system given by the set of active constraints, and gradients of matrix inversion are easy to compute. These are also the gradients given by CVXPy -- even though admittedly their LP solver is currently very slow.\n\nMost interesting is the case of $dx^*/dc$ as the true gradient is zero (dependence is piecewise constant) but this zero gradient is useless for optimization. Constructing meaningful gradient proxies is precisely the point addressed by (Berthet, 2020) but also earlier by [M5, M6, M7]. This is an ongoing research direction with competing methods.\n\n3) The LP solver is the last layer in the neural net and the ground truth solution $x^*_{\\textrm{true}}$ is available for supervision.\n\nIn this special case, which is actually the most common one in practice -- and also occurs in this paper, the situation is a lot easier. In fact, most works have focused on this scenario [M8, M9, (Song 2016)]. Also (Berthet, 2020) addresses it with the YF loss. The takeaway is that there are good loss functions for which a) the gradient can be computed only from the forward pass information b) robustness of the solutions can be enforced (e.g. by noise in (Berthet, 2020) or by a margin in M8).\n\nIn summary, there are suitable existing methods for efficient and often blackbox (not tied to a specific solver) gradient computation can be achieved. I find it imperative that the authors acknowledge this and compare to some of them.\n\n### Inaccurate claims and omissions in related work\n\nClaims about (Berthet, 2020)\n- \"approximate gradient\" - The gradient is approximate for good reason (see above)\n- \"necessary to solve ~n perturbed LPs\" - This is not true for the YF loss that can be applied at the end of the network (see above). There, in fact, only forward pass information is enough to obtain a gradient. Also, in the more general case (middle of the network) there are methods [M6] that need only one additional call to LP solver. \n- \"using off-the-shelf solvers ... can severely slow down\" - a) SOTA LP solvers are extremely fast, b) if there is a faster \"solver\" for a concrete LP, such as quicksort for the LP formulation of ranking, (Berthet, 2020) allows using it. In summary, (Berthet, 2020 and M5) allow the usage of the fastest possible algorithm and by definition cannot slow down runtime.\n\nMissing recent related work that also addresses nondecomposable metrics:\n[M10, M11, M12, M13]\n\n### Insufficient experimental evaluation\n\nAs far as I am concerned, the ideal experimental section for this paper would look like this:\n\na) more comparisons against recent literature on non-decomposable metrics (one recent baseline on one dataset is simply not enough)\nb) Identifying that the source of the performance is the selected method for solving LPs -- by comparing to alternatives described above. In that case,  additionally producing qualitative explanation on why this is happening.\nc) Find a fair way to report runtime and scaling with increasing instance sizes. Also comparing it to alternative approaches.\n\n# Summary\n\nMy intention is not to dismiss the method, it seems to work reasonably well and has a formal backing, after all. But if this paper is to be of long-term value to the community, I have to insist on a major revision for the reasons described above.\n\n### References\n\n[M1] Gal, 1975, Rim Multiparametric Linear Programming\n\n[M2] Freund, 1985, Postoptimal analysis of a linear program under simultaneous changes in matrix coefficients\n\n[M3] De Wolf, 2000, Generalized derivatives of the optimal value of a linear program with respect to matrix coefficients\n\n[M4] Gao, 2020, Differentiable Combinatorial Losses through Generalized Gradients of Linear Programs\n\n[M5] Vlastelica, 2020, Differentiation of Blackbox Combinatorial Solvers\n\n[M6] Ferber, 2020, MIPaaL: Mixed integer program as a layer\n\n[M7] Wilder, 2019,  Melding the data-decisions pipeline: Decision-focused learning for combinatorial optimization\n\n[M8] Tsochantaridis, 2005, Large Margin Methods for Structured and Interdependent Output Variables\n\n[M9] Elmachtoub, 2017, Smart \"predict then optimize\"\n\n[M10] Rolinek, 2020, Optimizing Rank-based Metrics with Blackbox Differentiation\n\n[M11] Brown, 2020, Smooth-AP: Smoothing the Path Towards Large-Scale Image Retrieval\n\n[M12] Fathony, 2020, AP-Perf: Incorporating Generic Performance Metrics in Differentiable Learning\n\n[M13] Khim, 2020, Multiclass Classification via Class-Weighted Nearest Neighbors\n\n", "extra_scores": {}, "extra_text": {}}, {"review_id": "Qtrs2DSgilf", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This paper shows how some nondecomposable functions can be\ninterpreted as solving combinatorial optimization problems\nthat can be relaxed to linear programs.\nAs this paper brings some new insights and directions here\nI recommend for a weak accept, although some of the experimental\nsettings and baselines feel incomplete (more details below).\n\n# Strengths\nOptimizing these performance metrics is useful for settings\nwhere these are more important than optimizing for the accuracy,\nand to the best of my knowledge these relaxations and ablations\nmeaningfully contribute to this direction.\n\nThe AUC experiments in Table 1 outperform [Liu 2019],\nwhich also optimizes for the AUC.\n\nThe derivative computation can be computationally expensive\nand ablating the impact of unrolling is interesting to\nbetter-understand tradeoffs.\n\n# Weaknesses\nThe biggest weakness I see is that the F-score experiment in\nTable 3 has no baseline that also directly optimizes the F-score,\nsuch as in [Fathony 2020] and other methods cited in the paper,\nand the non-negative matrix factorization experiments are\nlacking quantitative results.\n\nThis paper relies on relaxing the integer domain to be continuous\nand it's not clear how much this approximation impacts the derivatives.\n\nSome of the design choices seem arbitrary and unjustified, such as\nfocusing on the fast exterior penalty optimization at the start\nof Section 3 and then doing the forward pass with Newton's algorithm\non the unconstrained problem from [Mangasarian 2004] in Section 3.1\nto solve the LP.\n\n# Other comments and questions\nOne recent related work is omitted: [Fathony 2020].\n\nIs there any intuition behind why this approach sometimes outperforms\n[Liu 2019] in Table 1, but sometimes doesn't? My interpretation without\nthis is that both of these approaches are approximations and it's hard\nto know a-priori which will perform better.\n\nRemark 2 says that in previous work, the LP constraints do not\ndepend on the data, but the formulations of [Amos 2017, Agrawal 2019]\nallow constraints that depend on the data.\n\n# References\nFathony, R. & Kolter, Z.. (2020). AP-Perf: Incorporating Generic\nPerformance Metrics in Differentiable Learning. Proceedings of the\nTwenty Third International Conference on Artificial Intelligence\nand Statistics, in PMLR 108:4130-4140", "extra_scores": {}, "extra_text": {}}, {"review_id": "PJ4cSQMHOM", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "A. Summary:\nThis paper approximates several nondecomposable functions (AUC and F1-score) as linear programmings and uses them as loss functions for network training. In the linear programmings, the constraints are indeterministic at each mini-batch, the number of constraints increases quadratically to the number of training samples, so some previous works are inapplicable here. So does the primal-dual based forward pass and the corresponding implicit differentiation for a backward pass.\nInstead, the authors propose to select Newton's method for the forward pass and an adaptive online adjustment for the backward pass, where gradients are back-propagated through unrolled optimization before a proper accuracy is found, and then through complementary slackness afterward.\nFinally, the experiments emphasize the condition when the positive/negative examples are imbalanced, and demonstrates the superiority of the proposed methods to cross-entropy loss and one previous AUC loss.\n\nB. Strength:\n1. This paper gives a clear discussion about the difference between the proposed method and previous works, which motivates the author to explore an alternative solution to integrate linear programming as a differentiable component for training.\n2. It is a good choice to use Newton's method to tackle down the problem for training. For training, we care less about whether an optimization algorithm will achieve the global optimal solution, instead, we care more about whether the algorithm can achieve a local optimal rapidly, therefore the training will be feasible for larger-scale problems while the gradient can still be derived through complementary slackness.\n\nC. Weakness:\n1. Determining the accuracy \\epsilon seems to be critical for the proposed algorithm. It would be good if the author could have some more explanation about this. Besides, is that possible to make it a learnable variable through back-propagation at the first stage? Since the gradient is backpropagated via unrolling at first, it should be doable and save some manual effort to decide how to increase it.\n2.  In the experiments, the authors constructed an imbalanced dataset from the existing dataset, In addition, some other recent works, such as Lin, TY, Goyal, P, Girshick, R, He, K, Dollar, P, \"Focal Loss for Dense Object Detection\", ICCV'2017, also focus on deal with the class imbalance by other tricks. So it would be better to compare with these methods as well but not only vanilla cross-entropy. \n3.  The application of nonnegative matrix factorization should be further improved. As indicated in the paper, it only works when the number of channels is limited. The author may remove this application and focus on the AUC and F1-score. The contribution will be sufficient enough if these two widely used metrics can be discussed more solidly.\n4.  What if we fine-tune a pre-trained model using the proposed losses and the dataset is not manually selected to make it imbalance? Will the performance achieve consistent improvements?\n\nD. Justification of the score:\nIn general, this paper focus on an important problem, and the algorithm is discussed comprehensively. My major concern is about the experiments, where the experimenting settings are biased to the proposed algorithm. I will raise my score if these concerns can be addressed during rebuttal.\n", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper shows that various discrete loss functions can be formulated as an LP. It proposes to relax the constraint Ax = b, x >= 0 using a soft constraint and following Mangasarian, proposes to solve the relaxed problem using Newton's method. Backpropagation through these iterations is further proposed. The main motivation is that this results in a GPU-friendly implementation.

I think the proposed approach is novel. However, as pointed out by reviewers, the current writing lacks clarity and the experiments are quite weak. There is now a wealth of methods for differentiating through an LP using implicit differentiation, smoothing (which the present paper is a form of, see below) and perturbations. It is important to compare to these methods. The paper also ignores a large literature on convex surrogates for ranking metrics.

I recommend the authors to strengthen the writing and experiments, and to resubmit to a top-conference.

Additional comments by the AC
-------------------------------------------

As the sentences "Hence, solving such LPs using off-the-shelf solvers may slow down the training process" or "Often, this would involve running the solver on the CPU, which introduces overhead" indicate, the authors seem to imply that LPs need to be solved in canonical LP form, min_x <c,x> s.t. Ax <= b, x >=0, using an off-the-shelf LP solver. This is not how many LPs are solved in practice. For every loss, there will always be an ad-hoc solver for the corresponding LP. For instance, the Hungarian algorithm for the Birkhoff polytope.

The paper is missing an important reference: SparseMAP (https://arxiv.org/abs/1802.04223). In this paper, the authors add regularization to the primal LP and use Frank-Wolfe or active set methods to solve the problem.

Equation (6) corresponds to relaxing the hard constraint Ax=b, x>=0 with a soft one. This approach is in a sense opposite to SparseMAP. Indeed, relaxing the constraints in the primal is equivalent to adding regularization in the dual LP (see, e.g., https://papers.nips.cc/paper/2012/hash/bad5f33780c42f2588878a9d07405083-Abstract.html).
Speaking of (6), the authors should clarify that it's a convex objective.

In section 2, the authors review a number of losses which can be written as an LP. It would be better to explicitly state what are A, b and c for each loss (or g, h, E, F, p, B, G, q).

The matrix A could potentially be huge, depending on the LP. Do you need to materialize it in memory in practice? This would limit the approach to relatively small LPs.
Decision: Reject