Sinkhorn's algorithm is a simple yet powerful method that has found wide applications, most notably in computational optimal transport. More recently it has featured prominently in the training of large language models such as DeepSeek. Despite the many previous works that study its convergence behavior, a tight non-asymptotic local convergence analysis remains absent. In this paper, we solve this open problem. We show that the local convergence rate Sinkhorn algorithm is $O(1/\lambda_2 \log(1/\varepsilon))$, where $\lambda_2$ is the second largest eigenvalue of a matrix associated with the original problem data. Extensive numerical experiments validate the tightness of our result. Our analysis also shows an implicit connection between the Sinkhorn algorithm and the preconditioned gradient descent, and with this connection, we show that an accelerated $O(1/\sqrt{\lambda_2} \log(1/\varepsilon))$ local convergence is achievable through acceleration.

In electricity grid capacity expansion planning, representative days are used to model daily and hourly operating conditions for long-term investment planning. Common approaches to representative day selection rely on clustering algorithms (e.g., k-means). However, these methods may fail to capture operating conditions on “extreme days”, which are the days when peak loads occur. Failure to model these conditions can negatively impact system resilience. In particular, there is a need to effectively capture the operational and structural impact to the grid of compound hazards, as they may be more severe than individual hazards in isolation. We present a probabilistic approach for selecting extreme days for capacity expansion planning that explicitly considers the spatio-temporal relationships among compound hazards.

In several scientific and engineering fields, well-defined models exist on a spectrum of fidelity with an inherent zero sum trade-off between accuracy and problem complexity. In many stochastic programming (SP) applications, computational tractability is a major challenge that in practice limits either the space of exploration, the fidelity of model used, or both. Decomposition methods such as progressive hedging can reduce CPU cost through parallelization but still requires choosing a single model fidelity. In this work, we detail our framework for creating multifidelity representations of stochastic problems and novel bundling strategies of scenarios that can accelerate search and convergence in large-scale SPs, while leveraging both low and high fidelity model representations. We demonstrate our methods on a diverse set of SP applications, including scheduling, grid operations and planning, and chemical process control in order to explore the effects of problem linearity, discrete/continuous variables, bundling strategies, and algorithmic hyperparameters.

Two-stage stochastic power planning problems are central to modeling uncertainty in generation capacity expansion and investment decisions. These problems are often large-scale and potentially nonconvex, motivating the use of decomposition-based global optimization methods. In this work, we introduce a priority branch-and-bound (B&B) framework for decomposable nonlinear optimization, implemented in our open-source Python package, SNoGloDe (Structured Nonlinear Global Decomposition). The algorithm exploits the natural scenario-based separability of stochastic formulations: each realization of uncertain parameters defines a scenario subproblem, and solving these subproblems separately (ideally in parallel) provides a relaxation that serves as a valid bound for the original coupled problem.

Unlike conventional B&B algorithms that minimize node evaluations to close the optimality gap, our approach identifies diverse, high-quality alternative solutions within the search tree. In the context of power systems planning, such alternative solutions offer operators valuable insights into trade-offs among cost, reliability, and resilience under uncertainty.

SNoGloDe provides a flexible software framework for embedding these algorithms, integrating parallel computation, global optimization solvers, and customizable branching logic. This enables systematic exploration of the feasible space and supports the generation of multiple near-optimal solutions. Our results demonstrate that structured global decomposition can both accelerate solution times and enhance interpretability in stochastic power planning optimization.

A fundamental problem in optimization is finding an $\epsilon$-first-order stationary point of a smooth nonconvex function using only gradient information. The best-known gradient complexity for this task, assuming both the gradient and Hessian of the objective function are Lipschitz continuous, is $O(\epsilon^{-1.75})$. In this talk, I will present a new method with a gradient complexity of $O(d^{0.25}\epsilon^{-1.625})$, where $d$ is the problem dimension, yielding an improvement when $d = O(\epsilon^{-0.5})$.

Our key idea is to design a quasi-Newton method that operates through two levels of online learning. Specifically, we first reformulate the task of finding a stationary point for a nonconvex problem as minimizing the regret in an online convex optimization problem, where the loss is determined by the gradient of the objective function. Then, we introduce a novel optimistic quasi-Newton method to solve this online learning problem, with the Hessian approximation update itself framed as an online learning problem in the space of matrices. Beyond improving the known complexity bounds, our result provides the first guarantee that quasi-Newton methods can outperform gradient descent-type methods in nonconvex optimization.

There is emerging evidence that trust-region (TR) algorithms are very effective at solving derivative-free nonconvex stochastic optimization problems in which the objective function is a Monte Carlo (MC) estimate. A recent strand of methodologies adaptively adjusts the sample size of the MC estimates by keeping the estimation error below a measure of stationarity induced from the TR radius. In this work we explore stratified adaptive sampling strategies to equip the TR framework with accurate estimates of the objective function, thus optimizing the required number of MC calls to reach a given ϵ-accuracy of the solution. We prove a reduced sample complexity and explore inexpensive implementations in high dimension.

This study presents a derivative-free solver for multi-region stochastic equilibrium problems motivated by insurance market dynamics under climate change risks. The proposed framework models insurers competing across multiple regions and transferring a portion of their catastrophe risk to reinsurers, all under uncertain hurricane losses. Since the underlying equilibrium conditions are simulation-based and lack analytic gradients, a derivative-free approach is developed to efficiently search for equilibria. Numerical experiments using coastal North Carolina data demonstrate that the proposed method reliably computes equilibrium outcomes under climate-induced uncertainty within a reasonable computational budget.

Despite significant theoretical and computational advances in two-stage stochastic programming (SP), most existing research assumes that the uncertain parameters (or their probability distributions) are known a priori or belong to a predefined set. In particular, these studies typically assume that the distribution of random variables is non-adaptive to decisions. While such assumptions simplify mathematical formulation and computation, they seldom hold in practice. In this talk, we present a novel Bayesian learning framework for two-stage SP, in which the underlying distribution of random parameters is adaptive and can be progressively learned from the observed effects of decisions. We demonstrate this paradigm through a capacitated facility location problem under demand uncertainty. Furthermore, we propose a solution methodology for the resulting model and provide a theoretical analysis on the sensitivity of the optimal learning policy with respect to sampling costs and prior information.

We study a contextual stochastic optimization problem, where uncertain parameters and covariates are observed jointly, and the goal is to choose decisions minimizing expected cost given new covariates. Existing regression-based approaches assume that the prediction model fully captures the relationship between uncertainty, covariates, and decisions—yet this relationship is fundamentally unknown, a challenge known as model uncertainty. To address this, we propose a model uncertainty–aware residuals-based SAA framework (MUR-SAA), formulated as a distributionally robust optimization problem. By treating the regression function as a random element and optimizing against the worst-case expected cost over all plausible distributions of the unknown regression function, MUR-SAA produces solutions that are robust to model misspecification. We establish, for the first time, efficient mechanisms for quantifying the stability of general residual-based stochastic programs and our MUR-SAA model. Moreover, we identify conditions under which our MUR-SAA model enjoys desirable asymptotic convergence and finite-sample properties. We derive equivalent tractable reformulations of MUR-SAA and develop a decomposition algorithm to solve them. We illustrate our theoretical results and demonstrate the advantages of our MUR-SAA model using a contextual portfolio optimization problem.

Recent advances in machine learning have enabled data-driven proxies that approximate optimization models using historical or synthetic data. These models offer rapid, near-optimal solutions, making them attractive for time-sensitive applications. However, they often lack explicit representations of system constraints, resulting in infeasible or unreliable decisions in safety-critical contexts.

Tuned low-fidelity optimization models have emerged as a promising alternative by combining historical information with known system dynamics and constraints. While more interpretable and reliable, their training typically requires solving expensive bilevel optimization problems to identify optimal tuning parameters.

This work introduces a self-supervised learning framework for tuning low-fidelity optimization models using differentiable optimization layers. The approach integrates three components: a learning model that predicts tunable parameters, a differentiable optimization layer mapping these parameters to approximate decisions, and a high-fidelity model that implicitly evaluates solution quality. Applied to the AC Optimal Power Flow problem, with the DC-OPF model as its low-fidelity counterpart, the proposed method achieves superior data efficiency, training speed, and solution accuracy compared to conventional supervised approaches.

We design a machine learning (ML)-guided enhancement of the Undercover primal heuristic for mixed-integer nonlinear programs (MINLPs). Undercover identifies a subset of variables to fix and assigns them values, thereby restricting the MINLP to a mixed-integer linear program (MILP) that can be solved using off-the-shelf solvers to potentially obtain a feasible solution. While Undercover performs well on mixed-integer quadratically constrained quadratic programs (MIQCQPs), its effectiveness on general MINLPs is limited by the difficulty of choosing good fixing values. To address this challenge, we propose a graph convolutional network (GCN) that predicts high-quality fixing values for selected variables. The model is pretrained using best-known feasible solutions and further optimized through an end-to-end loss function by leveraging LP sensitivity theory. We integrate our ML-guided Undercover heuristic within the global solver SCIP and evaluate it on benchmark instances from QPLIB and MINLPLIB. Computational results demonstrate that the proposed approach substantially improves the performance of the default Undercover implementation in SCIP.

Physics-informed neural networks (PINNs) do not always guarantee strict constraint satisfaction. This is problematic in engineering systems where minor violations of governing laws can degrade the reliability and consistency of model predictions. To overcome these challenges, we introduce KKT-Hardnet (https://github.com/SOULS-TAMU/kkt-hardnet), a neural network architecture that enforces linear and nonlinear equality and inequality constraints to machine precision. It leverages a differentiable projection onto the feasible region by solving Karush-Kuhn-Tucker (KKT) conditions of a distance minimization problem. Our major contributions are twofold (1) we provided a unique strategy to include the inequality constraints by considering the complementary conditions through Fischer-Burmeister reformulation to ensure nonnegativity in corresponding Lagrange multipliers, (2) we developed the framework for hard nonlinear constraint satisfaction for both equality and inequality. Compared to multilayer perceptrons (MLPs) and PINNs, KKT-Hardnet achieves strict constraint satisfaction. It also circumvents the need to balance data and physics residuals in PINN training. This enables the integration of domain knowledge into machine learning towards reliable hybrid modeling of complex systems. Moreover, our framework allows us to train neural networks to learn and predict solutions to nonlinear optimization problems, offering a new avenue for integrating machine learning with optimization.

We introduce a sequential version of the minimum s-t cut problem, defined by a directed graph with source s and sink t, and nonnegative random variables for each arc representing arc weights. Starting with S = {s} , we observe weight realizations for outgoing arcs from S, and choose to either stop, or to add a node to S, with the objective of minimizing the expected weight of the chosen cut. When the graph is a directed s-t path, the model reduces to the classical prophet inequality problem in minimization form.

We formulate the problem as a dynamic program (DP) and show that computing the expected cost of an optimal policy is #P-hard even on simple series-parallel graphs. We leverage a linear programming (LP) formulation of the DP to construct a separable cost-to-go function approximation (CFA), yielding a lower bound. We design a heuristic policy based on the CFA and provide a guarantee on its worst-case performance that depends on graph topology but is independent of arc weight distributions.

As political redistricting in the United States continues to suffer from partisan gerrymandering, reform proposals proliferate. Several reform ideas fall into the category of two-player redistricting games. Inspired by the famous "I cut, you choose" protocol for fair division of a cake between two agents, these proposals turn drawing a state's districts into a two-player (i.e., two-party) extensive form game. Throughout such a game, each party controls some of the final voting district boundaries. The hope is that Nash equilibria yield fairer district plans than those proffered by either party. To find analytical solutions of Nash equilibria, proposers of redistricting games often relax geographic constraints, switching from a balanced, connected graph partitioning problem to one of rearranging voters arbitrarily into bins as if red and blue marbles.

To what degree will fairness results from this geography-free setting translate to reality? In this talk, we show that several redistricting games have pathological examples (despite strong restrictions) in which the problem’s geography gives one player an outsized advantage. We then apply MIP-based heuristics to simulate redistricting games for more realistic instances. Applying standard partisan fairness measures, we compare our resulting district plans to those found by other heuristics. We conclude with policy implications and possible extensions.

In this talk, we present new findings on facial reduction for semidefinite relaxations of combinatorial optimization problems. In semidefinite programming (SDP), Slater’s condition is crucial for both theoretical convergence guarantees and the practical performance of optimization algorithms. When Slater’s condition fails, facial reduction can restore it through a finite sequence of reformulations. However, these reformulations often involve solving auxiliary optimization problems that can be as challenging as the original. Recent research has therefore focused on developing more efficient strategies for performing facial reduction.

In our work, we specifically consider SDP problems that arise as relaxations of combinatorial optimization problems. This perspective enables us to exploit the underlying combinatorial structure, allowing the development of novel and highly efficient facial reduction techniques. We establish theoretical results demonstrating the effectiveness of our approach. Numerical experiments further show that applying our specialized facial reduction method significantly improves both the speed and accuracy of solving SDP problems.

Sums-of-squares (SOS) and other popular convex relaxations of polynomial minimization problems are often solved using off-the-shelf numerical methods which do not immediately yield rigorously provable lower bounds on the global minimum. For instance, the moment-SOS hierarchy yields rigorous upper (instead of lower) bounds on the sought lower bound, while the numerical solution of the SOS SDP (the dual of the moment-SOS SDP) yields a rigorous nonnegativity certificate of a slightly different polynomial than the one at hand, due to the inexactly satisfied equality constraints in the SDP.

More generally, for conic programs in standard form, numerically computed primal "feasible" solutions usually aren't feasible, and do not provide rigorous upper bounds on the optimal value.

The first talk of this two-part series presents a theoretical framework (termed "dual certificates") for obtaining rigorous upper bounds from suitable dual feasible solutions, complementing the standard duality theory in which the same dual solutions yield lower bounds. When applied to tractable subcones of nonnegative polynomials commonly encountered in polynomial optimization (SOS, SONC, SAGE and SDSOS polynomials, etc.), dual certificates enable the computation of rigorous nonnegativity certificates using inexact, numerical methods.

Dual membership certificates are easily verifiable, rigorous certificates of membership in a cone represented by vectors in the dual cone. This construction complements the classical interpretation of the dual cone, in which such vectors serve as certificates of non-membership. An important feature of the approach is that exact certificates can be computed using standard numerical methods.

The second talk of this two-part series is concerned with the computational aspects of dual certificates in optimization problems. First, we show that dual certificates can be obtained entirely numerically, for example by using primal-dual interior-point algorithms (IPAs). We then show that, using dual certificates, we can recover primal feasible solutions from suitable dual feasible ones. This enables the computation of both rigorously certifiable lower and upper bounds from the same dual strictly feasible solution. The gap between these bounds tends to zero as the dual solution approaches optimality.

Dual certificates, in turn, can be used to enhance existing IPAs: one may, e.g., "retrofit" dual methods to return exact feasible solutions, or certify the progress of IPAs. We present numerical examples for the theoretical results.

We present a systematic benchmarking study of state-of-the-art derivative-free optimization algorithms on a broad suite of test problems. We compare a number of global and local solvers, including customary stochastic methods and modern deterministic methods like the recently proposed branch-and-model algorithm. We evaluate optimization performance using several metrics such as scalability and sensitivity to initialization. The comparison demonstrates the strengths and limitations of each algorithm across diverse problem classes.

Electoral districts can be modeled as partitions of a network representing the underlying geographic units, which allows for analysis using discrete optimization and sampling methods. Spanning tree approaches have become popular for generating large samples of comparison plans in this setting, which can introduce a strong preference for partitions with many spanning trees in the target distribution. In this talk I will discuss recent applications and extensions of these techniques, including for court cases, support of line-drawing commissions, and evaluating nonpartisan justifications for proposed plans. Along the way I will present related open problems and proposals based on the cycle basis walk and classic network partitioning methods including spectral clustering and modularity.

When auditing a redistricting plan, one often wants to compare what *would* have happened under non-partisan plans. In accessing this counterfactual, one first places a probability distribution on the space of redistricting plans and then uses Monte Carlo sampling to integrate the space and determine features of interest (e.g., how many party members does any given map elect?)

Specifically, a class of methods called Recombination (or Recom) has become a popular Monte Carlo method to use in practice, however, it is only (empirically) capable of efficiently sampling from certain classes of measures. In this talk we will use spectral methods to understand these barriers to mixing on an extended class of measures. I will conclude with some discussion on a survey of alternative methods that extend the class of measures that we can efficiently sample and avoid the barriers found in Recom.

In recent months, politicians across the spectrum have complained about maps drawn by their political opponents that slice-and-dice communities to gerrymander seats for their party. One proposed reform is to explicitly require mapmakers to seek proportionality: the idea that if a party earns X% of the votes, then they should win roughly X% of the seats. However, some have expressed concern that achieving proportionality will require non-compact districts that shred counties. In this talk, we propose optimization methods to test whether this is true. Specifically, we ask whether it is possible to draw maps that achieve partisan proportionality and minimize the extent to which counties are divided across districts. For many states, we find that the answer is yes, i.e., there exist maps that achieve ideal performance on both objectives, simultaneously. Further, the generated maps are more compact than those currently in use. We accomplish this with a combination of heuristic and exact integer programming methods (e.g., Benders decomposition). Our results are practically relevant for proposed anti-gerrymandering legislation and inform independent redistricting commissions and special masters who must delicately balance competing criteria.

Thompson Sampling is one of the most widely used and studied bandit algorithms, known for its simple structure, low regret performance, and solid theoretical guarantees. Yet, in stark contrast to most other families of bandit algorithms, the exact mechanism through which posterior sampling (as introduced by Thompson) is able to "properly" balance exploration and exploitation, remains a mystery. In this talk we show that the core insight to address this question stems from recasting Thompson Sampling as an online optimization algorithm. To distill this, a key conceptual tool is introduced, which we refer to as "faithful" stationarization of the regret formulation. Essentially, the finite horizon dynamic optimization problem is converted into a stationary counterpart which "closely resembles" the original objective (in contrast, the classical infinite horizon discounted formulation, that leads to the Gittins index, alters the problem and objective in too significant a manner). The newly crafted time invariant objective can be studied using Bellman's principle which leads to a time invariant optimal policy. When viewed through this lens, Thompson Sampling admits a simple online optimization form that mimics the structure of the Bellman-optimal policy, and where greediness is regularized by a measure of residual uncertainty based on point-biserial correlation. This answers the question of how Thompson Sampling balances exploration-exploitation, and moreover, provides a principled framework to study and further improve Thompson's original idea.

We present the recent development of GPU-accelerated optimization solvers and optimization differentiators. Optimization layers in neural networks—nonlinear layers where the parameters of a convex optimization problem are learned by differentiating a function of the optimal solution—were developed at the end of the 2010s. At that time, their real-world applications were limited by the performance of CPU-bound optimization solvers and differentiators. GPU-acceleration of convex optimization problems had been solved only for special cases like dense QPs and cVar-constrained optimization. In the last two years, general conic interior point methods began to achieve speedups on GPUs, thanks to advancements in sparse matrix factorization. Coupled with new approaches to argmin differentiating, CVXPYlayers---a modeling language for differentiable optimization---developed support for GPU-acceleration, unlocking new performance opportunities for optimization layers in neural networks. We present recent advancements in the space and discuss applications unlocked by the newly achieved high-performance solvers.

Multistage robust multiobjective optimization addresses multistage decision-making problems under multiple conflicting objectives and uncertainty that results from the volatility of the environment. Although robust multiobjective optimization has been extensively studied, difficulties and complexities of solving multistage robust multiobjective optimization problems remain underexplored.
We propose a two-stage robust multiobjective optimization problem (TSRMOP) with uncertainty in general locations. Relying on policy optimization and scalarization of the second-stage problem, we develop a max-min property and reformulate the TSRMOP into a single-stage, single-objective problem whose optimal solutions are two-stage robust ε-properly efficient solutions to the original TSRMOP. By employing an achievement scalarizing function, the proposed approach is readily applicable to interactive decision-support.
We apply this solution approach to a real-life forest management problem that is solved interactively. By iteratively updating preferences and solving the scalarized problem, we find a two-stage robust ε-properly efficient solution that is the most preferred by the decision maker. A comparison of the computed two-stage solutions against the solutions computed with a single-stage robust multiobjective model on the same data reveals the advantage of two-stage modelling.

Multistage multiobjective optimization addresses decision problems under conflict and uncertainty allowing for decision making in stages to suit changing environments or currently available information. A two-stage robust biobjective linear program (TSRBOLP) models two-stage decision processes with strategic decisions being made first, and operational decisions being made at later time after more information about the decision problem becomes known. The goal is to compute the first-stage feasible solutions that are Pareto-efficient with respect to the first and second-stage objectives and account for the worst-case uncertainty scenario.

The TSRBOLP is studied under discrete or continuous uncertainty in the right-hand-side of the constraints and with one or two second-stage objectives. The weighted-sum scalarization transforms the TSRBOLP into a single-objective optimization problem (SOP) whose optimal solutions are Pareto-efficient to the TSRBOLP. The SOP assumes the form of a parametric linear program with bilinear terms in the constraints. A parametric Benders’ decomposition algorithm is designed to compute the exact or approximate parametric optimal solutions to the SOP. The algorithm relies on an adaptation of Benders’ decomposition into a parametric optimization environment. The implementation of the algorithm makes use of MATLAB Multiparametric Toolbox and the Biconvex Approximate Simplex Method (Pangia, 2024). Numerical examples are included.

We study a two-stage stochastic multi-objective linear program (TSSMOLP) under risk neutral formulation, which is a natural generalization of a two-stage stochastic linear program (TSSLP). The proposed TSSMOLP formulation allows modelers to specify multiple objectives in each stage. The second-stage recourse decision is modeled by a random multi-objective linear program (MOLP) whose solution maps to a random second-stage nondominated set. The TSSMOLP then comprises the objective, which is the Minkowski sum of a linear term plus the expected value of second-stage nondominated set, and the constraints, which are linear. We prove properties of TSSMOLPs, which are stochastic set optimization problems, and the multifunctions that arise therein. To aid theoretical analysis, we propose nondominance-equivalent reformulations to the original problem and show that (i) the global Pareto set of a TSSMOLP is cone-convex, (ii) solving a parameterized TSSLP produces a weakly efficient point for the TSSMOLP, and (iii) a TSSMOLP in which the uncertainty is governed by an atomic probability measure with finite support can be reformulated as an MOLP. Leveraging the above properties we demonstrate that solving a TSSMOLP under a sample average approximation framework is equivalent to solving a large scale MOLP, and validate our results through a numerical example.

Least-squares regression is typically formulated as a quadratic program. This talk presents a novel approach for reducing it to a piecewise linear convex minimization problem within the Risk Quadrangle (RQ) Framework. Evidently, this problem can be reduced to linear programming. Crucially, this is not a heuristic step: the linearized formulation is statistically justified and shown to be equivalent to quantile regression, estimating the specific quantile corresponding to the mean. As a result, it inherits many benefits of quantile-based approaches, such as effective handling of heteroskedastic data and robustness to outliers.

We conduct extensive numerical experiments to justify our theoretical claims and highlight the computational benefits of our method. Our results highlight the versatility of the RQ Framework in unifying and extending classical regression methods, opening the door to broader applications and deeper theoretical understanding.

This work presents a formal framework for analyzing iterative optimization algorithms. We develop a general method for deriving bounds on sequences defined by recursive inequalities, based on constructing a majorizing system of equalities. This approach provides a foundation for approximating the divergence rates of Lipschitz functions. We illustrate its effectiveness by applying it to establish convergence bounds for various momentum-based methods, including SPSTORM algorithm and Nesterov’s acceleration.

This study presents a relatively robust pricing framework for multiproduct monopolies under distribution-free demand uncertainty characterized by Frobenius-norm constraints. Compared to extant models based on interval uncertainty, our approach employs ambiguity sets that align naturally with standard multivariate estimation techniques. The optimal robust price vector, obtained as the solution to a quasi-concave programming problem, maximizes the worst-case ratio of achieved profit to the ex-post optimal profit under full information. A key structural result is that the performance index is strictly quasi-concave in price, allowing the unique optimum to be characterized by a first-order condition. An embedded fixed-point problem can be solved efficiently, and closed-form expressions are derived for symmetric linear, semi-log, and log-log demand systems with associated parameter uncertainty. In general nonlinear settings, the model remains applicable through local linearization and regression-based confidence regions. To evaluate the approach, we conduct both extensive simulation experiments and empirical analyses using large-scale datasets from a major online retailer. The results demonstrate the superior performance of the relatively robust pricing method compared to alternative robustness criteria such as maximin payoff, minimax regret, and Laplacian uncertainty.

We consider optimization problems containing nonconvex quadratic functions for which semidefinite programming (SDP) relaxations often yield strong bounds. We investigate linear inequalities that outer approximate the positive semidefinite cone and are sparse in the sense that they are supported only on the variables corresponding to products of variables that are present in quadratic functions. We show that these sparse linear inequalities yield an LP relaxation that gives the same bound as the SDP relaxation. We demonstrate how to identify these inequalities via a separation procedure that involves solving a structured ``projection'' SDP. In a computational study, we find that the sparse LP relaxations defined by these inequalities can significantly accelerate branch-and-bound methods for globally solving nonconvex optimization problems.

This work presents a new convexification framework that reformulates the single-row relaxations of chance-constrained programs as bilinear sets over a simplex in a lifted space. Using an aggregation-based procedure, we derive facet-defining inequalities in the original space of variables, unifying and extending known families of valid inequality while introducing new ones that capture larger portions of the convex hull. Main contributions include a general aggregation-based convexification method for bilinear sets over a simplex, a bilinear reformulation of mixing sets with knapsack constraints, and a unified description of existing and new inequality families. Computational results show that the proposed inequalities recover over 90% of convex hull facets, yielding stronger relaxations and improved efficiency.

Transit Network Design is a well-studied problem in the field of transportation, typically addressed by solving optimization models under fixed demand assumptions. Considering the limitations of these assumptions, this paper proposes a new framework, namely the Two-Level Rider Choice Transit Network Design (2LRC-TND), that leverages machine learning and contextual stochastic optimization (CSO) to incorporate two layers of demand uncertainties into the network design process. The first level identifies travelers who rely on public transit (core demand), while the second level captures the conditional adoption behavior of those who do not (latent demand), based on the availability and quality of transit services. To capture these two types of uncertainties, 2LRC-TND relies on two travel mode choice models, that use multiple machine learning models. To design a network, 2LRC-TND integrates the resulting choice models into a CSO, formulated using two approaches: Mixed-Integer Programming (MIP) and Constraint Programming (CP).

2LRC-TND is evaluated through a case study involving over 6,600 travel arcs and 38,000 trips in the Atlanta metropolitan area. Under different budget and ridership scenarios, the computational results demonstrate the effectiveness of the 2LRC-TND in designing transit networks that account for demand uncertainties and contextual information, offering a more realistic alternative to fixed-demand models. Additionally, under different experimental scenarios, not all MIP and CP models can be solved to optimality. This work also compares the effectiveness of MIP and CP in addressing the problem and discusses the ongoing development of customized solution algorithms.

We present a hybrid low-rank augmented Lagrangian method for large-scale linear programs. Each outer iteration solves a convex quadratic subproblem with an Enhanced Frank–Wolfe routine that alternates a support-restricted accelerated composite gradient refinement with a global Frank–Wolfe step, warm-started from the previous outer iterate. The subproblem is the augmented Lagrangian, with full multiplier updates and a doubling schedule for $\beta_t$ until a joint residual-and-gap test terminates. The theory gives an $O(1/\bar\varepsilon)$ inner rate and a logarithmic bound on outer iterations via a potential $\Phi(\varepsilon_p,\varepsilon_g)$ (where $\varepsilon_p$ is the primal-feasibility tolerance and $\varepsilon_g$ is the primal-dual gap tolerance). For general LPs over the nonnegative orthant, an adaptive-radius Phase-1/Phase-2 scheme doubles the radius until the complementarity surrogate is small, guaranteeing finite-phase convergence and explicit bounds on total Frank–Wolfe steps. Empirically, the method is competitive with the simplex method, the interior-point method, PDLP, and other first order methods on synthetic and benchmark datasets, with both CPU and GPU implementations compared.

We introduce a verification framework to exactly verify the worst-case performance of sequential convex programming (SCP) algorithms for parametric non-convex optimization. The verification problem is formulated as an optimization problem that maximizes a performance metric (e.g., the suboptimality after a given number of iterations) over parameters constrained to be in a parameter set and iterate sequences consistent with the SCP update rules. Our framework is general, extending the notion of SCP to include both conventional variants such as trust-region, convex-concave, and prox-linear methods, and algorithms that combine convex subproblems with rounding steps, as in relaxing and rounding schemes. Unlike existing analyses that may only provide local guarantees under limited conditions, our framework delivers global worst-case guarantees—quantifying how well an SCP algorithm performs across all problem instances in the specified family. Applications in control, signal processing, and operations research demonstrate that our framework provides, for the first time, global worst-case guarantees for SCP algorithms in the parametric setting.

As modern decision-making problems continue to grow in scale and structural complexity, there is a pressing need for optimization frameworks that remain both computationally efficient and reliable. Classical model-based distributed optimization is a powerful paradigm for large-scale decision-making, yet it typically requires extensive tuning, faces computational and communication limitations, and has limited ability to exploit data or generalize across problem instances. On the other hand, black-box deep neural networks have witnessed widespread success in various domains, but they often lack interpretability and performance guarantees. To leverage their complementary strengths, we introduce a new class of deep learning-aided distributed architectures for large-scale convex optimization. The proposed frameworks rely on unrolling the iterations of distributed optimizers as layers of a deep learning network, and then incorporating learnable components. Extensive experiments demonstrate substantial improvements over traditional optimization algorithms in wall-clock time across diverse applications including optimal control, linear regression and transportation networks. Critically, they exhibit strong generalization: trained on small problem instances, they scale to solve significantly larger ones using the same learned policies. We also provide certifiable performance guarantees on the outputs, ensuring high-quality solutions for unseen problem instances.

Conventional statistical inference methods are typically developed for models simple enough to admit tractable estimators through carefully designed iterative algorithms. In contrast, modern deep learning models are enormously complex, yet are trained by simple gradient-descent-type algorithms, often without any provable guarantee of algorithmic convergence to global/local optima.

Can we reconcile classical inference principles with these highly complicated, modern learning paradigms? In this talk, we will present a new inference framework addressing this question, by showing that valid statistical inference can be performed along the entire gradient descent trajectory, iteration by iteration, without requiring convexity of the loss landscape or convergence of the algorithm.

To illustrate this concept, we begin with a single-index (one-layer neural network) regression model and demonstrate how gradient descent iterates can be "debiased", at each iteration, to yield valid confidence intervals for the underlying signal and consistent estimates of generalization errors. We then extend this paradigm to the much more challenging setting of learning with general multi-layer neural networks in their full complexity, where the loss landscape can be arbitrarily complex. Crucially, the proposed method remains valid without requiring either algorithmic convergence or oracle knowledge of the unknowns, and may therefore inform practical decisions such as early stopping and hyperparameter tuning.

The key technical ingredient underlying this new inference paradigm is a recent entrywise dynamical mean-field theory for a broad class of first-order algorithms developed by the speaker.

This talk is partially based on joint with Masaaki Imaizumi and Xiaocong Xu.

I will discuss the natural gradient descent (NGD) optimization technique in the context of PDE-based optimization problems. In particular, I will discuss the impact of the norm used for generating the “information matrix” in connection with the regularity theory of the PDE under consideration. I will also discuss efficient computational methods for performing NGD using a streamlined least squares formulation and numerical linear algebra techniques.

Fluid models are often used to simplify large stochastic control problems. Instead of tracking stochastic transitions directly, they replace the transition dynamics with their expected evolution, resulting in a deterministic system that is much easier to analyze and optimize. While fluid models can reveal useful structural properties of good policies, they smooth out the randomness in the system and therefore lose information about the long-run state-action frequencies that are important for designing effective policies in the original stochastic problem.

We introduce a doubly randomized fluid model that keeps the simplicity of fluid approximations while preserving the original model’s key stochastic structure. Our starting point is a linear-program representation of long-run behavior, which directly captures how often each state-action pair is used under a policy. This formulation is typically intractable because it involves infinitely many variables and must satisfy a flow-balance condition for every state.

Our method addresses both difficulties using two layers of randomization. We approximate long-run state-action frequencies using a finite set of randomly generated basis functions, and we enforce flow balance through a finite number of randomly selected surrogate conditions. The resulting tractable problem recovers the long-run state-action frequencies that characterize good policies - something classical fluid models cannot do. We show that as we increase the number of basis functions and surrogate constraints, the approximation improves and converges to the behavior of the true stochastic system with epsilon-error guarantees.

Motivated by the challenges from nonconvex AC power flow physics and binary unit-commitment decisions, we investigate how convexification strategies—originally developed for AC Optimal Power Flow (ACOPF)—can be extended and unified within a full AC Unit Commitment framework. Our work compares several relaxations, including Second-Order Cone (SOC) relaxations, McCormick relaxations, and outer approximations of SOC constraints. To further strengthen these relaxations, we incorporate tightening techniques such as valid inequalities and a feasibility-based bound-tightening (FBBT) procedure that iteratively refines variable bounds. A key feature of our approach is that FBBT is performed on a relaxed UC model with slackened demand, allowing the resulting bounds to be reused across different demand profiles within a prescribed tolerance. Together, these components provide a more scalable convex relaxations for Unit Commitment with AC Power Flow.

Solving mixed-integer nonlinear programs (MINLPs) typically relies on constructing relaxations that are easier to tackle than the original problem. Recently, global parabolic (PARA) relaxations were introduced, featuring separable quadratic functions — paraboloids — as global under- or overestimators of general nonlinear constraint functions. At its core, this approach computes PARA approximations a-priori for a given tolerance and function domain. The resulting paraboloids are used in the form of a lookup table to construct relaxations by replacing the original constraint.

Although this approach provides great flexibility in approximating multivariate functions, its computational tractability is limited, preventing a direct comparison with established methods. We therefore propose a new method to compute PARA approximations a-priori that is bound to univariate constraint functions but achieves orders-of-magnitude speedups. This improvement enables a computational comparison with piecewise linear (PWL) relaxations on general MINLP instances from the well-established MINLPLib.

In particular, relaxations with varying tolerance are constructed for either technique and then solved. The results show that PWL relaxations are faster to solve when the tolerance is high, whereas for tighter tolerances the advantage shifts toward the PARA approach. This effect is attributed to the large number of auxiliary variables introduced by PWL to identify the active piece, which grows with smaller tolerance, while PARA relaxations do not require additional variables.

We consider integer programming problems with bounded general-integer variables belonging to the general class of network flow problems. For those, we computationally investigate the effect on mixed-integer linear programming (MIP) solvers of the different ways of producing extended formulations that replace a bounded general integer variable by a linear combination of a set of auxiliary binary variables linked by additional linear constraints. We show that MILP solvers perform very differently depending on which extended formulations is used and we interpret that different performance through the lens of cutting planes generation. Finally, we discuss a simple family of mixed-integer rounding inequalities that especially benefit from the reformulation, and we show its benefit within different MIP solvers. This provides methodological and practical guidelines for the use of those extended formulations in MIP and, to the best of our knowledge, this is the first extensive computational analysis of the topic.

This presentation will start with an overview of the Google OR-Tools' CP-SAT discrete optimization solver. Classical discrete optimization approaches, such as Mixed Integer Programming (MIP), SAT formulations, or Constraint Programming (CP), are often presented as competing techniques from which a practitioner must choose depending on the specific problem solved. CP-SAT achieves state-of-the-art performance by combining these diverse methodologies into a hybrid solver. Crucially, we combine the different approaches so they reinforce each other, resulting in a stronger solver than any single approach alone. In the second part of the talk, I will focus on the latest developments and new heuristics. Finally, I will discuss, from a practical point of view, how a hybrid solver architecture can help or hinder the integration of new heuristics.

Optimization problems that involve topology optimization in scenarios with large scale outages, such as post-disaster restoration or public safety power shutoff planning, are very challenging to solve. Using simple power flow representations such as DC power flow or network flow models results in low quality solutions which requires significantly higher-than-predicted load shed to become AC feasible. Recent work has shown that formulations based on the Second Order Cone (SOC) power flow formulation find very high-quality solutions with low load shed, but the computational burden of these formulations remains a significant challenge. With the aim of reducing computational time while maintaining high solution quality, this work explores formulations which replace the conic constraints with a small number of linear cuts. The goal of this approach is not to find an exact power flow solution, but rather to identify good binary decisions, where the power flow can be resolved after the binary variables are fixed. We find that a simple reformulation of the Second Order Cone Optimal Power Shutoff problem can greatly improve the solution speed, but that a full linearization of the SOC voltage cone equation results in an overestimation of the amount of power that can be delivered to loads.

In this work, we address the challenge of designing an algorithm that achieves the optimal convergence rate for optimizing unknown piecewise smooth (PWS) functions. PWS optimization refers to a special case of nonsmooth optimization in which the domain is partitioned into multiple regions, and the objective function is smooth within each region. The "unknown" setting means we have no prior knowledge of the locations or structure of these smooth regions; we treat the problem entirely as a black box.

Our approach builds upon the bundle level method. We rigorously show that when the number of cutting planes (cuts) exceeds the number of smooth pieces in the PWS function, the algorithm attains a convergence rate of $k\sqrt{L/\mu^*}\log(1/\varepsilon)$, where $k$ is the number of pieces, $L$ is the Lipschitz constant, $\mu^*$ is the quadratic growth (QG) modulus, and $\varepsilon$ is the target accuracy. Furthermore, we establish that this rate matches the corresponding lower complexity bound, thereby proving the optimality of our method. Notably, our algorithm is nearly parameter-free: apart from setting the number of cuts greater than the number of pieces, no additional problem-specific information is required.

Computed tomography (CT) is an imaging modality with widespread use in biomedical imaging, non-destructive materials testing, archaeology, etc. Standard approaches for inverting CT measurements, which rely on a numerically poorly conditioned preprocessing step, often struggle to reconstruct signals with high dynamic range (e.g., X-ray imaging of tissue with embedded metal). In this talk, I will present an iterative method based on nonsmooth optimization that sidesteps the aforementioned preprocessing step to obtain faster and more accurate reconstructions, and enjoys a linear rate of convergence under certain statistical assumptions.

Joint work with R. Willett.

Convexification is a classical technique in global optimization, but standard convex relaxations are typically one-shot: they provide a single convex surrogate that cannot be systematically tightened. As a result, convexification plays a limited role in nonlinear programming, where progressively improving approximations are essential in the latter area. In this talk, we introduce a parametric relaxation framework for a broad class of nonconvex optimization problems. The relaxation is governed by a parameter \gamma that explicitly trades off between the convexity of the surrogate model and its tightness to the original objective. A key feature of this framework is its hierarchical approximation property: as \gamma increases, the relaxed problem first touches the set of global minimizers, then all local minimizers, and eventually all stationary points of the original problem. We show that this metric provides a principled way to assess the tightness of stationary points and to guide algorithmic strategies for approaching globally optimal solutions. Algorithms and convergence rates to sequentially solve the \gamma-relaxation will also be discussed.

The Muon optimizer, a matrix-structured algorithm that leverages spectral orthogonalization of
gradients, is a milestone in the pretraining of large language models. However, the underlying mechanisms of Muon—particularly the role of gradient orthogonalization—remain poorly understood, with very few works providing end-to-end analyses that rigorously explain its advantages in concrete applications. We take a step by studying the effectiveness of a simplified variant of Muon through two case studies: matrix factorization, and in-context training of linear transformers. For both problems, we prove that Muon converges linearly with iteration complexities independent of the relevant condition number, provably outperforming gradient descent and Adam. Our analysis reveals that the Muon dynamics decouple into a
collection of independent scalar sequences in the spectral domain, each exhibiting similar convergence behavior. Our theory formalizes the preconditioning effect induced by spectral orthogonalization, offering insight into Muon’s effectiveness in these matrix optimization problems and potentially beyond.

Most existing rate and complexity guarantees for stochastic gradient methods in L-smooth settings impose restrictive requirements on steplength sequences, mandating that such sequences be non-adaptive, non-increasing, and abide by an upper bound given by 1/L. This severely limits the applicability and adversely affects the performance of such schemes for three key reasons: (i) requires knowing L; (ii) precludes larger steps when L is large; and (iii) employs non-increasing steplength sequences. Motivated by these shortcomings, we present an Armijo-enabled stochastic linesearch framework with standard stochastic zeroth- and first-order oracles. The resulting steplength sequence is non-monotonic. We provide the iteration and sample complexity of the proposed method and the sample complexity matches the standard SGD method. Preliminary numerical experiments are seen promising as the scheme is observed to compete well with de facto choices such as SGD and Adam across a range of problem types and instances.

We present efficient methods for solving monotone variational inequality problems when the operator is not uniformly Lipschitz continuous, and the feasible set may be unbounded. In recent years, optimization problems under generalized Lipschitz conditions have gained traction. However, there is limited literature on generalized conditions for VI problems. We introduce a new type of generalization of Lipschitz conditions for the VI problem that includes convex function-constrained optimization and variational inequality problems, average reward MDPs, and general robust optimization with semi-infinite (or robust) function constraints. In our problem setup, a uniform Lipschitz smoothness constant may not hold, especially since we allow the feasible set to be unbounded. Under the assumption that a VI solution exists and after carefully designing the algorithm, we show that our iterates remain bounded and establish state-of-the-art convergence rate guarantees for various cases, including monotone or strongly monotone operators with possibly generalized nonsmoothness. When the operator is stochastic, we show expected convergence guarantees. We also establish that the iterates are bounded in expectation. Our convergence criterion uses a Nesterov-type sufficiently large radius ball that contains the VI solution. However, unlike Nesterov, we do not need to know this radius to implement our methods for the abovementioned cases.

In this talk, we introduce a novel, first-order method for solving density-based multi-material topology optimization problems. A fundamental challenge in multi-material design is ensuring that the material indicator field satisfies complex local constraints—specifically, that material fractions are non-negative and sum up to one. For arbitrary material sets, these constraints define a high-dimensional convex polytope, making traditional projection-based algorithms computationally cumbersome and ill-suited for general interpolation domains.Our proposed method is a major extension of the single-material SiMPL algorithm. The key innovation is the generalization of the Bregman divergence -- a concept of generalized distance -- to implicitly enforce the constraints of any arbitrary polytope. This is achieved by deriving a new mirror map that establishes a bijection between the physical material distribution and an unconstrained latent variable. This map allows the optimization process to be performed entirely in the unconstrained latent space, leading to a simple, explicit update rule based only on the objective function's first-order derivative. Because we discretize this unconstrained latent field, and the mirror map guarantees feasibility, the method robustly produces a sequence of pointwise-feasible iterates regardless of the finite element order used, guaranteeing that all material mixtures respect the design space throughout the optimization. We will demonstrate the method's superior performance and robustness across several benchmark problems, showing rapid convergence and high-contrast, manufacturable designs.

We present a stochastic algorithm for Conditional Value-at-Risk (CVaR) optimization that tightly couples adaptive importance sampling with a Bregman proximal-point scheme. Starting from the dual representation of CVaR as a supremum of expectations over a risk envelope, the proposed algorithm alternates between (i) solving a stochastic primal subproblem and (ii) updating dual probabilities via a proximal step. The evolving dual distribution concentrates sampling mass in the tail that determines CVaR, dramatically improving efficiency at high confidence levels where naïve sampling is wasteful. A key design choice is the introduction of a Legendre-type function based on a generalized Fermi–Dirac entropy, which induces a Bregman divergence respecting the box constraints of the dual feasible set. This yields closed-form, logistic-style updates that keep iterates in the interior and provide a built-in, data-driven importance sampler whose likelihood ratio converges to the risk identifier of the optimal solution. Consequently, the primal stages can be solved with few informative samples while the dual stages steer sampling toward the active tail. We prove convergence in function values for general convex objectives and almost sure convergence of iterates under mild assumptions, and we derive rates for smooth and strongly convex subproblems. Numerical experiments demonstrate strong performance across representative tasks in financial portfolio selection, support vector regression, and risk-averse large-scale topology optimization, where the algorithm outperforms standard stochastic approximation baselines, particularly as the CVaR level approaches one. Overall, the work provides a practical approach to scalable CVaR optimization via entropy-regularized proximal updates and automatically tuned tail sampling.

We study a vertex deletion problem motivated by network robustness: remove a minimum-cost set of vertices to reduce network density below a given threshold. This measures the sensitivity of the densest subgraph to non-uniform vertex deletions. We show this problem unifies several classical problems in combinatorial optimization and present a simple bicriteria approximation algorithm.

We consider the k-edge connected spanning subgraph (kECSS) problem, where we are given an undirected graph with nonnegative edge costs, and we seek a minimum-cost subgraph that is k-edge connected, i.e., there are k edge-disjoint paths between every pair of nodes.

For even k, we present a polytime algorithm that computes a (k-2)-edge connected subgraph of cost at most the optimal value of the natural LP-relaxation for kECSS; for odd k, this leads to a (k-3)-edge connected subgraph of cost at most the LP optimum. Since kECSS is APX-hard for all k≥2, the tightest connectivity guarantee one can hope for without exceeding the optimal value is (k-1)-edge connectivity, so our results are nearly tight. They also significantly improve upon the recent work of Hershkowitz et al., both in terms of solution quality and the simplicity of algorithm and its analysis. Our result is obtained using iterative rounding, with additional insights involving uncrossing tight sets for kECSS. A slight variant of our approach yields (k-1)-edge connectivity incurring cost at most 1.5 times the LP optimum; with unit edge costs, the cost guarantee improves to (1+4/3k) times the LP optimum, which improves upon the state-of-the-art for unit edge costs, but with a unit loss in edge connectivity. The kECSS-result also yields results for the k-edge connected spanning multigraph (kECSM) problem, where multiple copies of an edge can be selected. Our techniques extend to the degree-bounded versions of kECSS and kECSM.

Joint work with Nikhil Kumar.

We present CANOPI, a novel algorithmic framework, for solving the Contingency-Aware Nodal Power Investments problem, a large-scale nonlinear optimization problem that jointly optimizes generation, storage, and transmission expansion. The underlying problem is nonlinear due to the impact of transmission upgrades on impedances, and the problem's large scale arises from the confluence of spatial and temporal resolutions. We propose algorithmic approaches to address these computational challenges. We pose a linear approximation of the overall nonlinear model, and develop a fixed-point algorithm to adjust for the nonlinear impedance feedback effect. We solve the large-scale linear expansion model with a specialized level-bundle method leveraging a novel interleaved approach to contingency constraint generation. We introduce a minimal cycle basis algorithm that improves the numerical sparsity of cycle-based DC power flow formulations, accelerating solve times for the operational subproblems. CANOPI is demonstrated on a 1493-bus Western Interconnection test system built from realistic-geography network data, with hourly operations spanning 52 week-long scenarios and a total possible set of 20 billion individual transmission contingency constraints. Numerical results quantify the reliability and economic benefits of fully incorporating transmission contingencies in integrated planning models and highlight the computational advantages of the proposed methods.

A variety of challenges face the U.S. electricity system: backlogged generator interconnection queues, a rapidly changing resource mix, increasingly impactful extreme weather events, and the prospect of increasing electricity demand after decades of little growth. Transmission expansion has the potential to address many of these issues, but traditional approaches to transmission planning—which have historically been regional and reactive in nature—do not capture the many benefits that transmission can provide. Here, we describe an integrated optimization approach to facilitate wide-area, multi-value transmission planning. The core of the approach uses two coupled models: a capacity expansion model (CEM) that cooptimizes the capacity and operations of generation, storage, and transmission at zonal resolution across the U.S., and a resource adequacy model that assesses the reliability of the resulting systems at high time resolution with Monte Carlo sampling of individual generator outages, iteratively identifying “stress periods” to guide capacity investments in the CEM. Three additional optimization approaches feed into the CEM: The selection of representative time periods is optimized to minimize regional errors in electricity demand and solar/wind capacity factors; zonal interface transfer limits are determined via a nodal power flow model that that optimizes interregional transfers; and a “transmission upgrade supply curve” is determined for each zonal interface, wherein congested transmission branches are sequentially upgraded (via reconductoring, rebuilds, voltage upgrades, and new parallel circuits/transformers) to minimize the cost per megawatt of interregional transfer capability. This integrated planning approach identifies generation/transmission portfolios that minimize system costs while meeting reliability targets nationwide.

Traditional power grid expansion planning studies assume that load growth is an exogenous parameter of the optimization model. The proliferation of datacenters supporting the AI revolution, as well as industrial facilities built for direct air capture of CO2 are challenging this view. These loads can be sited with some flexibility and are large enough to have an impact on the electricity network, even at the transmission level. In this work, we incorporate large-scale loads with flexible siting as part of the decision variables in a nodal capacity expansion planning model. The operational flexibility of large loads is modeled through tranches with differentiated reliability and incorporated into the model as expectation constraints. We develop an augmented Progressive Hedging algorithm that uncouples expectation constraints to solve this problem exploiting parallel computing. We show through some numerical testcases on realistic systems based on California and South Carolina the impact of proactive load siting on the need for transmission investments with different power flow formulations to highlight the importance of incorporating them into software planning tools.

Electric power infrastructure faces increasing risks of damage and disruption due to wildfires. The planning and operation of electric grids in wildfire-prone regions must consider the potential impact of unpredictable fires. However, traditional wildfire models do not easily answer how and when wildfires may disable some set of power grid elements. To address these questions, we propose an optimization model to characterize an adversarial wildfire that targets specific infrastructure while respecting physical fire spread dynamics. Using the Rothermel fire spread model, we design a feasible fire spread region that varies with wind velocity. We develop principled convex relaxations of the spread region, introducing new relaxations of the inner product over Euclidean balls. We formulate and solve a mixed-integer conic program to characterize the worst-case time to outage and load shed of a realistic wildfire, and present test cases derived from the recent Park, Eaton, and Palisades fires in California.

Nonconvexities in markets with discrete decisions and nonlinear constraints make efficient pricing challenging, often necessitating subsidies. A prime example is the unit commitment (UC) problem in electricity markets, where costly subsidies are commonly required. We propose a new pricing scheme for nonconvex markets with both discreteness and nonlinearity, by convexifying nonconvex structures through a semidefinite programming (SDP) relaxation and deriving prices from the relaxation's dual variables. When the choice set is bounded, we establish strong duality for the SDP, which allows us to extend the envelope theorem to the value function of the relaxation. This extension yields a marginal price signal for demand, which we use as our pricing mechanism. We demonstrate that under certain conditions - for instance, when the relaxation's right hand sides are linear in demand - the resulting lost opportunity cost is bounded by the relaxation's optimality gap. This result highlights the importance of achieving tight relaxations. The proposed framework applies to nonconvex electricity market problems, including for both direct current and alternating current UC. Our numerical experiments indicate that the SDP relaxations are often tight, reinforcing the effectiveness of the proposed pricing scheme. Across a suite of IEEE benchmark instances, the lost opportunity cost under our pricing scheme is, on average, 46% lower than that of the commonly used fixed-binary pricing scheme.

We study the Capacitated Multi-echelon Batch Inventory problem, which models batch inventory ordering in a two-echelon fulfillment network with capacity constraints at the lower echelon, as often seen in e-retailing. A centralized upper echelon places batch orders from a supplier and allocates units to the lower echelon to minimize total costs, including ordering, allocation, and expected fulfillment costs.

The hard capacity constraints make the problem structurally similar to capacitated facility location, which is well known to be challenging to approximate using linear relaxations. We address these challenges by exploiting key features of our setting: the allocation and fulfillment costs are metric, and the upper echelon is uncapacitated.

Leveraging these features, we develop a new clustering technique that uses information from the dual program to group nearby distribution centers in the lower echelon with enough inventory. This clustering enables us to round the LP solution while satisfying the capacity constraints and leads to a constant-factor approximation algorithm for the offline two-stage version of the problem.

We further extend our approach to an online setting in which demand arrives sequentially over a finite time horizon. To address this, we develop a framework that integrates our approximation algorithm with tools from the online decision-making literature to decide on the ordering, allocation, and online fulfillment. Our framework achieves an additive regret bound that is constant relative to the time horizon, making the theoretical guarantee of the online solution converge linearly to that of the offline solution as the horizon grows.

Recent advances in neural combinatorial optimization have shown that deep reinforcement learning can solve combinatorial optimization problems by training policies to generate near-optimal solutions. However, existing approaches struggle with constraints that require global reasoning on the entire decision space, such as time window constraints augmenting the traveling salesperson problem. Most existing approaches only handle simpler constraints, such as permutation structures, where one can ensure feasibility incrementally while generating the solution.

We address this limitation by integrating relaxed decision diagrams into the reinforcement learning framework. Our approach equips the learner with structured lookahead capabilities that enable it to anticipate and avoid decisions leading to downstream infeasibility during solution generation. We embed these decision diagrams directly into attention-based neural architectures and evaluate their effectiveness on the traveling salesperson with time windows problem and precedence constraints.

This paper develops a unified high-order accumulative regularization (AR) framework for convex and uniformly convex gradient norm minimization. Existing high-order methods often exhibit a gap: the function-value residual decreases fast, while the gradient norm converges much slower. To close this gap, we introduce AR that systematically transforms the fast function-value residual convergence rate into a fast (matching) gradient norm convergence rate.

Specifically, for composite convex problems, to compute an approximate solution such that the norm of its (sub)gradient does not exceed $\varepsilon,$ the proposed AR methods match the best corresponding convergence rate for the function-value residual. We further extend the framework to uniformly convex settings, establishing linear, superlinear, and sublinear convergence of the gradient norm under different lower curvature conditions. Moreover, we design parameter-free algorithms that require no input of problem parameters, e.g., the Lipschitz constant of the $p$-th-order gradient, the initial optimality gap and the uniform convexity parameter, and allow an inexact solution for each high-order step. To the best of our knowledge, no parameter-free methods can attain such a fast gradient norm convergence rate which matches that of the function-value residual in the convex case, and no such parameter-free methods for uniformly convex problems exist.These results substantially generalize existing parameter-free and inexact high-order methods and recover first-order algorithms as special cases, providing a unified approach for fast gradient minimization across a broad range of smoothness and curvature regimes.

First-order methods excel at delivering the medium-to-low accuracy solutions required for large-scale optimization. However, it is precisely in this regime that standard termination certificates, like the gradient norm or KKT conditions, prove unreliable. Their failure to measure true algorithmic progress leads to premature termination or wasted computation, and critically, prevents their use in designing adaptive algorithms.

We introduce a novel, verifiable termination certificate that resolves this shortcoming. Our certificate enables the use of the notoriously difficult-to-estimate quadratic growth (or weak convexity) parameter to develop parameter-free algorithms with nearly optimal complexity bounds. We demonstrate its power across a wide range of problems, including non-smooth, weakly smooth, piecewise smooth, and constrained optimization.

Many decision problems require acting after observing covariates, while predictions arrive from multiple models of varying reliability. We propose a trust-weighted ensemble that lifts stacking to the decision layer via an empirical residual-based sample average approximation (ER-SAA) method. Given T bases and a simplex weight w, each base contributes a residual pool from training data. At a new context x we form scenarios by pairing current prediction of each base and residuals from its own pool and solve a trust-weighted ER-SAA by allocating probability mass across bases according to wWe operate online in a “decide – observe – learn” protocol that updates w by either a prediction-driven or a decision-aware strategy. The framework adapts trust to source reliability and converts multi-source reference information into decisions that respect both central tendency and structured uncertainty.

Many optimization problems solved repeatedly in practice—such as in power systems or logistics—share a fixed structure but with different parameter values. Re-optimizing each instance from scratch can be computationally prohibitive, while adopting a single solution for all realizations is overly conservative. The K-adaptability paradigm offers a compromise by precomputing a small set of candidate solutions and selecting the best once uncertainty is revealed. When K is small, this approach enables real-time decision making and provides human decision-makers with a manageable set of implementable solutions. However, it is often impossible to identify a small number of representative solutions that remain feasible across all parameter realizations. We introduce K-adaptability with recourse, in which the goal is to precompute K candidate assignments for the discrete variables while allowing the continuous recourse decisions to be optimized online after uncertainty is observed. This structure preserves online computational tractability while substantially improving adaptability. We present formulations of this model and two additional solution approaches: a scenario-based heuristic inspired by scenario decomposition, and an exact column and constraint generation method. Preliminary results show that the heuristic consistently identifies solutions close to the optimum and demonstrate the computational efficiency of the column and constraint generation approach.

By modeling iterative algorithms as expression trees, we propose an MINLP model for automated discovery of optimization algorithms. The formulation seeks functions that satisfy a given convergence criterion while a desired objective, such as the complexity of the function, is optimized. We solve the MINLP with BARON and present results for the design of first-order algorithms and systems of equations.

This work introduces several novel mixed-integer linear programming (MILP) unit commitment models for thermal generators (e.g. bioenergy, natural gas, coal, etc.) featuring the ability to model adaptable power output trajectories during multi-period startups and shutdowns. These models capture various degrees of fidelity to fine-resolution dynamic models including variable power output and variable transition duration along with accurate representations of transition costs. We demonstrate that these models yield notably improved objective function values over existing state-of-the-art models while maintaining adequate computational performance. While the body of this work is performed within the context of electrical power grid operation, we demonstrate that these models are easily generalized to other sectors that involve nuanced state transitions such as chemical production scheduling and optimal control.

We study the effectiveness of cutting-plane (CP) methods relative to branch-and-cut (BC) techniques for convex 0/1 integer programming under variable disjunctions. A recent work of Basu et al. shows that for every valid inequality an optimal CP proof is at least as efficient as an BC proof, but possibly only up to an arbitrarily small epsilon slack in the bound. This raises a further question: can one obtain zero slack, at least in some structured settings, while keeping the proof size under control?

We give a positive answer in three dimensions. For convex 0/1 sets C contained in [0,1]^3 with variable disjunctions, we prove that every valid inequality that admits a BC proof of size N also admits a pure CP proof of size at most N, with zero slack. Thus, in 3D, branching does not provide additional proof power beyond what can be achieved with cutting planes alone.

The proof proceeds by eliminating maximal-depth branching nodes via structural induction, together with a geometric analysis of five types of interactions between the cube and the cutting planes. I will describe this construction, place the result within the Basu-et-al. landscape, and discuss obstacles and possible approaches to extending such zero-slack CP-BC efficiency guarantees beyond three dimensions.

We study convex optimization problems where disjoint blocks of variables are controlled by binary indicator variables, which themselves are subject to additional conditions such as cardinality constraints. Several classes of important examples can be formulated in such a way that both the objective and the constraints are block-separable convex functions. We introduce a new metric for sparsity, namely, the treewidth of the constraint-block intersection graph. We provide complexity results and a family of polynomial-time approximation algorithms for those types of problems.

We study alternative continuous cubic formulations for a class of fundamental cluster-detection problems in networks, including the maximum clique, s-defective clique, s-plex, and the independent union of cliques problems. While our broader work pertains to each of these problems, this talk will focus in particular on the maximum 2-plex problem. These combinatorial problems, central to graph theory and network analysis, have wide-ranging applications in social network analysis, bioinformatics, and communication systems. The considered formulations capture the structural properties of these discrete problems while enabling their treatment within a continuous optimization framework. Each problem is modeled as the maximization of a degree-three multi-linear polynomial objective over a polyhedral feasible set defined by either the standard simplex or unit hypercube.

These continuous models provide new perspectives for theoretical analysis and algorithmic development. In particular, we derive analytical bounds that offer insight into the quality of feasible solutions and can be used to guide solution refinement. Additionally, we leverage the structure of the cubic formulations to design local and global solution methods based on existing nonlinear optimization solvers.

To assess the practical utility of our approach, we conduct a comprehensive computational study on standard benchmark instances. We compare the performance of multiple solvers in terms of solution quality, computational efficiency, and scalability. Our numerical results demonstrate the effectiveness of the proposed formulations in yielding high-quality solutions and suggest that continuous polynomial models can serve as a powerful alternative to traditional discrete or integer programming techniques for cluster detection in complex networks.

In this talk, we propose a class of implementable, optimal-rate methods for nonconvex optimization over Riemannian manifolds when derivative oracles up to order p are available. Our analysis proves that, under this regime, Riemannian optimization is no harder than Euclidean optimization: the Riemannian convergence rates match the optimal Euclidean rates. Furthermore, we comprehensively characterize how the regularity of the retraction and objective function jointly impacts that of the pullback functions. The engine driving this characterization is a novel Faà di Bruno–type formula for pullback derivatives, constructed using the rarely applied higher-order covariant differential calculus. Moreover, we show that the entire class of efficient retractions of Gawlik et al. can be applied within these methods.

Our theoretical framework yields implementable third-order methods when equipped with suitable subproblem solvers. The mechanism enabling implementation is a symmetric Krylov tensor subspace method. Building on the Savas–Eldén approach, it can hybridize with recent techniques for minimizing quartically regularized third-order polynomials.

This work studies sample average approximation (SAA) in solving convex or strongly convex stochastic programming (SP) problems. In estimating SAA's sample efficiency, existing bounds often entail metric entropy terms, which often grow polynomially with problem dimensionality. Although metric entropy-free complexity results are known, they based on uniform Lipschitz conditions, which can be overly critical in many settings. Our work presents metric entropy-free sample complexity results for the SAA without uniform Lipschitz condition. The new results often lead to an $O(d)$-improvement in the complexity rate than the state-of-the-art, which also shows the matching sample complexity of SAA and SMD (stochastic mirror descent). Our numerical experiment results on SAA for solving a simulated SP problem align with our theoretical findings.

Constrained sampling --- the problem of sampling from a density corresponding to a potential function that is supported on a constrained subset of a normed space --- arises prominently in Bayesian inference, uncertainty quantification, statistical estimation, and modern machine learning. The derivative-free version of this problem stipulates that we have access only to a zero-order oracle of the potential function. A leading recent approach, Mirror MCMC, operates by (i) mapping points from the constrained (primal) space to an unconstrained dual space via the gradient of a mirror map, (ii) applying a standard unconstrained MCMC update in the dual space, and finally (iii) mapping back to the primal space via the inverse mirror map. In this talk, we develop the derivative-free version of Mirror MCMC. After showing that the Mirror MCMC iterates admit an equivalent primal representation where each update corresponds to the solution of a stochastic mirror-optimization subproblem, we construct an implementable inexact sequence using the zero-order oracle. We then establish sufficient conditions under which the resulting derivative-free sampler achieves condition-number–free mixing rates, and characterize how much inexactness can be tolerated in the subproblem solves without degrading these rates. For furthering intuition, we will illustrate with examples.

Wildfires pose a significant threat to the reliability and security of power distribution systems where wildfire-induced power line failures exhibit spatial correlations and are influenced by operational decisions. To address this challenge, we introduce a new joint chance-constrained stochastic programming model that accounts for both spatial dependency and decision-dependent uncertainty. Power line failures are affected by wildfires and their propagation (exogenous uncertainty) and by power flow decisions (decision-dependent uncertainty). The model simultaneously reconfigures the network topology in the face of wildfire ignition risks and minimizes the total costs associated with power generation, switching actions, and penalties for power imbalance while ensuring the operational reliability of the network. We incorporate an Archimedean copula to capture spatial dependency across power line failures, and introduce a decision-dependent distortion function to model how power flows impact line failure probabilities. We develop a scenario generation method providing scenarios representative of how wildfires affect power line failures and derive a deterministic and equivalent mixed-integer nonlinear programming reformulation of the chance-constrained model. Additionally, we approximate the truncated generator of the parameterized copula using a piecewise linear function in order to obtain a tractable mixed-integer linear programming formulation. The numerical tests showcase the applicability and effectiveness of the proposed reformulation and algorithmic approach, which outperforms alternative modeling approaches by (i) reducing the out-of sample loss of loads, (ii) increasing the reliability of the network, (iii) lowering the loading percentage of power lines in the most vulnerable areas, and (iv) decreasing the probability of electrically-induced wildfires.

We present a methodology to handle decision dependency via a residuals-based distributionally robust optimization model, where the underlying uncertainty depends on both covariate information and the decisions. We adopt regression models to learn the latent decision dependency and construct a nominal distribution (thereby ambiguity sets) around the learned model using empirical residuals from the regressions. Ambiguity sets can be formed via the Wasserstein distance, a sample robust approach, or with the same support as the nominal empirical distribution (e.g., phi-divergences), where both the nominal distribution and the radii of the ambiguity sets could be decision- and covariate-dependent. We establish conditions under which desired statistical properties, such as asymptotic optimality, rates of convergence, and finite sample guarantees, are satisfied. We present decomposition-based solution methods to efficiently solve the resulting (typically nonconvex) models under parametric and nonparametric regression models and devise data driven cross-validation schemes to size the ambiguity sets. Through numerical experiments, we illustrate the effectiveness of our approach and the benefits of integrating decision dependency into a residuals-based contextual DRO framework.

Given positive integers $m,n$, a fractional assignment $x\in [0,1{]}^{m\times n}$ and weights $d\in {\mathbb{R}}_{>0}^n$, we show that there exists an assignment $y\in \{0,1{\}}^{m\times n}$ so that for every $i\in [m]$ and $t\in [n]$,

$$\mid \sum _{j\in [t]}{d}_j(x_{ij}-y_{ij})\mid <\underset{j\in [n]}{\max }{d}_j\mathrm{.}$$

This generalizes a result of Tijdeman (1973) on the unweighted version, known as the chairman assignment problem. This also confirms a special case of the single-source unsplittable flow conjecture with arc-wise lower and upper bounds due to Morell and Skutella (IPCO 2020). As an application, we consider a scheduling problem where jobs have release times and machines have closing times, and a job can only be scheduled on a machine if it is released before the machine closes. We give a $3$3-approximation algorithm for maximum flow-time minimization.

Convex relaxations are of central interest in optimization, and it is typically challenging to determine whether a given convex relaxation will be tight for a given problem. We introduce a topological framework for analyzing situations in which a constrained optimization problem over a nonconvex set (such as a manifold) has a tight convex relaxation. In particular, we give a criterion for the existence of such a tight convex relaxation in terms of the existence of a continuous function of Lagrange multipliers for the constrained problem maximizing the corresponding Lagrangian. We call such a function a Lagrangian dual section, in reference to the topological notion of a section of a bundle.

As a corollary of this result, we will give new criteria for the exactness of SDP relaxations for Stiefel manifold optimization and inverse eigenvalue problems in terms of linear subspaces of matrices satisfying spectral properties such as being nonsingular. We will also illustrate a homotopy continuation style algorithm with global optimality guarantees with applications to the unbalanced procrustes problem.

Joint work with Venkat Chandrasekaran, Jose Israel Rodriguez, Tim Duff.

We consider the Electric Vehicle Routing Problem, where customers are served by an electrified fleet that recharges at charging stations. The set-partitioning formulation for the EVRP has infinitely many path-based variables due to continuous charging decisions. We develop a column generation algorithm, where the pricing problem is decomposed into two phases: (i) generating subpaths between two consecutive charging stations, and (ii) combining subpaths into paths. We formalize subpath-based domination properties to establish the finite convergence and exactness of the column generation algorithm. We design domination criteria and prove that the resulting two-phase column generation algorithm can solve several EVRP variants, including time windows, nonlinear charging functions, heterogenous charging costs, and relaxation tightening strategies (e.g., ng-routes, limited-memory subset-row inequalities). Real-world experiments show that the two-phase algorithm outperforms path-based label-setting benchmarks, and can scale to large problem instances.

Empirical risk minimization (ERM) is a cornerstone of modern machine learning (ML), supported by advances in optimization theory that ensure efficient solutions with provable algorithmic convergence rates and statistical learning rates. Privacy, memory, computational, and communications constraints increasingly necessitate data collection, processing, and storage across network-connected devices. In many applications, these networks operate in decentralized settings where a central server cannot be assumed, requiring decentralized ML algorithms that are both efficient and resilient. Decentralized learning, however, faces significant challenges, including an increased attack surface for adversarial interference during decentralized learning processes. This paper focuses on the man-in-the-middle (MITM) attack, wherein adversaries exploit communication vulnerabilities between devices, potentially causing models to deviate significantly from their intended ERM solutions. To address this challenge, we propose RESIST (Resilient dEcentralized learning using conSensus gradIent deScenT), an optimization algorithm designed resilient to compromised communication links. Unlike existing adversarially robust decentralized learning methods, which often (i) guarantee convergence only to a neighborhood of the solution, (ii) lack guarantees of linear convergence for strongly convex problems, or (iii) fail to ensure statistical consistency as sample sizes grow, RESIST overcomes all three limitations. It achieves algorithmic and statistical convergence for strongly convex, Polyak–Łojasiewicz, and nonconvex ERM problems by employing a multistep consensus gradient descent framework and robust statistics-based screening methods to mitigate the impact of MITM attacks. Experimental results demonstrate the robustness and scalability of RESIST across diverse attack strategies, screening methods, and loss functions, confirming its suitability for real-world decentralized optimization and learning in adversarial environments.

In computational science workflows, it is often the case that 1) simulations do not permit readily the algorithmic computation of derivatives, 2) objective functions for optimization involve multiple simulation outputs, and 3) those simulation outputs can be performed (at least partially) in parallel. While the Argonne-maintained derivative-free optimization (DFO) solvers in the Interpolation-Based Composite Derivative-Free Optimization (IBCDFO) repository are specifically designed to address problems exhibiting the first two of these characteristics, they do not support the third at the time of writing this abstract.

In this work, we reexamine our past work on a class of randomized algorithms, stochastic average model (SAM) methods. SAM methods are conceptually similar to stochastic average gradient methods, and effectively require that only randomized subsets of simulation outputs be locally modeled in each iteration of a model-based DFO method. From the perspective of variance reduction, given a batch size, the optimal sampling distribution for a batch within a SAM method permits a closed-form, but the closed-form distribution contains unknowable parameters.

However, clues concerning the values of the unknowable parameters are often available either through domain expertise, alternative lower-fidelity simulations, or simply regression models. Motivated by this, we consider the problem of determining sampling distributions for SAM methods as a contextual bandit problem and we apply the Exponential weights algorithm for Exploration and Exploitation with Experts (Exp4). Results are demonstrated on both synthetic problems and a problem of density functional theory calibration.

: I will explore a number of recent theory-driven advances in our understanding of learning rates and learning rate schedules, and show how they can be applied to training modern large scale language models.

Distributionally Robust Optimization (DRO) is a widely adopted framework for decision-making under uncertainty, but it may lead to overly conservative solutions. Ex-ante Distributionally Robust Regret Optimization (DRRO) has been raised to overcome this over-pessimism, where the goal is to minimize the gap between the ex-ante optimal solution and the decision. While even structured problems can lead to NP-hard formulations in Wasserstein DRRO, we identify formulations that are tractable, applicable in important practical problems and that lead to solutions that can be aggressive in the presence of uncertainty. For example, we examine several examples motivated by applications in bipartite matching, shortest-path problems, and other fields such as Reinforcement Learning from Human Feedback (RLHF).

Probing is an important technique for domain reduction in mixed-integer programming. This paper introduces a novel framework on constructing probing rules based on the sensitivity analysis of value functions for general mixed-binary convex optimization problems. In contrast to existing methods that mainly rely on computationally expensive reoptimization or on Lagrangian dual multipliers derived from a convex relaxation, our proposed method offers greater flexibility and achieves a balance between computational efficiency and probing strength. Computational results also demonstrate the strong performance of our new methods.

Schrage and Wolsey (1985) propose a sensitivity analysis for integer programs that provides dual bounds parameterized by right-hand side perturbations utilizing the structure of branch-and-bound trees. Taking a similar branch-and-bound tree-based perspective, we show how to develop dual bounds that consider simultaneous perturbations in right-hand sides, objective coefficients, and constraint matrices. We then demonstrate how these insights enable effective warm starts for sequences of perturbed mixed-integer programs, relevant in time-sensitive settings where models must be re-solved as operational conditions evolve. We present numerical results for unit commitment problems, where frequent successive data updates and tight decision timelines make rapid re-optimization highly beneficial.

Quantities such as escape probabilities in Markov chains can be recovered by solving linear systems in diagonally dominant matrices arising from graph Laplacians. However, these probabilities may vary exponentially across different nodes of the graph. As a result, under standard norm-based error guarantees, accurately recovering small entries requires exponentially small error parameters, leading to prohibitively large running times. For example, a direct application of existing Laplacian solvers or fast matrix multiplication algorithms requires $\Omega(mn^2)$ and $\Omega(n^{\omega+1})$ bit operations, respectively, where $m$ denotes the number of nonzero entries, $n$ is the matrix dimension, and $\omega$ is the matrix multiplication exponent.

In this work, we study faster algorithms for solving such linear systems with entrywise approximation guarantees. In particular, we present an almost-linear-time algorithm for solving symmetric diagonally dominant M-matrix (SDDM) linear systems with entrywise approximation guarantees.

The iteratively reweighted least squares (IRLS) method is a widely used practical approach to solving regression problems. When the quality of fit is measured by the mean squared error, the problem reduces to ordinary least squares, and hence to solving a linear system. When other measures of error are used, such as the ℓ_∞ norm of the vector of prediction errors, the problem becomes substantially harder and, at high accuracy, is equivalent to linear programming. IRLS addresses these settings by solving a sequence of weighted least squares subproblems with adaptively chosen weights. However, many variants lack robust convergence guarantees, and existing analyses often fail to explain their empirical performance.

We present a unified class of IRLS algorithms for ℓₚ regression for all p ∈ [1, ∞]. Our methods match the best known theoretical guarantees, are simple to implement, and perform well in practice. We further extend this framework to overconstrained regression with quasi self-concordant penalties, which capture a broad family of loss functions, and obtain iteration bounds that depend only on the feature dimension rather than the number of data points.

Based on joint work with Alina Ene and Ta Duy Nguyen.

This work investigates Fréchet regression and related barycenter-type optimization problems in general metric and Riemannian settings, with particular emphasis on the sphere. Fréchet regression models the relationship between scalar or vector-valued predictors and responses in a metric space by defining the regression output as the minimizer of an affine combination of distance-based loss functions. The weights in this affine combination are determined by the data, while the functions are constructed from the underlying metric. Existing formulations predominantly specify such cost functions but provide limited guidance on how to compute the minimizers or analyze the associated optimization problems, especially in non-Euclidean geometries.

In this study, the Fréchet regression loss is formulated as an affine combination of convex functions. Depending on the structure of the weights, the resulting problem may remain convex or become a genuine difference-of-convex (DC) problem. In the convex regime, established tools from convex optimization on manifolds yield direct solution methods. In the nonconvex regime, including signed-weight barycenter problems on hemispheres, the loss is represented as a DC function and the Difference-of-Convex Algorithm (DCA) on manifolds is employed. A rigorous existence theory is developed, establishing conditions under which minimizers of the regression loss exist in the manifold setting. In addition, a detailed and solid convergence analysis of manifold DCA is provided given the objective function, extending Euclidean convergence guarantees and connecting them with real-analytic structure and Kurdyka–Łojasiewicz-type inequalities. Together, these results offer a comprehensive theoretical and algorithmic framework for Fréchet regression in general positive constant curvature metric spaces.

Barycenter problems encode important geometric information about a metric space. While these problems are typically studied with positive weight coefficients associated to each distance term, more general signed Wasserstein barycenter problems have recently drawn a great deal of interest. These mixed sign problems have appeared in statistical inference setting as a way to generalize least squares regression to measure valued outputs and have appeared in numerical methods to improve the accuracy of Wasserstein gradient flow solvers. Unfortunately, the presence of negatively weighted distance terms destroys the L^2 convexity of the unsigned problem, resulting in a much more challenging optimization task. In this talk, I will discuss some theoretical properties of these mixed sign barycenter problems, focusing on sufficient conditions to guarantee the global optimality and uniqueness of a critical point.

Joint work with Bohan Zhou.

Fréchet regression is a flexible framework for modeling relationships in which the response variables are complex non-Euclidean objects, such as networks or probability distributions. In this work, we focus on the problem of estimating the conditional Fréchet mean on the Bures-Wasserstein manifold, the space of positive-definite matrices equipped with the Wasserstein metric. This estimation task leads to an optimization problem closely related to the Wasserstein barycenter problem, but with the notable extension of allowing negative weights. We establish a sufficient condition that guarantees the existence of a minimizer for this problem. To compute the estimator, we develop a Riemannian gradient descent algorithm and prove that it converges at a sublinear rate. Finally, we validate our approach through empirical experiments in a network regression setting.

We propose a new decomposition framework for continuous nonlinear constrained two-stage optimization, where both first- and second-stage problems can be nonconvex. A smoothing technique based on an interior-point formulation renders the optimal solution of the second-stage problem differentiable with respect to the first-stage parameters. As a consequence, efficient off-the-shelf optimization packages can be utilized. We show that the solution of the nonconvex second-stage problem behaves locally like a differentiable function so that existing proofs can be applied for the global convergence of the first-stage. We also prove fast local convergence of the algorithm as the barrier parameter is driven to zero. Numerical experiments for large-scale instances demonstrate the computational advantages of the decomposition framework.

Bilevel optimization can model hierarchical optimization problems for a multitude of applications, including environmental policy [1]. Bilevel optimization problems are very challenging to solve - they are NP-hard even for the linear continuous case [2]. Common analytical approaches include KKT-based reformulations [3], branch-and-bound/cut, and multi-parametric reformulation [4]. These methods are either computationally inefficient and struggle with medium to large scale problems, or cannot be used for bilevel problems with lower-level discrete variables.

In contrast, heuristic approaches (e.g. DOMINO [1]) use a hybrid strategy: a deterministic solver handles the lower-level problem, which is then treated as a black box by the upper-level metaheuristic solver. These methods guarantee bilevel feasibility for some classes of problems, but offer no finite-time guarantees regarding convergence to either a global or even a local optimum. We propose a critical region-evaluation strategy, adapted from multi-parametric optimization, to enhance black box methods and improve their local search efficiency and solution quality.

We validate the method, comparing its performance against both deterministic multi-parametric reformulation and established global heuristic optimization algorithms. These results demonstrate that combining global search heuristics and local parameter space analysis offers a balance of the metaheuristic and multiparametric approach solving challenging classes of bilevel optimization problems, opening new pathways for applications in engineering.

1. Beykal B., et al. Journal of Global Optimization, 2020.

2. Hansen B.P., et al. SIAM Journal on Scientific and Statistical Computing, 1992.

3. Fortuny-Amat J., McCarl B. Journal of the operational Research Society, 1981.

4. Avraamidou S., Pistikopoulos E.N. Computers & Chemical Engineering, 2019.

Alternating Current Unit Commitment (AC-UC) is a challenging large-scale mixed-integer nonlinear problem that has binary variables for generator commitment decisions with full AC network physics. Solving AC-UC directly is computationally expensive due to its non-convex power-flow equations and binary variables. In this work, we considered the Second-Order Cone (SOC) relaxation to Optimal Power Flow (OPF) problem and extended it to UC. Next, we structured methods to use historical data to decide on some set of commitment decisions to accelerate the solution of the full AC-UC problem. The SOC-UC model provides tight convex relaxations, yielding high-quality lower bounds. Using SOC-UC solutions from a collection of perturbed PGLib benchmark instances, we built a library of historical commitments and associated features describing system-level demand and generator costs. Two learning-assisted approaches are developed to predict commitment decisions for new AC-UC instances. The first, based on k-Nearest Neighbor (kNN) method, identifies the most similar SOC-UC training cases and selects their commitment schedules through majority voting. The second uses Support Vector Machines (SVMs) to learn hyperplanes that classify each binary variable into on/off or temporally consistent states based on demand and cost features. For both methods, predicted binary variables are fixed before solving AC-UC, reducing the decision space. Experiments on 50 training and 5 test instances demonstrate that kNN and SVM both recover high-quality commitment patterns while reducing AC-UC solve times. These results illustrate the potential of integrating convex relaxations with learning-based prediction to enable scalable global optimization for power-system scheduling.

Embedding neural networks into optimization models allows us to approximate complex system components while maintaining tractability. However, larger networks with more accurate predictions often make the optimization problem that includes the neural network much harder to solve, resulting in significantly slower optimization solver performance. Although previous studies have discussed individual factors that influence this trade-off, the insights are scattered and lack a unified evaluation.

In this work, we systematically study how different neural network design choices, including the number of neurons, pruning levels, and other architectural parameters, affect both prediction accuracy and the computational performance of the resulting optimization model. Using the same dataset, we train a collection of neural networks with varying configurations and construct an accuracy-size frontier. We then measure how these configurations impact optimization runtime and the quality of solutions relative to ground truth.

Our results highlight cases where smaller or more heavily pruned networks can preserve solution quality while reducing computational cost, as well as cases where additional complexity improves accuracy but slows down the optimization process. This study provides a clearer understanding of the accuracy-complexity trade-offs involved in embedding neural networks in optimization models and offers practical guidance for selecting appropriate neural network architectures under different computational budgets.

Large-scale non-convex mixed integer nonlinear programs (MINLPs) with a block-angular structure arise naturally in multi-scenario optimization contexts, including product and process family design, infrastructure planning, and other stochastic programming applications. While such problems are intractable at scale using off-the-shelf optimization solvers, their structure can be exploited through decomposition strategies to leverage parallel computation. These strategies transform the large-scale MINLPs into smaller subproblems while introducing non-anticipativity (coupling) constraints to ensure consistency of first stage (shared) decisions across all subproblems. To solve the decomposed problem, we employ progressive hedging (PH), implemented through the open-source package mpi-sppy, which enables efficient parallel solution of subproblems across distributed computing architectures. Although progressive hedging provides convergence guarantees for convex problems, its application to non-convex problems remains heuristic. To address this challenge, we introduce a tailored algorithmic design that incorporates dynamic gradient information at each PH iteration to adaptively adjust penalty parameters and accelerate convergence. Furthermore, we integrate domain reduction techniques into the algorithmic framework to tighten bounds on first-stage variables, substantially improving convergence in practice. Through computational experiments on process family design problems, we demonstrate the improvements in using our tailored algorithmic design compared to off-the-shelf solvers.

Multi-agent target tracking in the presence of nonlinear dynamics and agent heterogeneity, where state-space dimensions may differ, is a challenging problem that traditional graph Laplacian methods cannot easily address. This work leverages the framework of cellular sheaves, a mathematical generalization of graph theory, to natively model such heterogeneous systems. While existing coordination sheaf frameworks focus on cooperative problems like consensus, this work extends them to the non-cooperative target-tracking problem. The tracking of multiple, unknown targets is formulated as a harmonic extension problem on a cellular sheaf, accommodating nonlinear dynamics and external disturbances for all agents. A decentralized control law is developed using the sheaf Laplacian, and a corresponding Lyapunov-based stability analysis is provided to guarantee tracking error convergence, with results validated by simulation.

The gradient sliding methods were developed for problems that have separable structures: sum of convex smooth and nonsmooth functions, convex smooth objective and constraint functions, etc. The core idea of gradient sliding-type methods is that the complexity of first-order algorithms should be adaptable to different problem structures. While several gradient sliding methods have been proposed in the literature, few of them are parameter-free; most sliding algorithms require the knowledge of Lipschitz constants. In this talk, we show that gradient sliding methods can be designed without any knowledge of Lipschitz constants. Our results covers optimization problems with structure of the sum of either two convex smooth functions, or a convex smooth function and a convex nonsmooth function. In addition, we also propose novel gradient sliding algorithms for convex functional constrained minimization and monotone variational inequalities with convex functional constraints.

We propose a general optimization-based framework for computing differentially private M-estimators. We first show that robust statistics can be used in conjunction with noisy gradient descent to obtain optimal private estimators with global linear convergence. We establish local and global convergence guarantees, under local strong convexity, showing that our private estimators converge with high probability to a near-optimal neighborhood of the nonprivate M-estimators. We then extend this optimization framework to the setting where each user possesses multiple data points. We develop user-level differentially private algorithms that seek to protect the entire set of data points belonging to a user. Our main algorithm is a noisy gradient descent algorithm, combined with a user-level DP mean estimation procedure to privately compute the average gradient across users at each step. We will highlight the challenges associated with guaranteeing user-level DP and present finite sample global linear convergence guarantees for the iterates of our algorithm.

Constrained nonconvex optimization is central to robotics, optimal control, and machine learning, yet many first-order methods suffer from feasibility violations and brittle penalty or stepsize tuning. In this talk, I present Safe Sequential Quadratically Constrained Quadratic Programming (SS-QCQP), a first-order method for smooth constrained nonconvex optimization that guarantees feasibility at every iteration while ensuring monotonic objective descent. SS-QCQP is derived from a control-theoretic continuous-time dynamical system whose vector field is obtained by solving a convex QCQP that enforces forward invariance of the feasible set and descent of the objective. Under standard constraint qualification conditions, the resulting dynamics achieve an (O(1/t)) convergence rate to first-order stationary points. I then introduce a safeguarded Euler discretization with adaptive step-size selection that preserves these guarantees in discrete time, yielding an (O(1/k)) rate. To improve scalability, I present an active-set variant (SS-QCQP-AS) that enforces only near-active constraints. Experiments on a multi-agent nonlinear optimal control problem demonstrate strong feasibility, predicted convergence behavior, and solution quality comparable to SQP and IPOPT.

We propose a stochastic GDA method with backtracking (SGDA-B) to solve nonconvex-(strongly) concave minimax problems. SGDA-B is agnostic to the Lipschitz constant L, concavity modulus μ≥0, and the variance bound of the unbiased stochastic gradient estimator. Within O(1/ε⁴ log(1/p) and O(1/ε⁷ log(1/p)) stochastic gradient calls, SGDA-B can compute an ε-stationary point in terms of the gradient-map norm and the gradient norm of the Moreau envelope evaluated at the random output point with probability 1-p when μ>0 and μ=0, respectively.

Successful supply chain operations depend on the effective coordination of inventory management strategies and the information required for their implementation. This presentation explores the integration of time-series forecasting with inventory management to assess the value of information sharing in a two-tier supply chain. We analyze how different inventory policies influence supply chain costs under two information-sharing scenarios: (a) full information sharing and (b) no information sharing. We formulate an infinite-dimensional optimization problem where the decision variables are the MA coefficients characterizing a stationary ordering policy. By leveraging the canonical Smirnov-Beurling inner-outer factorization, we derive a fundamental mathematical identity that quantifies the value of information sharing, by viewing an inventory policy as an element of the Hardy space H2. Additionally, we show that an optimal ordering policy inherently incorporates the efficient delay of inventory replenishment in response to market demand. This delay is necessary to gradually release demand information, making the orders more predictable for upstream suppliers. We study how this delay is encoded in an inventory replenishment policy and how to optimally control it to minimize supply chain inventory costs.

We present some hardness results on finding the optimal policy for the static formulation of distributionally robust Markov decision processes. We construct problem instances such that when the considered policy class is Markovian and non-randomized, finding the optimal policy is NP-hard, and when the considered policy class is Markovian and randomized, the robust value function possesses sub-optimal strict local minima. The considered hard instances involve an ambiguity set with only two transition kernels.

In this talk, we study the convexification of sets defined by box constraints and bilinear inequality constraints in which all bilinear terms share a common variable. We present a complete characterization of the extreme points and the facet-defining inequalities of the resulting convex hull. In addition, we develop a separation algorithm for identifying the facet-defining inequalities.

Suppose there is a ground set of n points lying in a 1-dimension space, we want to box as many points as possible while minimizing the cost (size/width of the box). In this case the cost function is close to but not submodular, which makes the function hard to analyze. We show that the convex hull of the epigraph of the function under non-negative constraint can still be obtained by exploiting the combinatorial structure of the function.

The seminal work of Roth, Sonmez, & Unver shows that the Edmonds-Gallai structure theorem for non-bipartite matching can be leveraged to yield a randomized algorithm to match patient-donor pairs in kidney exchange with extraordinarily strong properties. This breakthrough led to randomized polynomial-time algorithms to find maximum cardinality matchings for which individual fairness, measured by the probability that nodes are matched, for example, has maximized Nash social welfare. The societal impact of these results can not be overstated. But the algorithms actually used go beyond cardinality matching, generalizing to weighted variants, as well as also allowing structures such as paths and 3-cycles to cover patient-donor pairs. We first show that strongly polynomial algorithms that guarantee the same fairness properties can be obtained in the weighted setting. The 3-cycle extension, even just for maximum cardinality coverage is NP-hard. We provide a general result that shows that any optimization subroutine for the coverage problem (for whichever structure is allowed) can be bootstrapped using a polynomial number of calls to yield a mechanism that has analogous fairness properties to those obtained for matching. We complement these theoretical results with computational results, both on well-studied synthetic data-sets and on samples drawn from real data, that demonstrate the striking advantages of adding these fairness elements to more general kidney-exchange mechanisms.

In this talk, we study a heterogeneous capacitated vehicle routing problem in which an infinite-capacity main vehicle and a finite-capacity sub vehicle cooperate to serve customers in a planar region. The main vehicle acts as a mobile depot, while the sub vehicle repeatedly departs, collects goods, and returns to rendezvous under a one-way cargo-transfer rule. We first develop a non-asymptotic upper bound on the expected routing cost under a simple traversing-strip policy, avoiding the regression-based tour-length estimators commonly used in the literature. The bound is explicit and captures both vehicle capacity and relative operating costs. Building on this analysis, we design two constructive heuristics: a TSP-based scheme that operates on grid cells and a Hamiltonian-path scheme that reduces redundant depot returns. Computational experiments on large synthetic instances show that both heuristics produce feasible tours whose costs consistently lie below our analytical bound and an adaptation of Daganzo (1984a)’s classical approximation, with the Hamiltonian-based approach achieving the strongest performance and clear economies of scale.

Machine learning is being increasingly used to accelerate branch-and-bound algorithms for mixed integer optimization. Studies in the literature demonstrate success in approximating branching criteria based on strong branching, a local and expensive strategy which empirically leads to small search trees. The state-of-the-art learning methods for this problem are based on neural networks, with many parameters and a need for a large amount of data for training. In this work, we demonstrate the effectiveness of sparse models in approximating strong branching scores, showing that our models remain simple yet powerful even when trained on limited data. Moreover, our approach achieves significant solver speedups without requiring GPU hardware. We demonstrate that our methodology speeds up SCIP when applied to commonly studied mixed-integer linear problems, and a challenging mixed-integer nonlinear problem from the energy sector.

We introduce new integer programming (IP) approaches for the resident-to-rotation assignment problem (RRAP). Given sets of residents, resident classes, and departments, a block structure for each class, staffing needs, rotation requirement, program rules, and resident vacation requests, the RRAP involves finding a feasible year-long rotation schedule that specifies resident assignments to rotations and vacation times. We first present an IP formulation for the RRAP, which can be easily implemented and efficiently solved using off-the-shelf solvers. However, it can lead to disparities in satisfying vacation requests among residents. To mitigate such disparities, we derive a fairness-promoting counterpart that finds an optimal rotation schedule, maximizing the number of satisfied vacation requests while minimizing a measure of disparity in satisfying these requests. Then, we design a computationally efficient Pareto Search Algorithm that finds the complete set of Pareto optimal solutions to the fairness-promoting IP quickly. Additionally, we present a user-friendly tool that implements the proposed models. Finally, we construct divers instances based on data from our collaborator and conduct extensive experiments to illustrate the potential practical benefits of our proposed approaches. Our results demonstrate the computational efficiency and implementability of our approaches and underscore their potential to enhance fairness rotation scheduling.

Multiarmed bandits (MABs) are a widely used framework for sequential decision-making under uncertainty, with significant applications in personalized healthcare. While traditional MAB methods typically assume stationary reward distributions, real-world settings are often nonstationary. For instance, in personalized interventions, treatment effectiveness may decline due to habituation and later recover, creating complex nonstationary behaviors that standard methods fail to capture.

In addition to adapting to these evolving reward structures, a critical challenge—especially in clinical trial contexts—is preserving statistical power for reliable inference. Existing bandit algorithms generally optimize cumulative rewards but often neglect adequate exploration which is necessary for detecting meaningful treatment effects.

We propose a novel Posterior Sampling algorithm explicitly designed to tackle this dual challenge. Our method adapts to nonstationary environments characterized by habituation and recovery dynamics, while simultaneously satisfying specified power constraints through controlled exploration.

We provide theoretical guarantees via sublinear regret bounds and validate the approach on two real-world datasets by simulation: a bipolar disorder clinical trial and a microrandomized trial promoting workplace standing. Results demonstrate the method’s effectiveness in personalizing interventions while maintaining accurate statistical inference.

Effective invasive species management requires strategies that balance ecological impact with limited resources. We present a discrete-time mathematical framework to model the spatiotemporal spread of Chinese privet (Ligustrum spp.), incorporating seed- and root-mediated dispersal across heterogeneous landscapes. The model captures nonlinear spatial interactions, dynamic invasion fronts, and varying habitat suitability. To inform practical management, we integrate the mathematical model with combinatorial optimization techniques, identifying adaptive, cost-effective removal strategies that minimize reinfestation risk while adhering to operational and budget constraints. Our approach accounts for treatment thresholds, spatial prioritization, and multi-period planning, enabling managers to assess trade-offs between short-term actions and long-term ecological outcomes. By coupling invasion dynamics with optimization-based decision-making, this framework provides a rigorous, data-driven tool for prioritizing interventions and designing scalable, ecologically effective invasive-species management plans.

In many real-world decision-making problems, data arrives from sources that differ widely in fidelity and acquisition cost. Effectively leveraging such heterogeneous observations is crucial for both rapid learning and optimization. On one hand, presence of low fidelity data lead to suboptimal decisions, and on the other hand, the cost of high-fidelity data may not justify its impact on convergence. In this talk, I will present two results illustrating how incorporating data quality yields provable and practical gains.

First, in a stochastic optimization setting, we consider two gradient oracles of different costs and noise levels. For e.g., the cost of querying foundational models depends on the number of tokens specified in their input. We develop an adaptive algorithm that decides which oracle to query based on the state of the optimization process. Under a fixed budget, we show how to obtain strictly faster convergence rates compared to strategies that rely on a single oracle throughout.

Next, in a multi-armed bandit setting, we consider the feedback from N arms, which is dependent on K mechanisms used to query the arms. For example, in a call center, sales depend not only on customer features (arms), but also on the quality of the agents (mechanisms) who make the calls. We first show that ignoring the mechanism quality can lead to linear regret. We next show how to optimally use low-fidelity data samples in this framework, and discuss the impact on a real-world case study at Blackthread Analytics.

This work investigates the resilient consensus problem over undirected graphs, incorporating stepsizes to mitigate the influence of malicious agents. Using trust observations, legitimate agents identify malicious agents and exclude them from their communication. For this purpose, agents form their trusted neighbors by comparing their cumulative trust values at each time step. Compared to prior works, we introduce stepsizes into the consensus process that progressively slow the consensus updates. This is achieved by using stepsizes to combine the standard consensus update with the agent's current value, thus attenuating the malicious effect. In the absence of finite-time detection guarantees, we show that legitimate agents reach agreement in expectation. This result holds provided the stepsizes are non-summable and the limit probability of legitimate agents correctly identifying their neighbors' types converges to one. We further derive bounds on the probability of deviation from the nominal consensus process, which is the ideal case featuring no malicious agents and no stepsizes. This bound is a function of the sum of products of stepsizes and misclassification probabilities. Numerical experiments are presented to verify the theoretical results.

We consider the problem of solving high-dimensional statistical inference problems over a network of agents (with no coordinating agent) who have exclusive access to a fraction of the total available samples. In the high-dimensional setting, the problem dimension is much larger than the total number of available samples, making the problem ill-conditioned. Despite this, we empirically observe that obtaining a statistically meaning solution is possible with many existing decentralized schemes, given that underlying parameter to estimate lies in a low dimensional subspace. Our observations challenge the existing theories in two key ways: (i) most decentralized schemes do not break down as the problem dimensionality increases, and (ii) decentralized schemes that are expected to behave like one another behave very differently in high dimensions. To understand the behavior of decentralized optimization methods in high-dimensional inference we introduce a unified framework and analysis, allowing us to develop an understanding of the mechanisms enabling dimension independent performance of decentralized schemes.

Interpretable machine learning algorithms are prized for their ability to extract useful information from data and to audit predictions for trust and safety. Formulating interpretable models, such as sparse additive models or compact rule sets, using discrete optimization approaches leads to impressive modeling flexibility, offering exact control over sparsity and interaction structure. In this talk, we overview how recent advancements in discrete optimization frameworks for ensemble pruning and additive models make it possible to construct inherently transparent models that perform competitively with black-box ensembles. We discuss various formulations that encourage different behaviors in the constructed models, as well as specialized algorithms to solve the resulting problems efficiently.

Despite the empirical success of Low-Rank Adaptation (LoRA) in fine-tuning pretrained models, there is little theoretical understanding of how first-order methods with carefully crafted initialization adapt models to new tasks. In this work, we take the first step towards bridging this gap by theoretically analyzing the learning dynamics of LoRA for matrix factorization (MF) under gradient flow (GF), emphasizing the crucial role of initialization. For small initialization, we theoretically show that GF converges to a neighborhood of the optimal solution, with smaller initialization leading to lower final error. Our analysis shows that the final error is affected by the misalignment between the singular spaces of the pre-trained model and the target matrix, and reducing the initialization scale improves alignment. To address this misalignment, we propose a spectral initialization for LoRA in MF and theoretically prove that GF with small spectral initialization converges to the fine-tuning task with arbitrary precision. Numerical experiments from MF and image classification validate our findings.

In this study, we propose building an optimization agent based on mix of LLM and Neural Network architecture which first transforms NLP problems written in text format to the KKT conditions. In second step, the solution to these nonlinear KKT equations is posed as a residual minimization problem. This minimization problem is then encoded as a Neural Network (NN) and solved efficiently using the gradient descent based backpropagation algorithm. This step doesn’t require training of neural network using artificial dataset and also mitigates the need for custom initialization associated with the classical Newton-Raphson method. We test our approach on benchmark NLP problems to showcase and compare it’s performance with conventional NLP solvers.

We study a distributionally robust generalization of the classic discrete-time, finite-horizon Linear-Quadratic-Gaussian (LQG) control problem. The state and observation noise distributions are chosen adversarially from divergence-based ambiguity sets centered at a nominal distribution that is zero-mean and Gaussian. The distribution of noise terms may have nonzero means and must have finite second moments satisfying an orthogonality condition (which reduces to uncorrelatedness when means are zero). Under a mild structural assumption on the divergence – satisfied by several relevant examples – we prove that a policy depending affinely on the history of observations is optimal, and the adversary’s worst-case distribution is Gaussian, and with an additional weak assumption, we prove that the adversary sets the mean to zero and the optimal policy becomes linear. We develop a Frank-Wolfe algorithm based on solving standard LQG subproblems via Kalman filtering and dynamic programming, and we show empirically that it consistently outperforms semidefinite-programming reformulations. We show that the results hold when the nominal distribution is elliptical – rather than Gaussian – if the divergence is the 2-Wasserstein metric.

This talk presents an infinite horizon, average cost formulation of distributionally robust linear quadratic Gaussian control with time invariant system matrices and a time invariant nominal disturbance model. We study a zero sum game between a controller and an adversary who selects disturbance distributions from divergence based ambiguity sets centered at a Gaussian nominal model. Importantly, we do not impose stationarity a priori on the control policies or on all distributions in the ambiguity set, which introduces additional technical challenges and requires slightly stronger assumptions than in the finite horizon case. Under mild regularity conditions on the divergence, stabilizability and detectability that mirror classical LQG assumptions for infinite horizon models, we establish the existence of a Nash equilibrium with a tractable and interpretable structure: the controller employs an optimal stationary linear output feedback policy, while the adversary selects a worst case time invariant Gaussian distribution. The analysis relies on Fan type minimax duality and shows that classical LQG separation principles extend to this distributionally robust setting, thereby generalizing our finite horizon results and aligning the robust model with the structural insights of traditional infinite horizon LQG with known noise distributions.

Joint work with Huiyi Cao.

In chemical engineering applications such as safety verification, nonconvex optimization problems must be solved to guaranteed global optimality, to within a specified tolerance. Typical methods for deterministic global optimization compute crucial bounding information by minimizing convex relaxations of objective functions. However, generating useful convex relaxations can be difficult if process models include embedded systems of ordinary differential equations (ODEs), since established reachable-set methods for ODEs are impressive yet still somewhat limited in scope and fidelity. These dynamic relaxation methods compute convex relaxations for ODE solutions together with interval bounds, by solving an auxiliary relaxation ODE system. By analogy, in the well-known McCormick relaxation approach for constructing convex relaxations of composite functions, these relaxations are also computed in tandem with constant interval bounds, which provide useful range estimates for intermediate variables.

We present a framework that extends state-of-the-art relaxation approaches for ODE solutions, with the ultimate goal of aiding methods for dynamic global optimization. This framework makes use of the relationship between ODE solution relaxations and their related interval bounds. We take mathematical techniques that were used to establish correctness of ODE solution relaxations, and we adapt these to develop several new types of interval bounds on ODE solutions, with guaranteed inclusion monotonicity properties. (We also derive several established methods in a new way.) For example, we obtain new effective interval bounds based on the edge-concave relaxations of Hasan (2018), which significantly outperform convex relaxations in some cases. Numerical examples are presented, and implications are discussed.

We present a new relaxation based on separable edge-concave underestimator (SEC) for general n-dimensional signomial problems. While the vertex polyhedral convex envelope of an edge-concave underestimator is constructed with convex combination of function evaluations at all 2ⁿ domain vertices, the SEC underestimator uses a single linear hyperplane constructed at only (n+1) linearly independent vertices. We show that there exists a single n-dimensional hyperplane passing through all 2ⁿ vertices of any separable function defined over a box. Moreover, we provide an explicit form of the linear hyperplane representing the convex envelope of the SEC underestimator. Having an explicit form of the hyperplane leads to an efficient linear programming relaxation of the original problem with high-dimensional signomials. We test our proposed relaxation technique to obtain root node lower bounds for several test functions with dimensions up to 100. We find that our method of relaxation, on average, provides tighter and faster lower bounds for higher dimensional signomials when compared with the performance of state-of-the-art commercial solvers. Future work entails to extending the separable edge-concave underestimators for general C² nonconvex problems.

Specialized computing architectures offer opportunities to dramatically boost the computational performance of algorithms for scientific computing applications. One dominant strategy is the use of graphics processing units (GPUs) for algorithms with highly repetitive floating-point operations and high data throughput, such as machine learning model training and data analysis. Such computing architectures are classified as “single-instruction, multiple-thread” or SIMT, and effectively execute a single instruction across many computing elements simultaneously. In contrast to conventional computing architectures (e.g., x86 and x86-64), which are classified as “multiple-instruction, multiple-data” or MIMD, that allow multiple independent and different instructions to be executed simultaneously on different data. Until recently, the acceleration of mathematical programming algorithms was not possible with SIMT architectures. Recently, MadNLP.jl was developed to dramatically accelerate the solution of large-scale nonlinear programs to local optimality. However, for applications where guaranteed optimality is paramount, such as in robust design, deterministic global optimization is required. In this talk, we discuss our recent advances in the development of a novel deterministic global optimization algorithm and software implementation for SIMT architectures. Our approach utilizes McCormick-based relaxations and their subgradients to construct linear programming lower-bounding problems. A novel implementation of PDLP was developed to solve LP lower-bounding problems entirely on the GPU with minimal data transfer between the host CPU and the GPU. Implemented within a spatial branch-and-bound framework, this approach achieves a significant speedup over an equivalent CPU implementation.

The Panama Canal is a vital hub for international shipping, and the scheduling of canal operations is critical for the efficiency and sustainability. As with other maritime bottlenecks, the Panama Canal currently requires most scheduled ships to arrive hours ahead of their transit, creating queues of ships to prevent canal idle time. Conversely, ship operators can reduce emissions and costs by adjusting sailing speed and arriving Just In Time (JIT). Aligning with the Canal’s goal of encouraging sustainable supply chains, JIT arrivals present an opportunity to reduce emissions. However, it remains unclear if emission-reducing policies will impact the canal’s revenue, which is primarily driven by canal throughput. In this study, we propose a combined ship scheduling and speed reduction model as a mixed integer linear program (MILP) for evaluating how changes to the canal’s scheduling policy can impact the sailing speed and thereby emission from ships transiting the canal. Our heuristic decomposition method ensures an acceptable computation time. We validate the model’s output applying historical data from the Panama Canal Authority (ACP) and the Automatic Identification System (AIS). The results demonstrate significant emission reduction potential, consistent with prior research. Furthermore, by comparing profit and emission reducing objectives, the results suggest that a change in policy is practically feasible without a significant difference in canal throughput, which is the primary driver of revenue. With acknowledged limitations, the results indicate that the proposed model is suitable for evaluating the impact of policy changes to the ship scheduling procedures at the Panama Canal.

Crowd-based bike delivery systems offer a scalable, cost-effective alternative to traditional restaurant delivery models, which rely on costly owned-fleet vehicles. This research evaluates the viability of such systems in Dublin through an agent-based simulation model integrating population characteristics, biker behavior, and restaurant operations.The framework combines a computational agent-based simulation with multiple optimization modules to guide agent interactions. A biased-randomization savings heuristic for Vehicle Routing Problems, adapted from traditional fleet-based systems, optimizes delivery routes within the simulation. The model explores diverse scenarios varying demand rates, biker adoption, and operational parameters. Results demonstrate promising performance across key performance indicators including economic efficiency, delivery reliability, and customer satisfaction. This work contributes to understanding the viability of sharing-economy delivery models and provides practical optimization insights for urban logistics networks operating at the intersection of data-driven decision-making and sustainable mobility solutions.

To enable interoceanic transit, the Panama Canal uses freshwater. Because of the mountainous geography, water from a lake at a higher altitude than each entrance to the canal facilitates vessel passages. With each vessel crossing, millions of gallons of freshwater are discharged to the oceans. Thus, increasing transits means greater water loss. This is the fundamental trade-off of the canal’s design.

In the past, water availability was not a concern, due to Panama’s rainforest climate. But severe droughts are risking its service reliability, creating rippling economic effects on global supply chains. The canal requires precise water conservation strategies for every scheduled vessel, to maintain its watershed and attain its business goals.

Scheduling vessels through the canal is highly constrained; not only by its operational capacity but also by opposing transit directions, pre-emptive processing, and cargo-related navigational restrictions.

To find operational efficiencies while considering conflicting goals, we propose a multi-objective bidirectional hybrid flow shop model that considers two water conservation modes. We define a Pareto frontier of non-dominated solutions that balances makespan and water utilization. As an alternative to drastically restricting transits during droughts, the Panama Canal Authority’s watershed-level policy can be represented in this model, and water conservation proactively implemented.

Child trafficking remains a critical social challenge marked by complex spatial and contextual heterogeneity. This study employs a Bayesian nonparametric modeling framework to investigate trafficking vulnerability patterns across Sierra Leonean chiefdoms, administrative units analogous to cities in the United States, within three high-burden districts. Incorporating areal information from the chiefdoms as random effects, a generalized linear mixed model (GLMM) is developed to assess the relationship between trafficking prevalence and potential risk factors. By applying a nonparametric Bayesian clustering approach to the random effects, the study innovatively identifies latent risk structures and spatial dependencies. Specifically, covariate-dependent models such as the Weighted Dirichlet Process Mixture (WDPM) and Dependent Dirichlet Process (DDP) models are leveraged to integrate socioeconomic and cultural covariates, thereby enhancing clustering precision and estimation accuracy. Unlike traditional clustering techniques, the Bayesian nonparametric framework adaptively determines cluster configurations without requiring a predefined number of groups. Clustering outcomes are rigorously validated through community-level hypothesis testing to ensure robustness and reliability. This novel modeling framework provides a comprehensive and data-driven approach to understanding trafficking risks in resource-limited settings. By uncovering hidden vulnerability patterns, the study delivers actionable insights for designing targeted, context-specific interventions. The findings underscore the potential of advanced Bayesian modeling to inform evidence-based policy-making and optimize resource allocation in efforts to combat child trafficking effectively.

Illicit Massage Businesses (IMBs) are a covert and persistent form of organized exploitation, operating under the guise of wellness services while facilitating human trafficking, sexual exploitation, and coerced labor. Detecting IMBs is challenging because traffickers rely on coded digital advertisements, rapidly rotating staff, shifting business locations, and shared infrastructure such as phone numbers, addresses, and aliases. Traditional detection strategies—community tips, regulatory inspections, and manual reviews—are reactive and struggle to uncover the broader operational networks that enable IMBs to persist.

To address these limitations, we introduce IMBWatch, a spatio-temporal graph neural network (ST-GNN) framework designed to identify IMB activity at scale. IMBWatch builds dynamic graphs from open-source intelligence, where nodes represent entities such as businesses, individuals, phone numbers, and aliases, and edges encode relational, spatial, and temporal dependencies. The model integrates graph convolutional operations with temporal attention mechanisms to capture how IMB networks evolve across time and geography, enabling the discovery of subtle but consistent behavioral patterns.

Empirical evaluations on real-world data from multiple U.S. cities show that IMBWatch outperforms baseline approaches, achieving higher accuracy. The system produces actionable insights that support more targeted and proactive investigative strategies. IMBWatch is scalable, adaptable to other trafficking-related domains, and will be made publicly available with anonymized datasets and open-source code. This makes it a practical, data-driven tool for strengthening anti-trafficking efforts and advancing research on illicit network behavior in complex, evolving environments.

First-price auctions are one of the most popular mechanisms for selling goods and services, with applications ranging from display advertising to timber sales. Unlike their close cousin, the second-price auction, first-price auctions do not admit a dominant strategy. Instead, each buyer must design a bidding strategy that maps values to bids—a task that is often challenging due to the lack of prior knowledge about competing bids. To address this challenge, we conduct a principled analysis of prior-independent bidding strategies for first-price auctions using worst-case regret as the performance measure. First, we develop a technique to evaluate the worst-case regret for (almost) any given value distribution and bidding strategy, reducing the complex task of ascertaining the worst-case competing-bid distribution to a simple line search. Next, building on our evaluation technique, we minimize worst-case regret and characterize a minimax-optimal bidding strategy for every value distribution. We achieve it by explicitly constructing a bidding strategy as a solution to an ordinary differential equation, and by proving its optimality for the intricate infinite-dimensional minimax problem underlying worst-case regret minimization. Our construction provides a computationally-tractable procedure for deriving minimax-optimal bidding strategies. When the value distribution is continuous, it yields a deterministic strategy that maps each value to a single bid. Importantly, our result allows us to quantify, through minimax regret, the performance loss due to a lack of knowledge about competing bids. We leverage this to analyze the impact of the value distribution on the performance loss.

We study a basic but unresolved question in stochastic optimization: when should stochastic gradient descent (SGD) be stopped based only on its observed trajectory? We develop anytime-valid confidence sequences for stochastic gradient methods that remain valid under continuous monitoring and directly yield statistically valid stopping rules. In convex optimization, they certify weighted suboptimality under general stepsize schedules; in nonconvex optimization, they certify weighted first-order stationarity. The result is a unified framework for online stopping of SGD with provable complexity guarantees

While recent work in OR/OM increasingly explores /how/ to better train AI algorithms for decision-making, less attention has been given to a complementary question: once AI predictors are successfully trained, what rules has the AI actually learned? This talk demonstrates how the geometry of a network’s learned embeddings can shed light on this issue. Using representative OR decision-making problems, such as multi-class scheduling, we show that after an AI learns to solve these problems, its internal vector representations can be mapped to classically understood models. This mapping allows for interpretability in the AI's decision-making procedure.

First-order methods are the workhorse of modern machine learning, but they are notoriously inefficient in non-convex landscapes, frequently becoming trapped by local minima and saddle points. These classical limitations make them a prime target for quantum acceleration. In this talk, I will first establish a fundamental connection between gradient flows (i.e., the continuous-time limit of gradient descent) and the evolution of a continuous-variable quantum system under a corresponding Hamiltonian. This mapping yields a new class of quantum algorithms that not only inherit the efficiency of gradient descent but also leverage quantum tunneling to escape saddle points and local minima. I will present results demonstrating this quantum advantage, including provable speedups and superior empirical performance in finding global solutions. Finally, I will briefly survey some recent advances in generalizing this quantum dynamical approach to higher-order optimization methods.

Federated learning enables distributed model training without sharing raw data. In this framework, keeping the clients engaged remains challenging since participants may receive uneven benefits from their contributions. We consider nonconvex and L-smooth local objective functions for the clients. We present an Incentive-enabled Federated Averaging method that rewards participants based on their level of participation in each round, through allowing clients to control the size of their training data dynamically in each round of communication. In each communication round, clients exchange both local parameters with the server and their updated training dataset sizes, which are dynamically adjusted using a Nash equilibrium seeking update rule. Subsequently, the server transmits updated global parameters back to clients to initiate the following round. We theoretically analyze the proposed method, obtaining explicit iteration and communication complexity guarantees. We plan to conduct numerical experiments on MNIST and CIFAR-10 datasets to evaluate the effectiveness of our approach, with a focus on examining the convergence behavior of the loss function.

The integration of neural networks into mathematical optimization is well established for simple architectures but remains challenging for complex models such as the Transformer. The main difficulty stems from the self-attention mechanism, whose quadratic interactions and non-convex softmax function complicate direct optimization. This paper presents the first rigorous mathematical programming formulation that captures the computational structure of a single self-attention head.

We develop a reduced, abstract representation that explicitly models the key nonlinearities of self-attention. To ensure tractability, we introduce systematic piecewise-linear approximations for the softmax and interaction terms, enabling the formulation to balance fidelity and computational efficiency. The resulting model preserves the expressive behavior of attention while keeping the problem size manageable for modern solvers.

This formulation provides a critical building block for embedding high-fidelity, pre-trained Transformer components within Mixed-Integer Nonlinear Programming (MINLP) frameworks. By enabling optimization models to directly incorporate attention-based surrogate functions, this work advances the integration of deep learning and mathematical programming and supports the development of new hybrid AI–optimization methodologies.

We consider the problem of minimizing a difference-of-convex objective over a nonlinear conic constraint. We assume that the cone is closed, convex and pointed with a nonempty interior, and the support function of a compact base of its dual cone admits a so-called majorizing smoothing approximation. These conditions are satisfied by widely studied cones such as the nonnegative orthants and the cone of positive semidefinite matrices. We reformulate the conic constraint equivalently as a single inequality constraint involving the aforementioned support function, and adapt the moving balls approximation (MBA) method for its solution. In essence, in each iteration of our algorithm, we approximate the support function by a smooth approximation function and apply one MBA step. The subproblems that arise in our algorithm always involve only one single inequality constraint. We design explicit rules to evolve the smooth approximation functions from iteration to iteration and establish the global iteration complexity for obtaining an approximate Karush-Kuhn-Tucker point. In addition, in the convex setting, we establish convergence of the sequence generated, and study its local convergence rate under a standard Hölderian growth condition.

This is a joint work with Nung-sing Sze and Jiefeng Xu.

Knitro was originally developed in the 1990s as an interior-point algorithm for general nonlinear, non-convex optimization. Over the years, Knitro evolved into a more general optimization toolbox that includes an Active-Set LP-based solver, a Sequential Quadratic Programming (SQP) solver, specialized LP, QP, and SOCP solvers, a branch-and-bound solver for mixed integer programming, and multi-start heuristics for global optimization. To add to this toolbox of algorithms, we have recently started developing first-order methods that do not require any matrix factorizations. The hope is that these might be able to provide useful solutions to extremely large-scale models where the factorizations in interior-point methods become too expensive. In this talk we will present some of this work, primarily based on Augmented Lagrangian (AL) type methods.

We introduce the concept of a Lagrangian game as a general framework for solving constrained Markov games (CMGs). CMGs model multiagent systems operating in dynamic environments under global or local feasibility constraints. These constraints encode safety, resource, and regulatory requirements and serve as the foundation for safe multiagent reinforcement learning (safe MARL). In our approach, agents repeatedly construct and solve an unconstrained Lagrangian game based on their current multipliers, simulate finite state-action trajectories to estimate cumulative rewards and constraint violations, and update the multipliers using a gradient step. This iterative process generates a sequence of Lagrangian games and corresponding policy profiles. We show that the resulting sequence of policies forms a nonstationary, approximately Nash solution for the original CMG and achieves constraint feasibility with controlled accuracy. Our analysis covers both general-sum and cooperative CMGs, including constrained multiagent MDPs, where the method yields nonstationary solutions that outperform any stationary policy up to a quantifiable error. We further extend the framework to locally constrained Markov games, where agents have individual private constraints and update multipliers using local information, while still converging to feasible, nonstationary approximate Nash equilibria.

The exploration-exploitation dilemma in reinforcement learning is a fundamental challenge to the design of efficient algorithms. Existing policy gradient methods address this by either positing or explicitly injecting exploration, which can yield sub-optimal performance and result in non-implementable methods. This work introduces a new class of auto-explorative methods. These are methods that automatically explore both state and action spaces without any exploration strategies while adapting to the true, unknown problem parameters like the mixing time of the induced Markov chain. The resulting parameter-free methods do not estimate these unknown parameters and instead utilize a novel and easy-to-compute certificate. For finite state and action Markov decision processes, we show these methods achieve the same worst-case complexity as when the parameters are known. Moreover, we extend these auto-explorative methods to the case where linear function approximators are used and demonstrate a similar complexity. Such results can potentially be applied to general state and action spaces as well.

Physics-embedded chemical process flowsheets exhibit inherent nonconvexities, trapping local solvers and rendering global optimization computationally intractable. Traditional surrogates built using simulation data are often difficult to optimize and fail to ensure material and energy balance closure.

We propose an alternative framework, first solving the cost optimization problem for individual unit operations, generating optimal designs across diverse feed conditions and product targets. A surrogate model is then constructed to map these inputs (feed conditions, product targets) directly to the optimal cost. We compare surrogate types (e.g., polynomial regression, linear decision trees, neural networks), analyzing their accuracy versus the dimensionality of the resulting surrogate-embedded optimization. Furthermore, multi-unit flowsheets contain various connectivity constraints to direct material and energy flows. We exploit dependent connectivity constraints to reduce the degrees of freedom (DOF). A new surrogate using DOF as inputs predicts interconnecting stream properties. These predictions feed the unit-level cost surrogates, directly linking flowsheet connectivity DOF to total cost. This approach significantly reduces the optimization problem’s dimensionality, enabling the solution of larger-scale flowsheets.

This problem-agnostic hierarchical decomposition methodology provides a powerful framework for a broad class of nonconvex network optimization problems, such as water distribution networks, gas networks, and supply chain problems.

Over the last decade, there has been growing interest among researchers in studying sparse methods for high-dimensional data. At the heart of this task lies the Best Subset Selection (BSS) in linear regression, a classical problem in statistics that exemplifies the computational difficulty of solving sparsity-constrained optimization problems. The BSS problem is known to be NP-hard, and finding the global optimum becomes increasingly difficult as the problem dimension grows. Branch-and-bound (BB) methods provide a general framework for solving BSS. However, in the underdetermined case the bound-based pruning is less likely to occur early in the search tree of any BB method, making it more computationally extensive. We propose a novel BB method that employs a specialized enumeration to generate candidate solutions and guarantee convergence to the global optimum. As with any BB method, its efficiency critically depends on how memory usage is controlled and how strong candidate solutions are identified during the search. Our method leverages the echelon form of each parent node to efficiently compute lower bounds for a sequence of child nodes in lower dimensions, making it particularly effective in the underdetermined case. We compare our method with an existing BB method and with GUROBI’s quadratic mixed-integer solver.

Gaussian process (GP) models are widely used nonparametric surrogates across science and engineering due to their flexibility, ability to represent uncertainty in a calibrated way, and rigorous treatment of noisy observations. Many surrogate-based optimization methods embed GP posterior mean and variance terms, leading to challenging global problems, especially when minimizing the popular lower confidence bound (LCB) acquisition function (which systematically balances exploration and exploration). Spatial branch-and-bound (BnB) methods offer a principled framework for solving such problems deterministically, but their performance hinges on constructing relaxations that are both tight and computationally tractable. Existing approaches, such as those based on McCormick relaxations, often yield loose bounds for GP posteriors and scale poorly with the number of training points. We introduce a specialized BnB algorithm that combines closed-form analytic bounds with piecewise-linear relaxations for selected kernel contributions, resulting in mixed-integer (quadratic) relaxations that can be solved efficiently using modern MIP solvers. Our method applies to a broad class of stationary kernels and accommodates both the O(N) posterior mean and O(N²) variance terms. Benchmark results in 4-10 dimensions demonstrate significant reductions in CPU time and BnB tree size relative to BARON, SCIP, and MAiNGO.

Designing experiments is essential for drawing reliable conclusions in research. A key challenge lies in dividing participants into groups that are balanced and representative of the population. Traditional methods, such as randomization and re-randomization, are widely used to reduce bias, but they do not guarantee optimal group balance, particularly in studies with small samples or heterogeneous, multi-attribute populations. Optimization-based methods have emerged as a promising alternative, enabling researchers to directly target balance criteria across several covariates simultaneously. In this talk, we present new multi-objective optimization models that integrate diverse discrepancy metrics, including moment differences, Wasserstein distance, variation distance, and a new bin-based metric. Through computational experiments, we evaluate the strengths and limitations of each method in achieving group similarity and empirically show that our multi-objective approaches demonstrate promising performance.

The Bachillerato Inicia UC program offers a pathway for students from Chilean technical high schools to articulate into undergraduate programs at Pontificia Universidad Católica de Chile. Upon applying, candidates rank up to three preferred programs. However, articulation is determined only after one year, based on their academic ranking within the cohort - a value unknown at the time of admission. The challenge is to select a feasible subset of students to admit, respecting the order of applicant scores, while guaranteeing that, under any realization of final rankings, each admitted student can be matched to one of their declared preferences, respecting program capacities. We formalize this as a two-stage stable matching problem under ranking uncertainty and design an algorithm that characterizes the set of robustly admissible candidates. The model ensures articulation guarantees that are consistent with declared preferences and capacity constraints. We present implementation results from the 2024 and 2025 admission cycles.

We address a class of integer programs on graphs with total variation regularization. Total variation regularizers are common in problems such as image denoising and topology optimization for their ability to smooth solutions while still preserving sharp edges. Our methods are centered around augmentation algorithms which maintain a set of moves which repeatedly "augment" a given feasible point until it cannot be improved any further using such moves. In the context of total variation-regularized problems, we show the moves which update the solution on a connected subgraph are a Graver basis, meaning they act as an optimality certificate, and the augmentation algorithm produces a globally optimal solution in certain cases. Unfortunately, this set of moves is so large as to be computationally intractable to search completely, so we substitute an exhaustive search with a randomized approach which samples connected subgraphs. We demonstrate how to randomize the simplex algorithm applied to an associated linear program to produce random connected subgraphs and how to bias the sampler to produce moves which are more likely to improve the current iterate's objective. We then apply our method to some problems from topology optimization and show it is competitive with the standard approaches.

There exist unique challenges to reliably managing supply chains electric vehicle (EV) batteries. They require assembly from multiple raw materials sourced from across the world. Furthermore, raw materials are often processed at locations that are different than the final production stage, which adds more complexity. There is a high degree of uncertainty in the price for these technologies because of the ever-changing landscape of the energy sector. Moreover, the uncertainty from market and geopolitical factors impacts the ability to ship materials from one location to another and is exacerbated by the highly interconnected nature of battery supply chains. To address these challenges, we develop a risk-averse stochastic optimization methodology to provide secure operation of the energy storage technology supply chain. Our model considers not only the raw material flow, but also the intermediate processing and the final production stage. Our model also allows for each stage to take place in different locations. We additionally incorporate stochastic components such as uncertainties caused by pricing. Finally, our model optimizes the Conditional Value at Risk metric, which allows us to produce results that are more reliable compared to traditional stochastic optimization methods. We present experimental results that demonstrate our methodology’s efficiency and its ability to reliably provide cost-effective solutions. A major consequence of our methodology is that it can allow EV producers to securely obtain necessary technologies.

We introduce an adaptive two-stage stochastic optimization model for energy infrastructure expansion planning under uncertain demand. Unlike traditional two-stage models that are too rigid or multi-stage models that are overly flexible, this new model balances commitment and flexibility by allowing investment decisions to be revised at specific adaptation times. The model’s performance is evaluated using two metrics and compared to conventional approaches. To improve its practicality, this paper also presents five strategies for sharing adaptation times among related investments, which further reduces the number of policy revisions. A case study on Rwanda’s electrification plan demonstrates that this approach can lead to cost savings of up to 3.44% compared to the two-stage model and significantly fewer adaptations compared to the multi-stage model. The model results also provide valuable policy implications, such as the optimal adaptation frequency and timing, as well as the ideal planning horizon length.

Certifying the quality of solutions to risk-averse stochastic optimization problems remains a central challenge in data-driven decision making. This difficulty is especially pronounced in Conditional Value-at-Risk (CVaR) models, where the tail-focused nature of the objective makes standard asymptotic SAA (Sample Average Approximation) guidance unreliable for determining when a computed solution is sufficiently accurate. To address this limitation, we introduce a comprehensive finite-sample certification framework that provides explicit, high-confidence optimality guarantees for CVaR-based SAA solutions. The framework generates rigorous upper bounds on the optimality gap that can be incorporated into existing SAA workflows without altering model structure. To corroborate the effectiveness of this approach, we apply the method to risk-aware epidemic control problems, where mitigation policies must be selected under significant structural uncertainty. Our finite-sample certificates enable decision makers to reduce sampling effort while maintaining conservative performance guarantees, which is critical in time-sensitive public health contexts where rapid, reliable assessment of risk-averse policies is needed. Numerical studies demonstrate that the approach reduces sampling effort while maintaining conservative guarantees, thereby enabling reliable risk-averse decision support in time-sensitive epidemic settings.

Urban air mobility (UAM) provides a promising solution for alleviating ground traffic congestion and achieving zero direct aviation emissions. The successful and reliable implementation of UAM requires careful coordination between regulatory authorities, such as the Federal Aviation Administration (FAA), and airline operators, particularly in uncertain conditions of passenger demand and airspace capacity. In this paper, we introduce a bilevel risk-averse optimization framework that integrates airline routing decisions with FAA capacity management and delay assignment at critical airspace bottlenecks. To address the computational complexity of this integrated routing and resource allocation problem, we propose and evaluate both heuristic and exact optimization methods. Numerical experiments demonstrate the effectiveness of the proposed approaches and provide valuable insights for stakeholders seeking efficient, equitable, and robust UAM operations in uncertain environments.

Machine learning tools like classifiers are now ubiquitous in decision-making for many high-stakes applications like loan approvals, granting bail, admissions and hiring, etc. When agents receive unfavorable outcomes from the classifier, it naturally creates incentives for them to modify their features strategically to improve future outcomes. Such modifications may either be genuine improvement or gaming or a combination of both. However, this strategic behavior is problematic: the classifier is now deployed on data that does not match the training-time distribution, and this can significantly hurt its accuracy. The emerging area of ``strategic classification" studies this problem, with the goal of designing classifiers that are robust to distribution shifts induced by strategic behavior. However, an equally important yet under-explored goal is to ensure that the design of the classifier also incentivizes agents to perform desirable modifications in the first place. Optimization theory provides natural tools to model, formalize, and tackle such design problems. In this work, we borrow heavily from these tools to explore the complexity of designing "good" classifiers---accurate classifiers that also incentivize desirable agent behavior---, under varying levels of modification power and classifier information available to the agents. Our work demonstrates the enormous power of optimization in addressing emerging challenges in other disciplines.

For long-horizon sequential decision-making problems, standard policy-gradient methods typically operate at a single temporal scale, which leads to high variance and weak credit assignment, especially when returns are driven by delayed consequences of earlier choices. We study policy optimization for such problems through the lens of the options framework, a hierarchical architecture that introduces temporally extended decisions on top of primitive actions. Building on this structure, we derive coupled gradient estimators that decompose returns into high- and low-level advantages, where the high-level advantage captures coarse progress over extended time intervals and the low-level advantage refines primitive decisions toward the goal. We analyze the resulting gradient estimators from an optimization perspective, highlighting how temporal abstraction improves the conditioning of the effective objective and reduces estimator variance. Across environments, we observe more stable training dynamics and improved sample efficiency compared to flat actor–critic baselines with standard generalized advantage estimation. Overall, the results indicate that our proposed gradient estimators, instantiated under the hierarchical options framework, provide a principled and effective way to exploit temporal abstraction in long-horizon settings and suggest a path toward more scalable training schemes for complex decision-making systems.

This paper proposes a Benders decomposition framework for the Optimal Transmission Switching with De-energization (OTSD) problem, which ensures secure power grid operations by preventively disconnecting power lines. A key obstacle for solving OTSD at scale is the representation of potential loss of network connectivity ---and the resulting localized blackouts--- after individual line outages. State-of-the-art MILP solvers struggle to tackle the problem's scale and underlying combinatorial explosion, which has motivated the development of expert heuristics, albeit without optimality guarantees.

To overcome these limitations, we propose a scalable Benders decomposition framework for OTSD, combined with several performance improvements. At each iteration, a security analysis identifies the set of sub-problems that yield violated cuts, thus reducing the computational burden of generating Benders cuts. In addition, we improve the granularity of feasibility cuts by introducing additional variables in the master problem. Finally, using graph-theoretic arguments, we present a closed-form formula for separating optimality cuts without solving a linear problem.

The proposed approach is evaluated on test systems representative of medium-sized sub-transmission grids encountered in real-world operations. Numerical experiments demonstrate that, when combined with a local search mechanism, the proposed method improves best-known solutions by up to 30% within minutes of runtime, significantly outperforming standard MILP solvers. These results highlight the potential of decomposition-based techniques to address large-scale OTSD instances and advance decision-making for secure power system operations.

We address the Selective One-to-One Pickup-and-Delivery Problem (SOPDP). The single-vehicle SOPDP, which is NP-hard, requires simultaneously solving two problems: (1) determining how to route an empty delivery vehicle back from its current location to its depot by a scheduled arrival time, and (2) selecting a profit-maximizing subset of spot-market delivery requests along the route subject to the vehicle’s capacity. For finding provably optimal solutions to the multi-vehicle SOPDP, we propose a compact mixed integer programming (MIP) formulation and simple, but effective, preprocessing routines. We also propose a Lagrangian Relaxation Matheuristic (LRM) that takes advantage of multithreading to find high quality solutions within practical time limits. Results from computational experiments on a set of 300 problem instances demonstrate the effectiveness of the LRM. For example, the minimum, median, mean, and maximum running times for Gurobi to solve the MIP for the largest instances are 5 seconds, 56 minutes, 5 hours, and 10 hours respectively. The corresponding times for the LRM are 1, 5, 8, and 30 minutes, respectively; and the minimum, median, and maximum gaps are 0%, 2.98%, and 4.82%, respectively.

We investigate tradeoffs between convergence rates of optimization algorithms and their memory usage for finding a feasible point within a ball of size \epsilon>0 with access to a separation oracle in dimension d. This feasibility problem is a foundational problem in optimization and naturally generalizes first-order convex optimization. Prior works identified three main convergence rate regimes for this problem. In low-accuracy regimes, memory does not play a critical role since gradient descent is already optimal; in moderate-accuracy regimes, prior works showed exponential lower bound tradeoffs between memory and query complexity; but in high-accuracy regimes much less is known. In particular, there are exponential gaps between (1) known memory-constrained algorithms which achieve (d\log(1/\epsilon))^d query complexity and (2) the classical query lower bound d\log(1/\epsilon) by Nemirovski and Yudin. Further, it is unknown whether memory plays any substantial role for solving feasibility problems in this regime. In this work, we show the first memory/query complexity lower bound tradeoff for high-accuracy regimes by proving that any algorithm with less than quadratic memory must have at least d (\log(1/\epsilon))^2 oracle complexity.

Data heterogeneity is a fundamental challenge in modern machine learning. Many gradient-based training algorithms are known to suffer from slow convergence when agents are heterogeneous. This paper is the first to show that, for machine learning problems with separable loss functions and multiple local agents each holding heterogeneous data, this negative effect can turn into an advantage for primal–dual optimization methods. Specifically, while primal-only optimization algorithms (e.g., multi-block PCG, distributed Newton) suffer from slow convergence under heterogeneity, primal–dual algorithms (e.g., EXTRA, consensus multi-block ADMM) can in fact benefit from data heterogeneity. By enforcing consensus constraints and leveraging dual updates, these methods extract additional “momentum” from heterogeneity and accelerate convergence. We further extend our analysis to federated learning, where data heterogeneity is typically harmful to primal-based methods such as FedAvg. In contrast, we demonstrate that heterogeneity can also accelerate primal–dual FL algorithms (e.g., Fed-ADMM). Our results suggest that the central challenge of heterogeneity in multi-agent learning is not to eliminate it, but to exploit it through the design of primal–dual methods, highlighting their potential to turn data heterogeneity from a liability into an advantage in distributed machine learning.

Sparse index tracking aims to replicate the performance of a benchmark index using a small subset of assets and is formulated as a mixed-integer quadratic optimization problem with a cardinality constraint for limiting the number of invested assets. This problem is NP-hard, making it computationally challenging to obtain a globally optimal solution to large-scale problems. Safe screening rules are known as an effective technique for accelerating sparse estimation by pre-eliminating irrelevant variables while ensuring global optimality. However, existing screening rules cannot be directly applied to budget-constrained portfolio optimization, where the sum of investment weights must be one. In this talk, we propose novel safe screening rules for index-tracking problems with both cardinality and budget constraints. We construct a convex perspective relaxation of the problem and derive screening scores by analyzing its Fenchel dual. Based on these scores, we establish theoretically guaranteed conditions for identifying assets that should be included in or excluded from the optimal portfolio. Computational results demonstrate that our screening rules can significantly speed up the solution of large-scale index-tracking problems without compromising global optimality.

Symbolic regression aims to uncover interpretable, closed-form expressions that map input variables to observed outputs. While state-of-the-art solvers—such as Operon and PySR—rely on genetic programming, recent work has approached symbolic regression as a mixed-integer nonlinear program (MINLP) and introduced STAR, a randomized rounding algorithm to obtain high quality solutions of the MINLP. In this work, we develop a novel GPU-accelerated implementation of STAR. We describe the algorithms utilized, including candidate generation, redundancy elimination, and constant optimization. Extensive experiments on SRBench demonstrate that our approach achieves higher symbolic solution rates than most top-performing SRBench methods.

Training Large Language Models (LLMs) typically relies on monolithic preconditioning strategies, where a single update rule—whether based on curvature approximation or spectral normalization—is applied uniformly across all parameters. In this talk, we demonstrate that a Hybrid Preconditioning approach, which applies distinct conditioners to different components of the parameter space, offers a more effective path to convergence.

We illustrate this concept through two recent algorithms: Cosmos and NoMuon.

First, we discuss Cosmos (arXiv:2502.17410), which operates via an orthogonal decoupling of the update direction. Instead of treating the entire space equally, Cosmos applies SOAP (a refined second-order preconditioner) specifically to the Principal Subspace to capture dominant curvature information, while utilizing Muon (a coarse second-order preconditioner) for the Complementary Subspace. This combination naturally balances aggressive curvature exploitation with spectral regularization.

Second, we introduce NoMuon (arXiv:2510.05491), which extends this hybrid philosophy to tensor structures. It applies Muon in the matrix space, while simultaneously utilizing neuron-wise adaptive scaling.

We argue that these combinations are complementary, and leverage the strengths of different conditioners where they are most needed. Our empirical analysis confirms that this hybrid strategy consistently outperforms methods relying on a single type of preconditioning, achieving superior efficiency in training LLMs.

We study linear chance-constrained optimization problems under a continuous-support Wasserstein ambiguity set through the lens of Gaussian mixture model (GMM). We cast the Wasserstein distance over the support of Gaussian distribution parameters using the Bures–Wasserstein metric, rather than the classical Euclidean metric–based form. The resulting GMM-based distributionally robust (DR) chance-constrained program (GMM-DR-CCP) is an infinite-dimensional problem. We establish duality of the worst-case chance satisfaction constraint using the newly cast ambiguity set that allows us to present the infinite-dimensional GMM-DR-CCP equivalently as a finite-dimensional semi-infinite problem. A cutting-surface algorithm is used to solve the semi-infinite problem where a piecewise-linear (PWL) approximation is employed to approximate the standard Gaussian CDF \Phi. While the master problem leverages our earlier results to relate approximation accuracy with the master-problem optimal objective under compact decision space, analogous properties for the subproblem follow from the continuous differentiability of a fractional composite map of \Phi on a compact joint mean–covariance space. The cutting‐surface algorithm is shown to attain any given optimality gap in finite iterations using a finite number of linear segments for PWL approximation. Numerical experiments with electric-vehicle charging station capacity planning and portfolio optimization demonstrate that a continuous-support DR model yields diversified and stable decisions without sacrificing objective values compared with its discrete-support counterpart. The computational result also addresses the limitations of discrete-support models based on the out-of-sample performance analysis, and their bias toward observed or sampled data. The EV application with a charging demand dataset further demonstrates the intrinsic robustness of GMM distribution representation.

In two-stage stochastic programming, it is typically assumed that the behavior of the uncertainty is independent of past decisions. However, in many real-world systems, current choices actively shape the probability of future outcomes. In this work, we relax this assumption by allowing the distribution of second-stage uncertain parameters to depend on first-stage decisions through a general distribution map. We propose an approximation scheme in which this distribution map is replaced by a piecewise-constant approximation over a partition of the first-stage feasible set. Under mild assumptions, we show that the original stochastic program arises as the limit of these approximations as the partition resolution increases. These convergence results establish the theoretical foundation for a practical algorithm that computes approximate solutions to decision-dependent stochastic programs within a prescribed error tolerance.

Cities worldwide face escalating climate risks, such as extreme heat and flooding, demanding innovative tools to proactively design resilient urban environments. Current urban digital twins often focus on single hazards or require extensive real-world data, limiting their use for proactive, multi-hazard scenario planning, especially in data-scarce regions. This project introduces Tessera, a novel digital twin framework that integrates Urban Heat Island (UHI) modeling, flood risk simulation, and synthetic city generation to create a virtual laboratory for climate policy planning. Tessera's core is a hybrid ConvGRU model that predicts spatiotemporal UHI effects by fusing dynamic climate data with static urban features. A key innovation is Tessera’s use of procedurally generated synthetic cities, which enables robust climate impact analysis and the testing of design interventions where detailed city data is unavailable.

The framework’s capability to rapidly simulate multi-hazard climate scenarios is validated by strong predictive performance, achieving a Root Mean Squared Error of 1.44°C for UHI prediction and an AUROC score of 0.73 for flood zone identification. This proven accuracy allows Tessera to serve as a scalable decision-support tool for urban planners and policymakers to stress-test future climate scenarios, representing a critical step towards the proactive, data-driven urban planning needed to navigate the complexities of climate change.

The pooling problem is an NP-hard nonconvex bilinear programming problem that has many industrial applications, including petrochemical refining, wastewater treatment, natural gas transportation, open-pit mining, and animal feed problems. In this talk, we develop a dive-and-fix algorithm which starts with the all zeros trivial solution and then iteratively augments the current solution by involving a new terminal node until all terminal nodes are considered. Two key advantages of the proposed dive-and-fix algorithm that make it particularly suitable for solving large-scale pooling problems are as follows. First, it only requires solving linear programming subproblems, each of which considers a single terminal node. Second, it guarantees to return a feasible solution (usually with high solution quality). As a byproduct of analysis, we prove that checking whether the optimal value of a pooling problem is zero can be accomplished in polynomial time, which closes an open question in the literature. Finally, computational results demonstrate the effectiveness and efficiency of the proposed dive-and-fix algorithm over the direct use of an off-the-shelf global solver and a state-of-the-art greedy construction heuristic algorithm in the literature.

Infinite-horizon risk-sensitive Markov decision processes (MDPs) with discounted costs are difficult to solve because optimal policies may be non-stationary. Standard approaches use state discretization or value function approximation to approximately solve the risk-sensitive problem using (risk-neutral) MDP techniques. For these methods, only limited performance guarantees vis-à-vis explicit stopping criteria or approximation error bounds exist. In this work, we develop variants of value iteration, linear programming, and policy iteration that operate directly on the risk-sensitive MDP. Our analysis establishes that the approximate value functions produced by value iteration and linear programming converge to the optimal risk-sensitive value function, and the resulting policies are provably near-optimal. We further show that the policies generated by our policy-iteration method converge in value to an optimal policy, despite the inherent non-stationarity of risk-sensitive optimal solutions. In addition to asymptotic guarantees, we derive explicit stopping conditions that ensure arbitrarily accurate approximations of the optimal value function.

Modern naval warfare requires autonomous systems capable of persistent maritime reconnaissance and search under uncertainty, resource constraints, and adversarial behavior. We present a safety-constrained, heterogeneous multi-agent reinforcement learning (MARL) framework that jointly optimizes agents’ tactical decisions and refueling logistics under partial observability. The method combines centralized training with decentralized execution and Proximal Policy Optimization, augmented by Lagrangian constrained policy optimization to ensure safety feasibility and mission continuity under strict fuel limitations.

A central feature is an intelligent, adaptive target whose motion induces evolving spatial uncertainty. To capture this, we introduce an Area-of-Uncertainty (AOU) field that updates dynamically with agents’ sensing and target behavior. The AOU provides an information-theoretic reward for reducing uncertainty, enabling coordinated exploration and pursuit driven by information value rather than static search heuristics.

We evaluate the framework in a maritime reconnaissance scenario with stochastic, adaptive target motion, partial observability, and operational constraints. AOU-guided policies achieve high mission completion, robust constraint satisfaction, and efficient fuel and sensing usage. This work advances optimization-informed MARL for dynamic, adversarial naval environments and supports next-generation autonomous systems for contested maritime operations. This research was funded by the Office of Naval Research Mathematical and Resource Optimization program under Grants N000142412280 and N0001425GI01857.

In public school transportation, eligibility policies are often changed in response to community feedback without a full understanding of the impact of these changes on operation costs. Saanich School District in British Columbia, Canada, services nearly 8000 students and currently faces such policy decisions without a clear way to compare costs and benefits of various policies.

This is not an issue unique to Saanich; almost every public school system provides transportation to its students. With limited funding, schools are forced to make difficult budgeting decisions every year, and transportation is a major portion of this annual budget. Each school must therefore choose an eligibility policy that correctly weighs the costs of increasing transportation availability with the benefits to students, their families, and the wider community.

In order to find a better policy, we have begun designing optimization models that maximize the utility to all parties involved while keeping costs below a given threshold. By adding nuance to this model and varying its parameters, we will be able to provide Saanich with a tool to more easily understand the impact on various policy decisions and ultimately decide on the policy that works best for their district.

Prostate cancer represents the second leading cause of cancer mortality in American men, with androgen deprivation therapy (ADT) serving as standard treatment for biochemically recurrent disease. While intermittent ADT (IADT) can reduce cumulative treatment exposure and toxicity compared to continuous therapy, current protocols rely on fixed, heuristic rules derived from clinical trials. These one-size-fits-all approaches fail to account for patient-specific tumor dynamics and mechanisms of treatment resistance. We introduce a novel mixed-integer nonlinear programming (MINLP) framework to derive optimal, personalized IADT schedules that maximize time to progression.

We formulate the optimization problem using a biologically motivated model of prostate-specific antigen (PSA) dynamics that incorporates cancer stem-like cell self-renewal as a mechanism underlying resistance evolution. The model detects clinical progression through a multi-day consecutive PSA increase detection framework using binary indicator variables. We develop both a baseline simulation model replicating the Bruchovsky et al. clinical trial protocol and an optimization model treating key protocol parameters as decision variables rather than fixed values.

Using patient archetypes parameterized with literature-derived cancer stem cell self-renewal rates, we demonstrate the optimization framework's ability to derive patient-specific treatment schedules. Computational experiments suggest that personalized optimization of protocol parameters can extend time to progression compared to fixed clinical trial protocols. This work exemplifies optimization at the interface of clinical data, treatment decisions, and improved patient outcomes in oncology.

Diffusion models have achieved state-of-the-art performance in generative modeling, with their power significantly enhanced by fine-tuning toward task-specific objectives. However, existing fine-tuning methods often lack a principled foundation and offer limited performance guarantees. In this talk, we present a mathematical framework, termed as generative optimization, for understanding guidance-based diffusion model fine-tuning, providing a systematic perspective on its optimization properties and algorithmic design. We abstract task-specific objectives as a reward function and fine-tuned diffusion models aim to maximize the reward by generating solutions. In the offline (one-step) setting, we show that guidance enables high-quality sample generation, achieving optimality akin to off-policy bandit algorithms. In the online (iterative) setting with objective function queries, we establish a strong connection between guided diffusion and numerical optimization.

In this talk, we will present a highly parallel and derivative-free martingale neural network method, based on the probability theory of Varadhan’s martingale formulation of PDEs, to solve Hamilton-Jacobi-Bellman (HJB) equations arising from stochastic optimal control problems (SOCPs), as well as general quasilinear parabolic partial differential equations (PDEs). In both cases, the PDEs are reformulated into a martingale problem such that loss functions will not require the computation of the gradient or Hessian matrix of the PDE solution, and can be computed in parallel in both time and spatial domains. Moreover, the martingale conditions for the PDEs are enforced using a Galerkin method realized with adversarial learning techniques, eliminating the need for direct computation of the conditional expectations associated with the martingale property. For SOCPs, a derivative-free implementation of the maximum principle for optimal controls is also introduced. The numerical results demonstrate the effectiveness and efficiency of the proposed method, which is capable of solving HJB and quasilinear parabolic PDEs accurately and fast in dimensions as high as 10,000.

Discrete diffusion models have emerged as a powerful modeling method for categorical data on finite state spaces. Driven by Continuous-time Markov Chains (CTMC), a pure jump process natively defined on the state space, discrete diffusion models face unique challenges in Reinforcement Learning (RL) post-training due to a lack of trajectory differentiability. In this talk, I will introduce an effective, scalable framework for addressing this challenge and performing RL finetuning with discrete diffusion models in a theoretically-grounded way, leveraging Stochastic Optimal Control (SOC) techniques applied to CTMC. I will discuss two applications of this framework.

In the first part of the talk, I will introduce the Masked Diffusion Neural Sampler (MDNS), an effective neural sampler training strategy for high-dimensional discrete distributions, which achieves state-of-the-art performance in learning sophisticated distributions, such as Ising and Potts models.

In the second part of the talk, I will present Distribution Matching Policy Optimization (DMPO), a novel RL strategy tailored to finetuning diffusion Large Language Models (diffusion LLMs) and grounded in SOC techniques, which significantly outperforms previous baselines on multiple math reasoning and planning benchmarks.

Most models of networked epidemic processes focus exclusively on person-to-person transmission, and existing mitigation strategies typically rely on centralized optimization methods such as spectral optimization or minimizing the basic reproduction number. However, many diseases also spread through underlying infrastructure systems—including transportation networks, water distribution networks, and hospitals—which can significantly influence epidemic dynamics. In this work, we introduce a layered epidemic model that couples a population network with an infrastructure network carrying a pathogen-borne disease. We propose a new metric, the local basic reproduction number (LBRN), which provides sufficient local conditions for predicting outbreaks. Since LBRNs depend only on local information, they naturally enable distributed decision-making. Leveraging this property, we develop an LBRN-based optimization framework for both budget and rate constrained optimal allocation problems with preventive and corrective control actions. Under mild assumptions imposed on the allocation cost functions, we show that these problems can be formulated as convex geometric programs. The nature of the LBRNs allows the resulting optimization to be solved through standard distributed dual-decomposition algorithms. We conclude with numerical experiments demonstrating the scalability of our method and comparing its performance with existing centralized approaches as the network size increases.