Differentiable Optimization of Generalized Nondecomposable Functions using Linear Programs
Year: 2021

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 performance measures such as AUC, multi-class AUC, F -measure and others. A 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 characterize the constraints or the feasible set. The constraint matrix is not only large but the constraints are also modified at each iteration. We show how adopting a set of ingenious ideas proposed by Mangasarian for 1-norm SVMs – which advocates for solving LPs with a generalized Newton method – provides a simple and effective solution that can be run on the GPU. In particular, this strategy needs little unrolling, which makes it more efficient during the backward pass. Further, even when the constraint matrix is too large to fit on the GPU memory (say large minibatch settings), we show that running the Newton method in a lower dimensional space yields accurate gradients for training, by utilizing a statistical concept called sufficient dimension reduction. While a number of specialized algorithms have been proposed for the models that we describe 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.

heading: 1 Introduction
text: Commonly used losses such as cross-entropy used in deep neural network (DNN) models can be expressed as a 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 many use cases. But in applications ranging from ranking/retrieval 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 a fixed recall (P@R) and others. Optimizing such measures in a scalable manner can pose challenges even in the shallow setting.
For AUC maximization, we now know that convex surrogate losses can be used in a linear model Liu et al. [2018], Natole et al. [2018] in the so-called ERM framework. These ideas have been incorporated within deep neural network models and solved using SGD type schemes in Liu et al. [2019]. Such results on stochastic and online data models have also been explored in Ataman et al. [2006], Cortes and Mohri [2004], Gao et al. [2013]. There are also available strategies for measures other
35th Conference on Neural Information Processing Systems (NeurIPS 2021), virtual.
than the AUC: Nan et al. [2012], Dembczynski et al. [2011] give exact algorithms for optimizing F -score and Eban et al. [2017], Ravi et al. [2020] proposes scalable methods for non-decomposable objectives which utilizes Lagrange multipliers to construct the proxy objectives. The authors in Mohapatra et al. [2018] discuss using a function that upper bounds (structured) hinge-loss to optimize average precision. Recently, Fathony and Kolter [2020] presented an adversarial prediction formulation for such nondecomposable measures, and showed that it is indeed possible to incorporate such measures within differentiable pipelines.
Our work utilizes the simple observation that a number of these non-decomposable objectives can be expressed in the form of an integer program that can be relaxed to a linear program (LP). 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-built support for specialized GPU hardware in modern libraries such as PyTorch, is desirable. First, as long as a suitable LP formulation for an objective is 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-decomposable objective based on simple calculations (e.g., number of constraints, primal-dual gap). The current tools do not entirely address all these requirements.
A characteristic of the LPs that arise from the nondecomposable losses mentioned above is that the constraints (including the mini-batch of samples themselves) are modified at each iteration – as a function of the updates to the representations of the data in the upstream layers. In Section 3, we provide LP formulations of widely used nondecomposable terms, which fall squarely within the capabilities of our solver. In Section 4, we show that the modified Newton’s algorithm in Mangasarian [2004] can be used for deep neural network (DNN) training in an end-to-end manner without requiring external solvers, where support for GPUs currently remains limited. Specifically, by exploiting self-concordance of the objective, we show that the algorithm can converge globally without line search. We then analyze the gradient properties of our approach, and some modifications to improve stability during backpropagation. Our experiments in Section 5 show that this scheme based on Mangasarian’s parametric exterior penalty formulation of the primal LP is a computationally effective and scalable strategy to solve LPs with a large number of constraints. On the practical side, we provide two ways to deal with the scaling issue when the constraint matrix is large. On the one hand, we show that sufficient dimension reduction can be used in our solver to solve the problem in a lower dimension space. On the other hand, when the matrix is sparse, with our new sparse implementation, we show that we can train Resnet-18 with 10× mini-batch size (memory savings) on a 2080TI Nvidia GPU. Our code is available, see Fig 1 for an overview of our solver for differentiating many popular nondecomposable objectives.

heading: 2 Related Works
text: Differentiable optimization of nondecomposable objectives. One closely related result to ours is Liu et al. [2019] which explored optimizing a surrogate loss of AUC using SGD type methods that made the stochastic AUC optimization more practical for deep neural networks. Compared with Liu et al. [2019], our work enables optimizing many more types of nondecomposable objectives for deep neural networks, and relies on the LP reformulation instead of a surrogate loss. In Fathony and Kolter [2020], (the expected values of) such metrics are approximated by LPs defined by predictor and adversary marginal distributions (represented as square matrices). The ADMM solver in Fathony and Kolter [2020], as we will see, inherits some drawbacks of ADMMs including higher memory footprint and slower convergence. Another body of work is the one involving the use of blackbox solvers Pogančić et al. [2019], Berthet et al. [2020] which can be adapted to solve the nondecomposable objectives (e.g., reformulating as LPs). For example, one idea is to use appropriate (infimal or gaussian) convolution to construct approximations of the gradients through the LP. Song et al. [2016] proposes a using a task loss to perturb the loss and is similar in spirit to blackbox solvers – they both solve a sequence of perturbed argmax/argmin problems first, and then use carefully designed procedures for the backward gradient that often require at least one more solve of the optimization problem. We will discuss these properties later in the paper, and in our experiments.
Differentiable LP solver. In principle, of course, backpropagating through a convex optimization model (and in particular, LPs) is quite possible. 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 a powerful external solver. Often, this would involve the overhead of running
the solver on the CPU. Solvers within CVXPY, are effective but due to their general-purpose nature, rely on interior point methods. OptNet Amos and Kolter [2017] 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 to LPs due to the polyhedral geometry in LPs. The ideas in Meng et al. [2020] are only applicable when the number of constraints is approximately equal to the number of variables – an invalid assumption for the models that we study here.

heading: 3 Nondecomposable Functions and corresponding LP models
text: We first present a standard LP form and then reparameterize several generalized nondecomposable objectives in this way, summarized in a table in the appendix. We start with the binary AUC, extend it to multi-class AUC, and then later, discuss a ratio objective, F -score. The appendix also includes a discussion of other objectives that are amenable to our method.

heading: 3.1 Notations and Generalized LP formulation
text: Notations. We use the following notations: (i) n: number of samples used in training. (ii) X ∈ Rn×dimx : the explanatory features fed to a classifier (e.g., parameterized by w); (iii) f(xi) (or f(i)): a score function for the classifier where xi ∈ X; (iv) Y ∈ {0, 1}: target label and Ŷ ∈ {0, 1}: predicted label for binary classification, both in Rn; (v) A⊗B: Kronecker product of matrices A and B. (vi) Ir: Identity matrix of size r and 1 is the indicator function. (vii) Bk, p (and B|,k′ ) gives the k-th row (and k′-th) column of B.
LP formulation. We consider a general linear program (LP) that contains nonnegative variables as well as inequality and equality constraints. The form of the LP is given as
max ũ,ṽ
gT ũ+ hT ṽ subject to Eũ+ F ṽ ≤ p, Bũ+Gṽ = q ũ, ṽ ≥ 0
We can write it more succinctly as Variable z = [ũ ṽ] ; Coefficient c = [−g − h] ; Constraints A = [
E B −B F G −G
]T ; Constants b = [p q − q]T
The corresponding primal LP can be written as, minz cT z subject to Az ≤ b, z ≥ 0.

heading: 3.2 Maximizing AUC
text: The Area under the ROC Curve (AUC) calculates the probability that a classifier f(·) will rank a randomly chosen positive sample higher than a randomly chosen negative sample. AUC varies
between 0 and 1, where 1 represents all positives being ranked above the negatives. AUC may be estimated using the Wilcoxon-Mann-Whitney (WMW) Statistic Hanley and McNeil [1982], as
Definition 3.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+|+ |X−| = n where | · | is the cardinality of the set. Then, AUC is given as (∑|X+| i=1 ∑|X−| i=1 1f(xi)>f(xj) )/( |X+| |X−| ) for xi : i ∈ {1, · · · , n}.
Here, we follow Ataman et al. [2006] to calculate the AUC based on the WMW statistic as follows.
min zij |X+|∑ i=1 |X−|∑ j=1 zij s.t. f(xi)− f(xj) ≥ ϵ− zij where xi ∈ X+, xj ∈ X−; zij ≥ 0, (1)
where ϵ is a given constant. Problem (1) computes AUC indirectly by minimizing the number of pairs (one each from the positive and negative classes) where the positive sample is not ranked higher than the negative sample: so, the number of zero entries in z equals the number of pairs where this condition is not true. For a given z, we have that,
AUC = (n− ∥z∥0) / (|X+| |X−|) = ( n− ∑ i ∑ j ϵ−1relu(0,−zij + ϵ) )/ (|X+| |X−|) .
If zij is 0, then ϵ−1relu(0,−zij + ϵ) equals 1. Otherwise zij > 0, it follows from the first constraint in (1), that zij ≥ ϵ, so ϵ−1relu(0,−zij + ϵ) equals 0. Observe that in (1), the number of constraints is |X||X−|, which is O ( n2 ) .

heading: 3.3 Maximizing Multi-class AUC
text: An extension of AUC to the multi-class setting, AUCµ, is defined in Kleiman and Page [2019]. The AUCµ objective optimizes an indicator matrix calculated on the orientation function, Oi,j defined as,
Definition 3.2 (Orientation Function; Kleiman and Page [2019]). Assume we have K classes {y1, · · · , yK}. Let f(xi, p) ∈ RK (extension of f(.) to the multi-class setting) indicate the model’s prediction on xi for each of the K classes (class-specific probability). Let x∗i provide the index of xi’s true class label. Let P ∈ Rk×k be a partition matrix where Pk,k′ is the cost of classifying a sample as class k when it should be k′. Define,
vkk′ = Pk, p −Pk′, p and ṽ = vx∗i x∗j ∈ R K ; Oi,j = (ṽx∗i − ṽx∗j )(⟨ṽ, f(xi, p)⟩ − ⟨ṽ, f(xj , p)⟩).
The LP formulation of AUCµ is fully characterized by P. We can set P(k, k) = 0 ∀k and 1 for all other entries. We can also define a P with arbitrary entries or formulate AUC in a one-vs-all setting (denoted as AUCova). Here, for presentation, use the simplified 0/1 partition matrix P. Let f̃(i, j) = f(xi, x ∗ i )− f(xj , x∗i ) + f(xj , x∗j )− f(xi, x∗i ). Then the LP formulation is given by,
AUCbinµ : min zij n∑ i=1 n∑ j=1:x∗i <x ∗ j zij s.t. d̃ij f̃(i, j) ≥ ϵ− zij , ∀i, j : x∗i < x∗j zij ≥ 0 (2)
where d̃ij = ṽx∗i − ṽx∗j . Problem 2 can be seen as an extension of our binary AUC model, where zij is the ranking between a pair of points defined for multiple classes.

heading: 3.4 Maximizing F -score
text: The F -score (or F -measure) is a representative of objectives expressed as ratios of some combination of True positives (TP), False positives (FP), True negatives (TN) and False negatives (FN). We use the result in Dembczynski et al. [2011] to express the F -score in the ratio form and further relax it into the LP form. The detailed LP formulation of F -score and other ratio functions is in the appendix.

heading: 4 Backpropagation via Fast Exterior Penalty Optimization
text: Unlike traditional feedforward networks, where the output of each layer is a relatively simple (albeit non-linear) function of the previous layer, a LP layer must solve a constrained problem, therefore implementing scalable/efficient backpropagation schemes that minimizes overhead requires more care and 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. [2017] or use projections Zeng et al. [2019]. Recently 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 1. Some previous works Zeng et al. [2019] have considered LPs where the constraints are deterministic (for a fixed input dimension), i.e., do not depend on the data X , which is different from the LPs in Sections 3.2–3.4.
Note that perturbation techniques in Berthet et al. [2020] are applicable to our LPs as well. The Fenchel Young losses in Berthet et al. [2020] is attractive because there is no need to compute the Jacobian. Implementation-wise, one could think of the backward pass as a function given the input and output of the forward pass. But the gradient expressions of the losses involves an expectation and hence may require multiple calls to a LP solver in order to approximate the expectation. Parallelization and warm starts were shown to alleviate this dependency to some extent by sampling in parallel.
Rationale for our approach. Consider a LP with a large m number of constraints in fixed dimensions n (n ≪ m). This assumption holds in all formulations in Section 3. This is because we assume that the architecture is fixed whereas minibatch size depends on the complexity of the task (stable gradient or when noise in gradient is high). Hence, solving such LPs using off-the-shelf solvers as in Berthet et al. [2020], in the general case, could slow down training. The strategy in Agrawal et al. [2019] does offer benefits over Amos and Kolter [2017] for sparse QPs. Our strategy is to run Mangasarian’s Newton’s method on an exterior penalty function of the LP. There are two advantages: (i) during forward pass, quadratic local convergence of Newton’s method indicates that unrolling the method may be reasonable; 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 results and some modifications to deal with discontinuous Hessian (and its inverse) that is required for Newton’s method. A similar approach is adopted in Amos and Kolter [2017] where Primal-Dual Interior Point methods with implicit differentiation is used for differentiation purposes. But the exterior penalty in (3) satisfies a nice property: primal and dual solutions are related by a closed form expression which can be exploited for efficient backpropagation.

heading: 4.1 Forward Pass using Newton’s Algorithm on a Sufficiently Reduced Space
text: Fast automatic (forward or reverse mode) differentiation requires performing the forward pass efficiently. In our setting, we 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 a primal LP, for a fixed accuracy ε > 0, Mangasarian [2004] solves an unconstrained problem,
min y
g(y) := 1
2 ∥σ(Ay − b)∥2 + εcT y, (3)
where σ(·) = max(·, 0) represents the elementwise relu function. A modified Newton’s method can be used to solve (3) that performs the following iterations:
y = y + λd (4)
where d = H̃(y)−1∇g(y) := (∇2g(y) + ρI)−1∇g(y) (5)
In large scale settings of A, b, such Newton methods are known to perform empirically better than gradient descent Mangasarian [2006], Keerthi et al. [2007]. We will evaluate 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 (even for convex problems), 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 above 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 a rate of convergence. In DNN training, such line search based convergence results can be very expensive. The key difficulty is handling the discontinuity in the Hessian. As a remedy, we use self concordance of (3) to guarantee global convergence of (4) iterations for the exterior penalty formulation in (3). To do so, we first show a result (proof in appendix) that characterizes the discrepancy between the actual Hessian of (3) and the modified one in (4) when the (feature) matrix A is randomly distributed. Lemma 1. Assume that A is a random matrix, and fix arbitrary y ∈ Rn. Then with probability one, g in (3) can be approximated by a quadratic function (given by H̃,∇g(y)) over a sufficiently small neighborhood of y.
Intuitively, Lemma 1 states that with probability one, each y has a 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. Assume that the primal LP has a unique optimal solution, and that the level set {z : Az ≤ b, cT z ≤ α} is bounded for all α (for dual feasibility). Then short step (no line search) Newton’s method converges globally at a linear rate with local quadratic convergence.
Proof. First, note that the objective function is piecewise quadratic since it is a sum of piecewise quadratic functions defined by coordinatewise relu function σ. In particular, g is self concordant since its third derivative is zero almost everywhere. Now setting ρ < ε, we see that an approximate solution of the problem with the modified Hessian is also an approximate solution to equation 3. Moreover, since the possible values of H̃ is finite, the local norm (also known as Newton’s decrement) ∇g(y)T H̃(y)−1∇g(y) is finite. Hence, we can choose ρ > 0 so that there is a descent direction d, that is, there exists a step size λ > 0 such that λ∇g(x)T d < 0. Finally, we use Theorem 4.1.12 in Nesterov [2013] to claim the desired result.
The assumptions in Theorem 2 are standard: 1. uniqueness can 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. Remark 2. Convergence in Thm. 2 is guaranteed under standard constraint qualification assumptions. Linear Independent Constraint Qualification (LICQ) is satisfied for AUC, and Multi-class AUC formulations in §3. 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).
Our analysis of the Newton’s method for LP so far indicate that we may be able to use a constant step size λ (avoid linesearch) to obtain fast convergence provided that we are able to choose ρ sufficiently small. For our purposes, we assume λ to be a hyperparameter and can be tuned by cross validation.
Sketching d using Sufficient Dimension Reduction (SDR). For training the backbone network, we have to compute d in (4) using A (and b) which could be very large within each training iteration. But each minibatch corresponds to an LP instantiation in (3), so it may suffice to compute an approximate d in each Newton’s iteration. To approximate d, we use sufficient dimension reduction in which a lower dimensional LP obtained by projecting A is solved Zhang et al. [2020], Kim et al. [2020]. That is, instead of the inverse of H̃ , we will compute a lower dimensional sketch of H̃ by using SA (and Sb) instead of A (and b) for some sampling matrix S and solve problem (3). During each iteration, the metrics are computed on the current minibatch as is done in Fathony and Kolter [2020]. Thus (without SDR) the memory cost is directly proportional to minibatch size viz., number of samples and feature dimensions. The advantage of using SDR (as opposed to naive sketching) is that the lower dimensional space can also be chosen using data driven cross validation or other techniques
Kim et al. [2016], making it ideal for training purposes. In essence, the size of minibatch does not matter – as long as the minibatch can be sufficiently reduced, our solver is directly applicable, especially for low resource, memory constrained settings. For a fixed (low) dimension parameter, our sketched d is,
d = H−1g where g = ∂
∂y (∥relu(Ay − b)∥22) ≈ (SA)TS · relu(sign(Ay − b)), (6)
H ≈ (diag(S · relu(sign(Ay − b))⊙ S · relu(sign(Ay − b))) · SA)T · SA. (7)
It is easy to see that the approximate update can be seen as equivalent to Iterative Hessian Sketch algorithm which has geometric convergence rate Pilanci and Wainwright [2016].

heading: 4.2 Backward Pass using Optimal Dual Variables Aided by Unrolling
text: The advantage of optimizing the exterior penalty in (3) is that given an iterate yt, accuracy ε, we can get the optimal dual solution vt by simple thresholding, i.e., vt = 1/ε(ATσ(Ay − b)). By complementarity slackness, nonzero coordinates of vt specify the set of active constraints in Az ≤ b. So, given an approximate yt such that ∇g(yt)H̃(yt)−1∇g(yt) ≤ ε, to get the primal solution z∗, we solve the “active” linear system given by Ãb̃, where Ã denotes the active rows of A and the corresponding subvector b̃. So, backpropagation through the layer reduces to computing derivatives of Ã−1b̃ (simple via automatic differentiation). The appendix gives a systematic way of choosing ε.
How to backpropagate through unrolled iterations? We assume that the chain rule is applicable up to this LP layer and is ∂L∂z ∂z ∂A (for one of the parameters A), and note that it is possible to find ∂L ∂z (either directly or using a chain rule). Therefore we focus our attention on ∂z ∂A , which involves the LP layer. Indeed, unrolling each iteration in (4) is equivalent to a “sublayer”. So, in 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 calculate ∂d∂A where d = Q −1q, Q = H̃ and w = ∇g(y). We can use the product rule to arrive at: ∂d = − ( wTQ−1 ⊗Q−1 ) ∂Q + Q−1∂w. 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 ∂Q∂var or ∂w∂var where var ∈ {c, A, b} which can also be done by another application of chain rule. Please see appendix for empirical verification of unrolled gradient and the one provided by Ã−1b̃.
Before proceeding, we should note an issue that comes up when differentiating each step of the unrolled algorithm because the Hessian is piecewise linear (constant) as a function of the input to that particular layer. Here, some possible numerical approximations are needed, as we describe below.
Remark 3. Note that the diagonal matrix term in ∂Q∂A is nondifferentiable due to the presence of 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, z ∈ S, if a ball (of radius r > 0,) Br(z) ⊆ S, then the Lebesgue measure of the set of nondifferentiable points on S is zero. Please see appendix for a formal justification where we show this by approximating the step function using a sequence of logistic functions with increasing slope parameter at the origin.
Therefore, in this setting, Remark 3 provides a way to compute an approximate sub-gradient when using Newtons method based LP layers. The function is a piece-wise quadratic function and differentiable everywhere, and the inverse of the Hessian acts as a preconditioner.
Summary. Our forward pass involved three steps: (1) finite steps of Newton’s method using which we (2) computed the dual variable by a thresholding operation, and (3) finally, to get the primal solution, these dual variables are first used to identify the active constraints followed by solving a linear system. To backpropagate through these three steps, we must differentiate through each layer of our procedure including Ã−1b̃: this is independent of whether we use unrolling or Danskin’s theorem. For instance, using Danskin’s theorem here would involve differentiating through the fixed point of the Newton’s iterations similar to (regularized) gradient descent iterations in iMAML Rajeswaran et al. [2019].

heading: 5 Experiments
text: 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. Nonnegative matrix factorization is discussed in the appendix.
Optimizing Binary AUC
We follow the 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 algorithms (PPD-SG and PPD-AdaGrad) from Liu et al. [2019].
Datasets: Cat&Dog, CIFAR10, CIFAR100, and STL10 (See appendix for dataset details). We follow Liu et al. [2019] to use 19k/1k, 45k/5k, 45k/5k, 4k/1k training/validation split on Cat&Dog, CIFAR10, CIFAR100, STL10 respectively.
Experimental setting. We follow Liu et al. [2019] to construct an imbalanced binary classification task by using half of the 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 (ratio of majority class to the minority class) in Table 1. Con-
structing imbalanced datasets by dropping samples was also used in Cui et al. [2019] to construct long tailed CIFAR-LT. We use the same random seed, learning rate and total number of iterations in all of our experiments including multi class AUC and F -score experiments. See appendix for model architecture, learning rate, etc.
Results. The results are shown in Table 1. We can see that our method slightly outperforms Liu et al. [2019] and outperforms cross-entropy loss by a large margin, especially on imbalanced datasets, where the AUC objective shows superiority over cross-entropy loss.
Per-iteration complexity. The wall clock run time per iteration of cross-entropy, our solver, PPDSG and PPD-AdaGrad are 0.018, 0.069, 0.122 and 0.130 respectively (in terms of seconds). Our method and PPD cost more than directly optimizing cross-entropy because of the additional cost of solving AUC and our method is more time efficient compared with PPD.
Influence of ϵ in our formulation 1. We test the influence of ϵ using Cat&Dog as an example. Results (see appendix for the table) show that ϵ = 0.1 gets slightly worse performance than ϵ = 0.01, while ϵ = 0.001 performs much worse. To choose ϵ, we follow the approach in Mangasarian [2004]. If for two successive values of ϵ1 > ϵ2, the value of the ϵ perturbed quadratic function is the same, then it is the least 2-norm solution of the dual. Therefore, we simply choose an ϵ that satisfies this property, which is chosen to be 0.01 in our experiments.
Optimizing Multiclass AUC
We also evaluate our method for multiclass AUC. Similar to binary AUC, we construct imbalanced multiclass datasets by dividing datasets into 3 classes and drop samples from 2 of them and report one-versus-all AUC (as AUCova) and AUCbinµ score
. For STL10, we group class 0 − 2, 3 − 5, 6 − 9 into the three classes, and drop samples from the
first two classes. For CIFR100, we group class 0− 32, 33− 65, 66− 99 into three classes, and also drop samples from the first two classes.
Results. Results are in Table 2. In addition to one-versus-all AUC metric, we also report the performance in terms of AUCµ Kleiman and Page [2019] which is specifically designed for measuring multiclass AUC and preserves nice properties of binary AUC such as being insensitive to class skew. Our method outperforms cross-entropy loss on all four datasets and under all different skewed ratios.
Optimizing F1-score
We show that by directly optimizing F1-score, we can achieve a better performance than when using cross entropy loss. In addition to cross entropy loss, we perform evaluations with two other methods that can also directly optimize the F1-score. First, we replace our solver with CVXPY-SCS Agrawal et al. [2019], which refers to the SCS solver used within CVXPY package; second, we perform comparisons with AP-Perf Fathony and Kolter [2020] which offers differentiable optimization of F1-score using an adversarial prediction framework. The datasets and our setup to group them into two classes remain the same as in binary AUC section. The results show that our method yields consistent improvement over cross-entropy loss in terms of F1-score. The accuracy of cross-entropy loss on Cat&Dog, CIFAR10, CIFAR100 and STL10 is 76.0%, 70.3%, 60.4%, 71.8%, respectively, while our method gives 77.8%, 72.6%, 63.4%, 72.7%. When we optimize F1-score directly, there exists a local optimal point where assigning all examples to the positive class leads to F1-score of 66.7%. Both CVXPY-SCS and AP-Perf fall into this local optimal on four datasets (except that CVXPY-SCS gets out of the local optimal point on Cat&Dog and gives 70.1% accuracy). Note that although our solver, CVXPY-SCS and AP-Perf seek to solve the same objective, the backward gradient will be different due to different approximations used by the methods.
Comparing with blackbox solvers. Recently, several results Pogančić et al. [2019], Berthet et al. [2020] use blackbox solvers for solving combinatorial problems, which can also be used as a LP solver for problems where we need dz∗/dc (c is the cost vector in standard form LP). These methods can utilize existing LP solvers and compute backward gradient by calling the LP solvers on perturbed LP problems. To our knowledge, the current state-of-the-art LP solvers are often CPU based, which means that one needs to transfer the LP parameters from GPU to CPU, then solve the LP on the CPU before getting the solution back. These CPU-based solvers will usually be slower than our GPU based solver, especially when the batch size is large. To evaluate this property, we randomly
construct LP problems using AUC as an example and make the cost vector the parameter to optimize. In Fig. 3b, we compare the runtime of our solver and the blackbox solver in Pogančić et al. [2019] (with the LP solver from scipy package). We can see that our GPU based solver is not sensitive to the increase in the batch size, while the CPU based solver is, due to the lack of support for mini-batch operations. Our solver has a more clear advantage in backward pass because the matrix needed for computing gradient is already computed during forward pass thus the backward pass is nearly free, while the blackbox solver Pogančić et al. [2019] needs at least the same time as forward pass to compute backward gradient. We also compare our solver with Berthet et al. [2020] in Table 3.
Dealing with large scale LPs. In practice we often encounter large scale LPs whose constraint matrix may be too large to fit into the GPU memory. We first show that we can utilize sufficient dimension reduction in Fig. 3c, which demonstrates a reasonable tradeoff between reduced size and accuracy. Then we show that when the constraint matrix is sparse, we can readily utilize sparse operations to save memory. We have incorporated functionality to use sparsity within our solver flaport [2020], rusty1s el al. [2020]. Consider AUC maximization in which the constraint matrix is extremely sparse. With our sparse implementation, we are able to run SGD on Resnet-18 with minibatch size of up to 200 on a 2080Ti GPU with just 1GB memory whereas the same problem can take upto 11GB using dense operations (≈ 10× memory savings) with some overhead. Limitations. Scaling is still a key limitation for differentiable LP solvers. We present two ways to mitigate these problems but there is an associated cost. Our proposed SDR sacrifices some accuracy. Partly due to support within available libraries sparse operations can sometimes be slow when the size of the matrices is large which should improve with better sparsity support in existing libraries.

heading: 6 Conclusions
text: We demonstrated that various non-decomposable objectives can be optimized within deep neural networks in a differentiable way under the same general framework of LPs using a modified Newton’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-dependent. Our code is available at https://github.com/zihangm/nondecomposable-as-lp.
Potential negative societal impacts. Optimizing non-decomposable measures such as AUC is a fundamental task in statistics/machine learning and not tied to specific applications, so we do not see a negative societal impact of the technical development described here. It is possible, however, that modules such as these can be used within general applications that are detrimental to society.

heading: 7 Acknowledgments
text: This work was supported by NIH RF1AG059312 and RF1AG062336, and NSF CCF #1918211. We thank Sameer Agarwal for Twitter discussions regarding implicit differentiation in deep networks, and later pointing to Olvi’s paper over email correspondence. Sathya Ravi was supported by UIC start-up funds. Yichao Wu was supported by NSF grant DMS-1821171.

Review summaries:
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.

This 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.

This work proposes to relax some commonly used discrete metrics into a Linear Program (LP). By using the Newton LP method, the algorithm is able to backpropagate through the LP and thus to train a neural network end-to-end while directly optimizing for the metric. In particular, this submission describes how the AUC, multi-class AUC and F-score can be casted as LPs, and how the resulting problems can be solved with a Newton LP method that can be differentiated through. Experiments are presented on Cat&Dog, CIFAR-10, CIFAR-100 and STL.

This paper present a novel approach to perform direct optimization of non-decomposable objective functions (such as AUC or F-measure) in the context of Deep Learning. They propose a general LP-based framework describing these objective functions and provide a general algorithm to solve this problem (and get the corresponding gradients). Finally, this new approach is compared to other methods trying to optimize these objectives directly and surrogate loss functions that are commonly used in practice.