The rapid advancement of autonomous vehicle (AV) technology is reshaping ride-hailing markets, where transportation network companies (TNCs) increasingly operate mixed AV and human-driven vehicle (HV) fleets. This transition demands an enhanced planning framework that is flexible, accommodating varying AV penetration rates, expressive, capturing operational heterogeneity of mixed autonomy, and unifying. Termed MAGE (Mixed-Autonomy General Equilibrium with Customer Patience), the resulting framework is structured around four coupled modules and distinguishes three behavioral classes: AVs, centrally coordinated by TNCs for profit-driven dispatch and fleet composition; HVs, which serve TNC demand but exhibit decentralized behavior due to human spontaneity and individual preferences; and solo drivers (SVs), who minimize individual travel costs. The behavioral asymmetry between AVs and HVs is explicitly preserved across both pickup and service stages--the two operational phases unique to TNC fleets, where deviations from the Wardrop principle arise through distinct mechanisms, while SVs adhere to this principle. Customer patience endogenously links supply-side heterogeneity to demands: travelers trade waiting for ride-hailing service within patience against solo driving, thus capping wait times and influencing fleet deployment within a flow-dependent congestion network. The overall system is mathematically formulated as a Nonlinear Complementarity Problem, which we show is equivalent to a Variational Inequality. We establish equilibrium existence through a new proof that relaxes a restrictive fleet-size assumption in prior work. Numerical results reveal that AV routing lever amplifies competitive advantages among TNCs but exhibits diminishing marginal returns, eventually saturating, while traveler patience stabilizes the market by implicitly disciplining excessive AV routing control.

We study online optimization for a broad class of structured non-convex non-smooth problems where each loss is a composition of a difference-of-convex function with a smooth mapping, and the feasible region is defined by constraint functions of the same kind. We propose a time-smoothed proximal linear algorithm and a local-regret measure based on a proximal residual mapping. We show that this residual is a proper stationarity measure for the original problem: its fixed-point condition implies first-order stationarity. Our analysis relies on a tangent-cone characterization for a feasible region described by composite difference-of-convex constraints, which is of independent interest and allows each update to be computed via a convex optimization oracle, despite the non-convexity of the problem. We establish a local-regret bound and a bound on the total number of inner convex subproblems. We also derive an error bound connecting the proximal residual to the distance to stationarity, providing a quantitative certificate of approximate stationarity.

Offline policy learning aims to use historical data to learn an optimal personalized decision rule. In the standard estimate-then-optimize framework, reweighting-based methods (e.g., inverse propensity weighting or doubly robust estimators) are widely used to produce unbiased estimates of policy values. However, when the propensity scores of some treatments are small, these reweighting-based methods suffer from high variance in policy value estimation, which may mislead the downstream policy optimization and yield a learned policy with inferior value. In this paper, we systematically develop an offline policy learning algorithm based on a weight-clipping estimator that truncates small propensity scores via a clipping threshold chosen to minimize the mean squared error (MSE) in policy value estimation. Focusing on linear policies, we address the bilevel and discontinuous objective induced by weight-clipping-based policy optimization by reformulating the problem as a Heaviside composite optimization problem, which provides a rigorous computational framework. The reformulated policy optimization problem is then solved efficiently using the progressive integer programming method, making practical policy learning tractable. We establish an upper bound for the suboptimality of the proposed algorithm, which reveals how the reduction in MSE of policy value estimation, enabled by our proposed weight-clipping estimator, leads to improved policy learning performance.

This paper discusses the computational resolution and presents numerical results for solving affine combinations of Heaviside composite optimization problems (abbreviated as A-HSCOPs) by a progressive integer programming (abbreviated as PIP) method. The characteristics of these problems are that the Heaviside functions, which appear in the objective and define the constraints, are discontinuous, and their mixed-signed combinations result in the overall objective lacking the matching semicontinuity needed for the optimization and in the feasible set being not necessarily closed. Added to these challenging properties is the nondifferentiability of the inner functions in the composition. In this paper, we propose resolutions to all these challenges by first an approximation to remedy the lack of semicontinuity in the objective and closedness in the constraints, followed by a progressive integer programming approach with successive decomposition to handle the intrinsically discrete nature of the Heaviside function. Convergence to the local optimizers of the given Heaviside optimization problem is established. The effectiveness of the overall solution strategy is supported by extensive computational experiments on the score-based and tree-based multiclass classification problems with precision constraints.

Indefinite quadratic programs (QPs) are known to be very difficult to be solved to global optimality, so are linear programs with linear complementarity constraints (LPCCs). It is a classic result that for a QP with an optimal solution, the QP has an equivalent formulation as a certain LPCC in terms of their globally optimal solutions. Thus it is natural to attempt to solve an (indefinite) QP as a LPCC. This paper presents a progressive mixed integer linear programming method for solving a general LPCC. Instead of solving the LPCC with a full set of integer variables expressing the complementarity conditions, the presented method solves a finite number of mixed integer subprograms by starting with a small fraction of integer variables and progressively increasing this fraction. After describing the PIP (for progressive integer programming) method and providing some details for its implementation and tuning possibilities, we demonstrate, via an extensive set of computational experiments, the superior performance of the progressive approach over the direct solution of the full-integer formulation of the LPCCs in obtaining high-quality solutions. It is also shown that the solution obtained at the termination of the PIP method is a local minimizer of the LPCC, a property that cannot be claimed by any known non-enumerative method for solving this nonconvex program. In all the experiments, the PIP method is initiated at a feasible solution of the LPCC obtained from a nonlinear programming solver, and with high likelihood, can successfully improve it. Thus, the PIP method can improve a stationary solution of an indefinite QP, something that is not likely to be achievable by a nonlinear programming method. Finally, some analysis is presented that provides a better understanding of the roles of the LPCC suboptimal solutions in the local optimality of the indefinite QP. This local aspect of the connection between a QP and its LPCC formulation has seemingly not been addressed in the literature.

This paper explores Logarithmic Integral Optimization (LIO) problems, providing a unified computational framework for various tasks in computational statistics. Key among these are Maximum Likelihood Estimation (MLE) and Maximum a Posteriori (MAP) inference for probabilistic models. Specifically, we investigate scenarios where the model consists of conditional density functions with intractable normalizers. This feature can pose substantial computational challenges for the associated LIO, especially when coupled with the growing prevalence of nonconvex and nondifferentiable modelings in contemporary applications. To address these challenges, we propose an efficient algorithm for LIO, termed Adaptive Importance Sampling-based Surrogation. This method is designed to simultaneously handle nonconvexity and nondifferentiability, while also improving the sampling approximation of the intractable integral term in LIO through variance reduction. The justification of this algorithm is supported by our analysis, which establishes an almost sure subsequential convergence to a necessary candidate for a local minimizer, referred to as a surrogation stationary point. Furthermore, we demonstrate the effectiveness of our algorithm through extensive numerical experiments, confirming its efficiency and stability in facilitating more advanced probabilistic models with intractable normalizers.

This paper proposes a Heaviside composite optimization approach and presents a progressive method for solving multi-treatment learning problems with non-convex constraints. A Heaviside composite function is a composite of a Heaviside function (i.e., the indicator function of either the open ( 0,∞ ) or closed [ 0,∞ ) interval) with a possibly nondifferentiable function. Modeling-wise, we show how Heaviside composite optimization provides a rigorous mathematical formulation for learning multi-treatment rules subject to Gini constraints. A Heaviside composite function has an equivalent discrete formulation and the resulting optimization problem can in principle be solved by integer programming (IP) methods. Nevertheless, for constrained treatment learning problems with large datasets, a straightforward application of off-the-shelf IP solvers is usually ineffective in achieving global optimality. To alleviate such a computational burden, our major contribution is the proposal of the progressive method by leveraging the effectiveness of state-of-the-art IP solvers for problems of modest sizes. We provide the theoretical advantage of the progressive method with the connection to continuous optimization and show that the computed solution is locally optimal for a broad class of Heaviside composite optimization problems. The superior numerical performance of the proposed method is demonstrated by extensive computational experimentation. A brief discussion of how score-based and tree-based multi-classification problems can also be formulated as Heaviside composite optimization problems and thus treated by the same progressive method is presented in an appendix.

Quasi-differentiable functions were introduced by Pshenichnyi in a 1969 monograph written in Russian and translated in an English version in 1971. This class of nonsmooth functions was studied extensively in two decades since but has not received much attention in today's wide optimization literature. This regrettable omission is in spite of the fact that many functions in modern day applications of optimization can be shown to be quasi-differentiable. In essence, a quasi-differentiable function is one whose directional derivative at an arbitrary reference vector, as a function of the direction, is the difference of two positively homogenous, convex functions. Thus, to bring quasi-differentiable functions closer to the class of difference-of-convex functions that has received fast growing attention in recent years in connection with many applied subjects, we propose to rename quasi-differentiable functions as quasi-difference-convex (quasi-dc) functions. Besides modernizing and advancing this class of nonconvex and nondifferentiable functions, our research aims to put together a unified treatment of iterative convex-programming based descent algorithms for solving a broad class of composite quasi-dc programs and to establish their subsequential convergence, sequential convergence, and rates of convergence; the latter two topics are in line with the modern focus of such analysis for convex programs and some extensions and are departures from the sole emphasis of subsequential convergence in the traditional studies of quasi-differentiable optimization. Through this research, we have gained significant new insights and understanding, advanced the fundamentals, and broadened the applications of this neglected yet pervasive class of nonconvex and nondifferentiable functions and their optimization.

This paper proposes a Heaviside composite optimization approach and presents a progressive (mixed) integer programming (PIP) method for solving multi-class classification and multi-action treatment problems with constraints. A Heaviside composite function is a composite of a Heaviside function (i.e., the indicator function of either the open $( \, 0,\infty )$ or closed $[ \, 0,\infty \, )$ interval) with a possibly nondifferentiable function. Modeling-wise, we show how Heaviside composite optimization provides a unified formulation for learning the optimal multi-class classification and multi-action treatment rules, subject to rule-dependent constraints stipulating a variety of domain restrictions. A Heaviside composite function has an equivalent discrete formulation, and the resulting optimization problem can in principle be solved by integer programming (IP) methods. Nevertheless, for constrained learning problems with large data sets, a straightforward application of off-the-shelf IP solvers is usually ineffective in achieving global optimality. To alleviate such a computational burden, our major contribution is the proposal of the PIP method by leveraging the effectiveness of state-of-the-art IP solvers for problems of modest sizes. We provide the theoretical advantage of the PIP method with the connection to continuous optimization and show that the computed solution is locally optimal for a broad class of Heaviside composite optimization problems. The numerical performance of the PIP method is demonstrated by extensive computational experimentation.

We develop a general equilibrium model to capture the complex interactions between different modes, such as solo driving, public transit, as well as rideshare and ride-hailing services such as Uber and Lyft, under a joint morning and evening commute framework. Formulated as a variational inequality (VI) and equivalently as a mixed complementarity problem (MiCP), the model allows (a) travelers to switch between different transportation modes and (b) passengers from different Origin-Destination (OD) pairs to share a ride together. The computational results on the Sioux-Falls network show that our model captures the possible mode switches and the coupling effects between morning and evening commutes. Furthermore, our numerical examples demonstrate that modelling morning and evening commutes separately tends to overestimate the travelers' disutility and the average Vehicle Miles Traveled (VMT) in the network.

Refining and extending works by Ye and Kitahara-Mizuno, this paper presents new results on the number of pivots of simplex-type methods for solving linear programs of the Leontief kind, certain linear complementarity problems of the P kind, and nonnegatively constrained convex quadratic programs. Our results contribute to the further understanding of the complexity and efficiency of simplex-type methods for solving these problems. Two applications of the quadratic programming results are presented.

The minimization of nonlower semicontinuous functions is a difficult topic that has been minimally studied. Among such functions is a Heaviside composite function that is the composition of a Heaviside function with a possibly nonsmooth multivariate function. Unifying a statistical estimation problem with hierarchical selection of variables and a sample average approximation of composite chance constrained stochastic programs, a Heaviside composite optimization problem is one whose objective and constraints are defined by sums of possibly nonlinear multiples of such composite functions. Via a pulled-out formulation, a pseudostationarity concept for a feasible point was introduced in an earlier work as a necessary condition for a local minimizer of a Heaviside composite optimization problem. The present paper extends this previous study in several directions: (a) showing that pseudostationarity is implied by (and thus, weaker than) a sharper subdifferential-based stationarity condition that we term epistationarity; (b) introducing a set-theoretic sufficient condition, which we term a local convexity-like property, under which an epistationary point of a possibly nonlower semicontinuous optimization problem is a local minimizer; (c) providing several classes of Heaviside composite functions satisfying this local convexity-like property; (d) extending the epigraphical formulation of a nonnegative multiple of a Heaviside composite function to a lifted formulation for arbitrarily signed multiples of the Heaviside composite function, based on which we show that an epistationary solution of the given Heaviside composite program with broad classes of B-differentiable component functions can in principle be approximately computed by surrogation methods. Funding: The work of Y. Cui was based on research supported by the National Science Foundation [Grants CCF-2153352, DMS-2309729, and CCF-2416172] and the National Institutes of Health [Grant 1R01CA287413-01]. The work of J.-S. Pang was based on research supported by the Air Force Office of Scientific Research [Grant FA9550-22-1-0045].

This paper studies the existence of a (Lipschitz) continuous (single -valued) solution function of parametric variational inequalities under functional and constraint perturbations. At the most elementary level, this issue can be explained from classical parametric linear programming and its resolution by the parametric simplex method, which computes a solution trajectory of the problem when the objective coefficients and the right-hand sides of the constraints are parameterized by a single scalar parameter. The computed optimal solution vector (and not the optimal objective value) is a continuous piecewise affine function in the parameter when the objective coefficients are kept constant, whereas the computed solution vector can be discontinuous when the right-hand constraint coefficients are kept fixed and there is a basis change at a critical value of the parameter in the objective. We investigate this issue more broadly first in the context of an affine variational inequality (AVI) and obtain results that go beyond those pertaining to the lower semicontinuity of the solution map with joint vector perturbations; the latter property is closely tied to a stability theory of a parametric AVI and in particular to Robinson's seminal concept of strong regularity. Extensions to nonlinear variational inequalities is also investigated without requiring solution uniqueness (and therefore applicable to nonstrongly regular problems). The role of solution uniqueness in this issue of continuous single -valued solution selection is further clarified.

This paper studies several solution paths of sparse quadratic minimization problems as a function of the weighing parameter of the bi-objective of estimation loss versus solution sparsity. Three such paths are considered: the “ ℓ _0 -path” where the discontinuous ℓ _0 -function provides the exact sparsity count; the “ ℓ _1 -path” where the ℓ _1 -function provides a convex surrogate of sparsity count; and the “capped ℓ _1 -path” where the nonconvex nondifferentiable capped ℓ _1 -function aims to enhance the ℓ _1 -approximation. Serving different purposes, each of these three formulations is different from each other, both analytically and computationally. Our results deepen the understanding of (old and new) properties of the associated paths, highlight the pros, cons, and tradeoffs of these sparse optimization models, and provide numerical evidence to support the practical superiority of the capped ℓ _1 -path. Our study of the capped ℓ _1 -path is interesting in its own right as the path pertains to computable directionally stationary (= strongly locally minimizing in this context, as opposed to globally optimal) solutions of a parametric nonconvex nondifferentiable optimization problem. Motivated by classical parametric quadratic programming theory and reinforced by modern statistical learning studies, both casting an exponential perspective in fully describing such solution paths, we also aim to address the question of whether some of them can be fully traced in strongly polynomial time in the problem dimensions. A major conclusion of this paper is that a path of directional stationary solutions of the capped ℓ _1 -regularized problem offers interesting theoretical properties and practical compromise between the ℓ _0 -path and the ℓ _1 -path. Indeed, while the ℓ _0 -path is computationally prohibitive and greatly handicapped by the repeated solution of mixed-integer nonlinear programs, the quality of ℓ _1 -path, in terms of the two criteria—loss and sparsity—in the estimation objective, is inferior to the capped ℓ _1 -path; the latter can be obtained efficiently by a combination of a parametric pivoting-like scheme supplemented by an algorithm that takes advantage of the Z-matrix structure of the loss function.

This paper investigates computational algorithms for Maximum a Posteriori (MAP) and Maximum Likelihood Estimation (MLE) inference of Bayesian Hierarchical Models (BHMs), via a unified formulation as a nonconvex and nondifferentiable logarithmic integral optimization problem. Specifically, we explore situations where a BHM comprises density functions with intractable normalizers, a feature that can present significant computational obstacles, particularly when combined with nonconvexity and nondifferentiability, which are increasingly prevalent in contemporary applications of computational statistics. To deal with these challenges, we propose an efficient algorithmic approach, termed Adaptive Importance Sampling-based Surrogation , to simultaneously handle nonconvexity and nondifferentiability while also improving the sampling approximation of the intractable normalizer through variance reduction. Performance of this al-gorithm is guaranteed by our analysis which establishes an almost sure subsequential convergence to a necessary candidate for a local minimizer, referred to as a surrogation stationary point . We also demonstrate the effectiveness of our algorithm through extensive numerical experiments, verifying its efficiency and stability in enabling more advanced BHMs where intractable normalizers arise as a result of enhanced modeling capability.

This paper studies several versions of the sparse optimization problem in statistical estimation defined by a pairwise separation objective. The sparsity (i.e., ℓ _0 ) function is approximated by a folded concave function; the pairwise separation gives rise to an objective of the Z-type. After presenting several realistic estimation problems to illustrate the Z-structure, we introduce a linear-step inner-outer loop algorithm for computing a directional stationary solution of the nonconvex nondifferentiable folded concave sparsity problem. When specialized to a quadratic loss function with a Z-matrix and a piecewise quadratic folded concave sparsity function, the overall complexity of the algorithm is a low-order polynomial in the number of variables of the problem; thus the algorithm is strongly polynomial in this quadratic case. We also consider the parametric version of the problem that has a weighted ℓ _1 -regularizer and a quadratic loss function with a (hidden) Z-matrix. We present a linear-step algorithm in two cases depending on whether the variables have prescribed signs or with unknown signs. In both cases, a parametric algorithm is presented and its strong polynomiality is established under suitable conditions on the weights. Such a parametric algorithm can be combined with an interval search scheme for choosing the parameter to optimize a secondary objective function in a bilevel setting. The analysis makes use of a least-element property of a Z-function, and, for the case of a quadratic loss function, the strongly polynomial solvability of a linear complementarity problem with a hidden Z-matrix. The origin of the latter class of matrices can be traced to an inspirational paper of Olvi Mangasarian to whom we dedicate our present work.

For the treatment of outliers, the paper “Risk-Based Robust Statistical Learning by Stochastic Difference-of-Convex Value-Function Optimization” by Junyi Liu and Jong-Shi Pang proposes a risk-based robust statistical learning model. Employing a variant of the conditional value-at-risk risk measure, called the interval conditional value-at-risk (In-CVaR), the model aims to exclude the risks associated with the left and right tails of the loss. The resulting nonsmooth and nonconvex model considers the population In-CVaR risk and distinguishes the upside and downside losses with asymmetric weights. For the solution of the model in both regression and classification, the authors show that the objective function is the difference of two convex functions each being the optimal objective value of a univariate convex stochastic program. A sampling and convex programming-based algorithm is developed with the appropriate control of incremental sample sizes, and its subsequential almost-sure convergence to a critical point is established. Numerical results illustrate the practical performance of the model and methodology.

This paper begins with a class of convex quadratic programs (QPs) with bounded variables solvable by the parametric principal pivoting algorithm with $\mathcal{O}(n^3)$ strongly polynomial complexity, where $n$ is the number of variables of the problem. Extension of the Hessian class is also discussed. Our research is motivated by a recent reference [7] wherein the efficient solution of a quadratic program with a tridiagonal Hessian matrix in the quadratic objective is needed for the construction of a polynomial-time algorithm for solving an associated sparse variable selection problem. With the tridiagonal structure, the complexity of the QP algorithm reduces to $\mathcal{O}(n^2)$. Our strongly polynomiality results extend previous works of some strongly polynomially solvable linear complementarity problems with a P-matrix [9]; special cases of the extended results include weakly quasi-diagonally dominant problems in addition to the tridiagonal ones.

There are many significant applied contexts that require the solution of discontinuous optimization problems in finite dimensions. Yet these problems are very difficult, both computationally and analytically. With the functions being discontinuous and a minimizer (local or global) of the problems, even if it exists, being impossible to verifiably compute, a foremost question is what kind of ''stationary solutions'' one can expect to obtain; these solutions provide promising candidates for minimizers; i.e., their defining conditions are necessary for optimality. Motivated by recent results on sparse optimization, we introduce in this paper such a kind of solution, termed ''pseudo B- (for Bouligand) stationary solution'', for a broad class of discontinuous piecewise continuous optimization problems with objective and constraint defined by indicator functions of the positive real axis composite with functions that are possibly nonsmooth. We present two approaches for computing such a solution. One approach is based on lifting the problem to a higher dimension via the epigraphical formulation of the indicator functions; this requires the addition of some auxiliary variables. The other approach is based on certain continuous (albeit not necessarily differentiable) piecewise approximations of the indicator functions and the convergence to a pseudo B-stationary solution of the original problem is established. The conditions for convergence are discussed and illustrated by an example.

This paper studies a structured compound stochastic program (SP) involving multiple expectations coupled by nonconvex and nonsmooth functions. We present a successive convex programming-based sampling algorithm and establish its subsequential convergence. We describe stationary properties of the limit points for several classes of the compound SP. We further discuss probabilistic stopping rules based on the computable error bound for the algorithm. We present several risk measure minimization problems that can be formulated as such a compound stochastic program; these include generalized deviation optimization problems based on the optimized certainty equivalent and buffered probability of exceedance (bPOE), a distributionally robust bPOE optimization problem, and a multiclass classification problem employing the cost-sensitive error criteria with bPOE.