Search

We prove that solution operators of elliptic obstacle-type variational inequalities (or, more generally, locally Lipschitz continuous functions possessing certain pointwise-a.e. convexity properties) are Newton differentiable when considered as maps between suitable Lebesgue spaces and equipped with the strong-weak Bouligand differential as a generalized set-valued derivative. It is shown that this Newton differentiability allows to solve optimal control problems with H1-cost terms and one-sided pointwise control constraints by means of a semismooth Newton method. The superlinear convergence of the resulting algorithm is proved in the infinite-dimensional setting and its mesh independence is demonstrated in numerical experiments. We expect that the findings of this paper are also helpful for the design of numerical solution procedures for quasi-variational inequalities and the optimal control of obstacle-type variational problems., Comment: Minor revisions. Published in SIAM Journal on Control and Optimization, Vol. 61, Iss. 3 (2023)

We study a tropical linear regression problem consisting in finding the best approximation of a set of points by a tropical hyperplane. We establish a strong duality theorem, showing that the value of this problem coincides with the maximal radius of a Hilbert's ball included in a tropical polyhedron. We also show that this regression problem is polynomial-time equivalent to mean payoff games. We illustrate our results by solving an inverse problem from auction theory. In this setting, a tropical hyperplane represents the set of equilibrium prices. Tropical linear regression allows us to quantify the distance of a market to the set of equilibria, and infer secret preferences of a decision maker.

In this work we consider stochastic gradient descent (SGD) for solving linear inverse problems in Banach spaces. SGD and its variants have been established as one of the most successful optimisation methods in machine learning, imaging and signal processing, etc. At each iteration SGD uses a single datum, or a small subset of data, resulting in highly scalable methods that are very attractive for large-scale inverse problems. Nonetheless, the theoretical analysis of SGD-based approaches for inverse problems has thus far been largely limited to Euclidean and Hilbert spaces. In this work we present a novel convergence analysis of SGD for linear inverse problems in general Banach spaces: we show the almost sure convergence of the iterates to the minimum norm solution and establish the regularising property for suitable a priori stopping criteria. Numerical results are also presented to illustrate features of the approach.

In this paper, we study optimal control problems of semilinear elliptic and parabolic equations. A tracking cost functional, quadratic in the control and state variables, is considered. No control constraints are imposed. We prove that the corresponding state equations are well posed for controls in L2. However, it is well known that in the L2 framework the mappings involved in the control problem are not Frechet differentiable in general, which makes any analysis of the optimality conditions challenging. Nevertheless, we prove that every L2 optimal control belongs to L∞, and consequently standard optimality conditions are available. The first author was supported by MCIN/AEI/10.13039/501100011033 under research project PID2020-114837GB-I00. The second author was partially supported by the German Research Foundation (DFG) under project grant Wa 3626/3-2.

Convex nonsmooth optimization problems, whose solutions live in very high dimensional spaces, have become ubiquitous. To solve them, the class of first-order algorithms known as proximal splitting algorithms is particularly adequate: they consist of simple operations, handling the terms in the objective function separately. In this overview, we demystify a selection of recent proximal splitting algorithms: we present them within a unified framework, which consists in applying splitting methods for monotone inclusions in primal-dual product spaces, with well-chosen metrics. Along the way, we easily derive new variants of the algorithms and revisit existing convergence results, extending the parameter ranges in several cases. In particular, we emphasize that when the smooth term in the objective function is quadratic, e.g., for least-squares problems, convergence is guaranteed with larger values of the relaxation parameter than previously known. Such larger values are usually beneficial for the convergence speed in practice., To appear in SIAM Review

When an iterative method is applied to solve the linear equation system in interior point methods (IPMs), the attention is usually placed on accelerating their convergence by designing appropriate preconditioners, but the linear solver is applied as a black box with a standard termination criterion which asks for a sufficient reduction of the residual in the linear system. Such an approach often leads to an unnecessary 'oversolving' of linear equations. In this paper, an IPM that relies on an inner termination criterion not based on the residual of the linear system is introduced and analyzed. Moreover, new indicators for the early termination of the inner iterations are derived from a deep understanding of IPM needs. The new technique has been adapted to the Conjugate Gradient (CG) and to the Minimum Residual method (MINRES) applied in the IPM context. The new criterion has been tested on a set of quadratic optimization problems including compressed sensing, image processing and instances with partial differential equation constraints, and it has been compared to standard residual tests with variable tolerance. Evidence gathered from these computational experiments shows that the new technique delivers significant improvements in terms of inner (linear) iterations and those translate into significant savings of the IPM solution time.

Optimal control problem is typically solved by first finding the value function through Hamilton-Jacobi equation (HJE) and then taking the minimizer of the Hamiltonian to obtain the control. In this work, instead of focusing on the value function, we propose a new formulation for the gradient of the value function (value-gradient) as a decoupled system of partial differential equations in the context of continuous-time deterministic discounted optimal control problem. We develop an efficient iterative scheme for this system of equations in parallel by utilizing the properties that they share the same characteristic curves as the HJE for the value function. For the theoretical part, we prove that this iterative scheme converges linearly in $L_\alpha^2$ sense for some suitable exponent $\alpha$ in a weight function. For the numerical method, we combine characteristic line method with machine learning techniques. Specifically, we generate multiple characteristic curves at each policy iteration from an ensemble of initial states, and compute both the value function and its gradient simultaneously on each curve as the labelled data. Then supervised machine learning is applied to minimize the weighted squared loss for both the value function and its gradients. Experimental results demonstrate that this new method not only significantly increases the accuracy but also improves the efficiency and robustness of the numerical estimates, particularly with less amount of characteristics data or fewer training steps.

Given a parameter dependent fixed point equation $x = F(x,u)$, we derive an abstract compactness principle for the fixed point map $u \mapsto x^*(u)$ under the assumptions that (i) the fixed point equation can be solved by the contraction principle and (ii) the map $u \mapsto F(x,u)$ is compact for fixed $x$. This result is applied to infinite-dimensional, semi-linear control systems and their reachable sets. More precisely, we extend a non-controllability result of Ball, Marsden, and Slemrod [1] to semi-linear systems. First we consider $L^p$-controls, $p>1$. Subsequently we analyze the case $p=1$.

Necessary conditions for asymptotic stability and stabilizability of subsets for dynamical and control systems are obtained. The main necessary condition is homotopical and is in turn used to obtain a homological one. A certain extension is ruled out. Questions are posed., Comment: Accepted to SIAM J Control and Optimization

A hyperbolicity cone is said to be rank-one generated (ROG) if all its extreme rays have rank one, where the rank is computed with respect to the underlying hyperbolic polynomial. This is a natural class of hyperbolicity cones which are strictly more general than the ROG spectrahedral cones. In this work, we present a study of the automorphisms of ROG hyperbolicity cones and their derivative relaxations. One of our main results states that the automorphisms of the derivative relaxations are exactly the automorphisms of the original cone fixing a certain direction. As an application, we completely determine the automorphisms of the derivative relaxations of the nonnegative orthant and of the cone of positive semidefinite matrices. More generally, we also prove relations between the automorphisms of a spectral cone and the underlying permutation-invariant set, which might be of independent interest., Comment: 25 pages. Some minor fixes and changes. To appear at the SIAM Journal on Applied Algebra and Geometry

In this paper, we consider an optimal reinsurance problem to minimize the probability of drawdown for the scaled Cram\'er-Lundberg risk model when the reinsurance premium is computed according to the mean-variance premium principle. We extend the work of Liang et al. [16] to the case of minimizing the probability of drawdown. By using the comparison method and the tool of adjustment coefficients, we show that the minimum probability of drawdown for the scaled classical risk model converges to the minimum probability for its diffusion approximation, and the rate of convergence is of order $O(n^{-1/2})$. We further show that using the optimal strategy from the diffusion approximation in the scaled classical risk model is $O(n^{-1/2})$-optimal.

The absolute value equations (AVE) problem is an algebraic problem of solving Ax+|x|=b. So far, most of the research focused on methods for solving AVEs, but we address the problem itself by analysing properties of AVE and the corresponding solution set. In particular, we investigate topological properties of the solution set, such as convexity, boundedness, connectedness, or whether it consists of finitely many solutions. Further, we address problems related to nonnegativity of solutions such as solvability or unique solvability. AVE can be formulated by means of different optimization problems, and in this regard we are interested in how the solutions of AVE are related with optima, Karush-Kuhn-Tucker points and feasible solutions of these optimization problems. We characterize the matrix classes associated with the above mentioned properties and inspect the computational complexity of the recognition problem; some of the classes are polynomially recognizable, but some others are proved to be NP-hard. For the intractable cases, we propose various sufficient conditions. We also post new challenging problems that raised during the investigation of the problem.

The paper is concerned with a zero-sum Stackelberg stochastic linear-quadratic (LQ, for short) differential game over finite horizons. Under a fairly weak condition, the Stackelberg equilibrium is explicitly obtained by first solving a forward stochastic LQ optimal control problem (SLQ problem, for short) and then a backward SLQ problem. Two Riccati equations are derived in constructing the Stackelberg equilibrium. An interesting finding is that the difference of these two Riccati equations coincides with the Riccati equation associated with the zero-sum Nash stochastic LQ differential game, which implies that the Stackelberg equilibrium and the Nash equilibrium are actually identical. Consequently, the Stackelberg equilibrium admits a linear state feedback representation, and the Nash game can be solved in a leader-follower manner.

We study non-autonomous observation systems \begin{align*} \dot{x}(t) = A(t) x(t),\quad y(t) = C(t) x(t),\quad x(0) = x_0\in X, \end{align*} where $(A(t))$ is a strongly measurable family of closed operators on a Banach space $X$ and $(C(t))$ is a family of bounded observation operators from $X$ to a Banach space $Y$. Based on an abstract uncertainty principle and a dissipation estimate, we prove that the observation system satisfies a final-state observability estimate in $\mathrm{L}^r(E; Y)$ for measurable subsets $E \subseteq [0,T], T > 0$. We present applications of the above result to families $(A(t))$ of uniformly strongly elliptic differential operators as well as non-autonomous Ornstein-Uhlenbeck operators $P(t)$ on $\mathrm{L}^p(\mathbb{R}^d)$ with observation operators $C(t)u = u|_{\Omega(t)}$. In the setting of non-autonomous strongly elliptic operators, we derive necessary and sufficient geometric conditions on the family of sets $(\Omega(t))$ such that the corresponding observation system satisfies a final-state observability estimate., Comment: minor corrections

The objective of this paper is to analyze the existence of equilibria for a class of deterministic mean field games of controls. The interaction between players is due to both a congestion term and a price function which depends on the distributions of the optimal strategies. Moreover, final state and mixed state-control constraints are considered, the dynamics being nonlinear and affine with respect to the control. The existence of equilibria is obtained by Kakutani's theorem, applied to a fixed point formulation of the problem. Finally, uniqueness results are shown under monotonicity assumptions.

We present the concept of a Generalized Feedback Nash Equilibrium (GFNE) in dynamic games, extending the Feedback Nash Equilibrium concept to games in which players are subject to state and input constraints. We formalize necessary and sufficient conditions for (local) GFNE solutions at the trajectory level, which enable the development of efficient numerical methods for their computation. Specifically, we propose a Newton-style method for finding game trajectories which satisfy necessary conditions for an equilibrium, which can then be checked against sufficiency conditions. We show that the evaluation of the necessary conditions in general requires computing a series of nested, implicitly-defined derivatives, which quickly becomes intractable. To this end, we introduce an approximation to the necessary conditions which is amenable to efficient evaluation, and in turn, computation of solutions. We term the solutions to the approximate necessary conditions Generalized Feedback Quasi-Nash Equilibria (GFQNE), and we introduce numerical methods for their computation. In particular, we develop a Sequential Linear-Quadratic Game approach, in which a LQ local approximation of the game is solved at each iteration. The development of this method relies on the ability to compute a GFNE to inequality- and equality-constrained LQ games, and therefore specific methods for the solution of these special cases are developed in detail. We demonstrate the effectiveness of the proposed solution approach on a dynamic game arising in an autonomous driving application.

This paper presents a computational approach for finding the optimal shapes of peristaltic pumps transporting rigid particles in Stokes flow. In particular, we consider shapes that minimize the rate of energy dissipation while pumping a prescribed volume of fluid, number of particles and/or distance traversed by the particles over a set time period. Our approach relies on a recently developed fast and accurate boundary integral solver for simulating multiphase flows through periodic geometries of arbitrary shapes. In order to fully capitalize on the dimensionality reduction feature of the boundary integral methods, shape sensitivities must ideally involve evaluating the physical variables on the particle or pump boundaries only. We show that this can indeed be accomplished owing to the linearity of Stokes flow. The forward problem solves for the particle motion in a slip-driven pipe flow while the adjoint problems in our construction solve quasi-static Dirichlet boundary value problems backwards in time, retracing the particle evolution. The shape sensitivities simply depend on the solution of one forward and one adjoint (for each shape functional) problems. We validate these analytic shape derivative formulas by comparing against finite-difference based gradients and present several examples showcasing optimal pump shapes under various constraints., Comment: 24 pages, 8 figures

In this paper we study and compare two multigrid relaxation schemes with coarsening by two, three, and four for solving elliptic sparse optimal control problems with control constraints. First, we perform a detailed local Fourier analysis (LFA) of a well-known collective Jacobi relaxation (CJR) scheme, where the optimal smoothing factors are derived. This insightful analysis reveals that the optimal relaxation parameters depend on mesh size and regularization parameters, which was not investigated in literature. Second, we propose and analyze a new mass-based Braess-Sarazin relaxation (BSR) scheme, which is proven to provide smaller smoothing factors than the CJR scheme when $\alpha\ge ch^4$ for a small constant $c$. Here $\alpha$ is the regularization parameter and $h$ is the spatial mesh step size. These schemes are successfully extended to control-constrained cases through the semi-smooth Newton method. Coarsening by three or four with BSR is competitive with coarsening by two. Numerical examples are presented to validate our theoretical outcomes. The proposed inexact BSR (IBSR) scheme, where two preconditioned conjugate gradients iterations are applied to the Schur complement system, yields better computational efficiency than the CJR scheme., Comment: 23 pages; 2 tables and 5 figures

In this paper, we extend a recently established subgradient method for the computation of Riemannian metrics that optimizes certain singular value functions associated with dynamical systems. This extension is threefold. First, we introduce a projected subgradient method which results in Riemannian metrics whose parameters are confined to a compact convex set and we can thus prove that a minimizer exists; second, we allow inexact subgradients and study the effect of the errors on the computed metrics; and third, we analyze the subgradient algorithm for three different choices of step sizes: constant, exogenous and Polyak. The new methods are illustrated by application to dimension and entropy estimation of the H\'enon map., Comment: 30 pages, 5 figures

Solving high-dimensional optimal control problems in real-time is an important but challenging problem, with applications to multi-agent path planning problems, which have drawn increased attention given the growing popularity of drones in recent years. In this paper, we propose a novel neural network method called SympOCnet that applies the Symplectic network to solve high-dimensional optimal control problems with state constraints. We present several numerical results on path planning problems in two-dimensional and three-dimensional spaces. Specifically, we demonstrate that our SympOCnet can solve a problem with more than 500 dimensions in 1.5 hours on a single GPU, which shows the effectiveness and efficiency of SympOCnet. The proposed method is scalable and has the potential to solve truly high-dimensional path planning problems in real-time.

Optimal Transport, a theory for optimal allocation of resources, is widely used in various fields such as astrophysics, machine learning, and imaging science. However, many applications impose elementwise constraints on the transport plan which traditional optimal transport cannot enforce. Here we introduce Supervised Optimal Transport (sOT) that formulates a constrained optimal transport problem where couplings between certain elements are prohibited according to specific applications. sOT is proved to be equivalent to an $l^1$ penalized optimization problem, from which efficient algorithms are designed to solve its entropy regularized formulation. We demonstrate the capability of sOT by comparing it to other variants and extensions of traditional OT in color transfer problem. We also study the barycenter problem in sOT formulation, where we discover and prove a unique reverse and portion selection (control) mechanism. Supervised optimal transport is broadly applicable to applications in which constrained transport plan is involved and the original unit should be preserved by avoiding normalization.

We study the problem of computing weighted sum-of-squares (WSOS) certificates for positive polynomials over a compact semialgebraic set. Building on the theory of interior-point methods for convex optimization, we introduce the concept of dual cone certificates, which allows us to interpret vectors from the dual of the sum-of-squares cone as rigorous nonnegativity certificates of a WSOS polynomial. Whereas conventional WSOS certificates are alternative representations of the polynomials they certify, dual certificates are distinct from the certified polynomials; moreover, each dual certificate certifies a full-dimensional convex cone of WSOS polynomials. As a result, rational WSOS certificates can be constructed from numerically computed dual certificates at little additional cost, without any rounding or projection steps applied to the numerical certificates. As an additional algorithmic application, we present an almost entirely numerical hybrid algorithm for computing the optimal WSOS lower bound of a given polynomial along with a rational dual certificate, with a polynomial-time computational cost per iteration and linear rate of convergence., Comment: Preprint accepted in SIAM Journal on Optimization. Major revision compared to v1

In this paper, we analyze and discuss the well-posedness of a new variant of the so-called sweeping process, introduced by J.J. Moreau in the early 70's \cite{More71} with motivation in plasticity theory. In this variant, the normal cone to the (mildly non-convex) prox-regular moving set $C(t)$, supposed to have an absolutely continuous variation, is perturbed by a sum of a Carath\'{e}odory mapping and an integral forcing term. The integrand of the forcing term depends on two time-variables, that is, we study a general integro-differential sweeping process of Volterra type. By setting up an appropriate semi-discretization method combined with a new Gronwall-like inequality (differential inequality), we show that the integro-differential sweeping process has one and only one absolutely continuous solution. We also establish the continuity of the solution with respect to the initial value. The results of the paper are applied to the study of nonlinear integro-differential complementarity systems which are combination of Volterra integro-differential equations with nonlinear complementarity constraints. Another application is concerned with non-regular electrical circuits containing time-varying capacitors and nonsmooth electronic device like diodes. Both applications represent an additional novelty of our paper., Comment: 36 pages, 1 figures

This work considers variational Bayesian inference as an inexpensive and scalable alternative to a fully Bayesian approach in the context of sparsity-promoting priors. In particular, the priors considered arise from scale mixtures of Normal distributions with a generalized inverse Gaussian mixing distribution. This includes the variational Bayesian LASSO as an inexpensive and scalable alternative to the Bayesian LASSO introduced in [65]. It also includes a family of priors which more strongly promote sparsity. For linear models the method requires only the iterative solution of deterministic least squares problems. Furthermore, for p unknown covariates the method can be implemented exactly online with a cost of $O(p^3)$ in computation and $O(p^2)$ in memory per iteration -- in other words, the cost per iteration is independent of n, and in principle infinite data can be considered. For large $p$ an approximation is able to achieve promising results for a cost of $O(p)$ per iteration, in both computation and memory. Strategies for hyper-parameter tuning are also considered. The method is implemented for real and simulated data. It is shown that the performance in terms of variable selection and uncertainty quantification of the variational Bayesian LASSO can be comparable to the Bayesian LASSO for problems which are tractable with that method, and for a fraction of the cost. The present method comfortably handles $n = 65536$, $p = 131073$ on a laptop in less than 30 minutes, and $n = 10^5$, $p = 2.1 \times 10^6$ overnight.

We prove that the viscosity solution to a Hamilton-Jacobi equation with a smooth convex Hamiltonian of the form $H(x,p)$ is differentiable with respect to the initial condition. Moreover, the directional G��teaux derivatives can be explicitly computed almost everywhere in $\mathbb{R}^N$ by means of the optimality system of the associated optimal control problem. We also prove that, in the one-dimensional case in space and in the quadratic case in any space dimension, these directional G��teaux derivatives actually correspond to the unique duality solution to the linear transport equation with discontinuous coefficient, resulting from the linearization of the Hamilton-Jacobi equation. The motivation behind these differentiability results arises from the following optimal inverse-design problem: given a time horizon $T>0$ and a target function $u_T$, construct an initial condition such that the corresponding viscosity solution at time $T$ minimizes the $L^2$-distance to $u_T$. Our differentiability results allow us to derive a necessary first-order optimality condition for this optimization problem, and the implementation of gradient-based methods to numerically approximate the optimal inverse design., 37 pages, 2 figures

Based on the stochastic maximum principle for the partially coupled forward-backward stochastic control system (FBSCS for short), a modified method of successive approximations (MSA for short) is established for stochastic recursive optimal control problems. The second-order adjoint processes are introduced in the augmented Hamiltonian minimization step since the control domain is not necessarily convex. Thanks to the theory of bounded mean oscillation martingales (BMO martingales for short), we give a delicate proof of the error estimate and then prove the convergence of the modified MSA algorithm. In a special case, we obtain a logarithmic convergence rate. When the control domain is convex and compact, a sufficient condition which makes the control returned from the MSA algorithm be a near-optimal control is given for a class of linear FBSCSs., Comment: 30 pages

We extend the method of controlled Lagrangians with kinetic shaping to those mechanical systems on semidirect product Lie groups with broken symmetry, more specifically to the Euler--Poincar\'e equations with advected parameters. We find a matching condition for the controlled Lagrangian for such systems whose configuration manifold is a general semidirect product Lie group $\mathsf{G} \ltimes V$. Our motivating examples are a bottom-heavy underwater vehicle and a top spinning on a movable base. Their configuration space is the special Euclidean group $\mathsf{SE}(3) = \mathsf{SO}(3) \ltimes \mathbb{R}^{3}$, where the $\mathsf{SE}(3)$-symmetry is broken by the gravity. The controls resulting from the matching condition stabilize unstable equilibria of these examples. Furthermore, the matching helps us find additional dissipative controls that asymptotically stabilize those unstable equilibria., Comment: 29 pages, 10 figures

Spanning trees are a representative example of linear matroid bases that are efficiently countable. Perfect matchings of Pfaffian bipartite graphs are a countable example of common bases of two matrices. Generalizing these two examples, Webb (2004) introduced the notion of Pfaffian pairs as a pair of matrices for which counting of their common bases is tractable via the Cauchy-Binet formula. This paper studies counting on linear matroid problems extending Webb's work. We first introduce "Pfaffian parities" as an extension of Pfaffian pairs to the linear matroid parity problem, which is a common generalization of the linear matroid intersection problem and the matching problem. We enumerate combinatorial examples of Pfaffian pairs and parities. The variety of the examples illustrates that Pfaffian pairs and parities serve as a unified framework of efficiently countable discrete structures. Based on this framework, we derive celebrated counting theorems, such as Kirchhoff's matrix-tree theorem, Tutte's directed matrix-tree theorem, the Pfaffian matrix-tree theorem, and the Lindstr\"om-Gessel-Viennot lemma. Our study then turns to algorithmic aspects. We observe that the fastest randomized algorithms for the linear matroid intersection and parity problems by Harvey (2009) and Cheung-Lau-Leung (2014) can be derandomized for Pfaffian pairs and parities. We further present polynomial-time algorithms to count the number of minimum-weight solutions on weighted Pfaffian pairs and parities. Our algorithms make use of Frank's weight splitting lemma for the weighted matroid intersection problem and the algebraic optimality criterion of the weighted linear matroid parity problem given by Iwata-Kobayashi (2017).

We study the lattice width of lattice-free polyhedra given by $\mathbf{A}\mathbf{x}\leq\mathbf{b}$ in terms of $\Delta(\mathbf{A})$, the maximal $n\times n$ minor in absolute value of $\mathbf{A}\in\mathbb{Z}^{m\times n}$. Our main contribution is to link the lattice width of lattice-free polyhedra to the height of Hilbert bases and to the diameter of finite abelian groups. This leads to a bound on the lattice width of lattice-free pyramids which solely depends on $\Delta(\mathbf{A})$ provided a conjecture regarding the height of Hilbert bases holds. Further, we exploit a combination of techniques to obtain novel bounds on the lattice width of simplices. A second part of the paper is devoted to a study of the above mentioned Hilbert basis conjecture. We give a complete characterization of the Hilbert basis if $\Delta(\mathbf{A}) = 2$ which implies the conjecture in that case and prove its validity for simplicial cones.

We generalize the well-known primal-dual algorithm proposed by Chambolle and Pock for saddle point problems, and improve the condition for ensuring its convergence. The improved convergence-guaranteeing condition is effective for the generic setting, and it is shown to be optimal. It also allows us to discern larger step sizes for the resulting subproblems, and thus provides a simple and universal way to improve numerical performance of the original primal-dual algorithm. In addition, we present a structure-exploring heuristic to further relax the convergence-guaranteeing condition for some specific saddle point problems, which could yield much larger step sizes and hence significantly better performance. Effectiveness of this heuristic is numerically illustrated by the classic assignment problem.

A Dirichlet $k$-partition of a domain is a collection of $k$ pairwise disjoint open subsets such that the sum of their first Laplace--Dirichlet eigenvalues is minimal. In this paper, we propose a new relaxation of the problem by introducing auxiliary indicator functions of domains and develop a simple and efficient diffusion generated method to compute Dirichlet $k$-partitions for arbitrary domains. The method only alternates three steps: 1. convolution, 2. thresholding, and 3. projection. The method is simple, easy to implement, insensitive to initial guesses and can be effectively applied to arbitrary domains without any special discretization. At each iteration, the computational complexity is linear in the discretization of the computational domain. Moreover, we theoretically prove the energy decaying property of the method. Experiments are performed to show the accuracy of approximation, efficiency and unconditional stability of the algorithm. We apply the proposed algorithms on both 2- and 3-dimensional flat tori, triangle, square, pentagon, hexagon, disk, three-fold star, five-fold star, cube, ball, and tetrahedron domains to compute Dirichlet $k$-partitions for different $k$ to show the effectiveness of the proposed method. Compared to previous work with reported computational time, the proposed method achieves hundreds of times acceleration., 24 pages, 17 figures

Material response of real, passive, linear, time-invariant media to external influences is described by complex analytic functions of frequency that can always be written in terms of Stieltjes functions -- a special class of analytic functions mapping complex upper half-plane into itself. Reconstructing such functions from their experimentally measured values at specific frequencies is one of the central problems that we address in this paper. A definitive reconstruction algorithm that produces a certificate of optimality as well as a graphical representation of the uncertainty of reconstruction is proposed. Its effectiveness is demonstrated in the context of the electrochemical impedance spectroscopy., 37 pages, 6 figures

In this paper we consider impulse control of continuous time Markov processes with average cost per unit time functional. This problem is approximated using impulse control problems stopped at the first exit time from increasing sequence of open sets. We find solution to Bellman equation corresponding to the original problem and show that stopped impulse control problems approximate optimal value of the cost functional.

We study a general convex optimization problem, which covers various classic problems in different areas and particularly includes many optimal transport related problems arising in recent years. To solve this problem, we revisit the classic Bregman proximal point algorithm (BPPA) and introduce a new inexact stopping condition for solving the subproblems, which can circumvent the underlying feasibility difficulty often appearing in existing inexact conditions when the problem has a complex feasible set. Our inexact condition also covers several existing inexact conditions as special cases and hence makes our inexact BPPA (iBPPA) more flexible to fit different scenarios in practice. Moreover, inspired by Nesterov's acceleration technique, we develop an inertial variant of our iBPPA, denoted by V-iBPPA, and establish the iteration complexity of $O(1/k^{\lambda})$, where $\lambda\geq1$ is a quadrangle scaling exponent of the kernel function. In particular, when the proximal parameter is a constant and the kernel function is strongly convex with Lipschitz continuous gradient (hence $\lambda=2$), our V-iBPPA achieves a faster rate of $O(1/k^2)$ just as existing accelerated inexact proximal point algorithms. Some preliminary numerical experiments for solving the standard OT problem are conducted to show the convergence behaviors of our iBPPA and V-iBPPA under different inexactness settings. The experiments also empirically verify the potential of our V-iBPPA on improving the convergence speed.

We present several characterizations, via some weak observability inequalities, on the complete stabilizability for a control system $[A,B]$, i.e., $y'(t)=Ay(t)+Bu(t)$, $t\geq 0$, where $A$ generates a $C_0$-semigroup on a Hilbert space $X$ and $B$ is a linear and bounded operator from another Hilbert space $U$ to $X$. We then extend the aforementioned characterizations in two directions: first, the control operator $B$ is unbounded; second, the control system is time-periodic. We also give some sufficient conditions, from the perspective of the spectral projections, to ensure the weak observability inequalities. As applications, we provide several examples, which are not null controllable, but can be verified, via the weak observability inequalities, to be completely stabilizable., 28 pages

This paper studies an infinite horizon optimal control problem for discrete-time linear systems and quadratic criteria, both with random parameters which are independent and identically distributed with respect to time. A classical approach is to solve an algebraic Riccati equation that involves mathematical expectations and requires certain statistical information of the parameters. In this paper, we propose an online iterative algorithm in the spirit of Q-learning for the situation where only one random sample of parameters emerges at each time step. The first theorem proves the equivalence of three properties: the convergence of the learning sequence, the well-posedness of the control problem, and the solvability of the algebraic Riccati equation. The second theorem shows that the adaptive feedback control in terms of the learning sequence stabilizes the system as long as the control problem is well-posed. Numerical examples are presented to illustrate our results., Comment: 24 pages, 3 figures

We consider a pointwise tracking optimal control problem for a semilinear elliptic partial differential equation. We derive the existence of optimal solutions and analyze first and, necessary and sufficient, second order optimality conditions. We devise two strategies of discretization to approximate a solution of the optimal control problem: a semidiscrete scheme where the control variable is not discretized -- the so-called variational discretization approach -- and a fully discrete scheme where the control variable is discretized with piecewise constant functions. For both solution techniques, we analyze convergence properties of discretizations and derive error estimates.

We introduce a simple, accurate, and extremely efficient method for numerically solving the multi-marginal optimal transport (MMOT) problems arising in density functional theory. The method relies on (i) the sparsity of optimal plans [for $N$ marginals discretized by $\ell$ gridpoints each, general Kantorovich plans require $\ell^N$ gridpoints but the support of optimizers is of size $O(\ell\cdot N)$ [FV18]], (ii) the method of column generation (CG) from discrete optimization which to our knowledge has not hitherto been used in MMOT, and (iii) ideas from machine learning. The well-known bottleneck in CG consists in generating new candidate columns efficiently; we prove that in our context, finding the best new column is an NP-complete problem. To overcome this bottleneck we use a genetic learning method tailormade for MMOT in which the dual state within CG plays the role of an "adversary", in loose similarity to Wasserstein GANs. On a sequence of benchmark problems with up to 120 gridpoints and up to 30 marginals, our method always found the exact optimizers. Moreover, empirically the number of computational steps needed to find them appears to scale only polynomially when both $N$ and $\ell$ are simultaneously increased (while keeping their ratio fixed to mimic a thermodynamic limit of the particle system).

Controlling systems of ordinary differential equations (ODEs) is ubiquitous in science and engineering. For finding an optimal feedback controller, the value function and associated fundamental equations such as the Bellman equation and the Hamilton-Jacobi-Bellman (HJB) equation are essential. The numerical treatment of these equations poses formidable challenges due to their non-linearity and their (possibly) high-dimensionality. In this paper we consider a finite horizon control system with associated Bellman equation. After a time-discretization, we obtain a sequence of short time horizon problems which we call local optimal control problems. For solving the local optimal control problems we apply two different methods, one being the well-known policy iteration, where a fixed-point iteration is required for every time step. The other algorithm borrows ideas from Model Predictive Control (MPC), by solving the local optimal control problem via open-loop control methods on a short time horizon, allowing us to replace the fixed-point iteration by an adjoint method. For high-dimensional systems we apply low rank hierarchical tensor product approximation/tree-based tensor formats, in particular tensor trains (TT tensors) and multi-polynomials, together with high-dimensional quadrature, e.g. Monte-Carlo. We prove a linear error propagation with respect to the time discretization and give numerical evidence by controlling a diffusion equation with unstable reaction term and an Allen-Kahn equation.

We investigate the feasibility of sample average approximation (SAA) for general stochastic optimization problems, including two-stage stochastic programming without the relatively complete recourse assumption. Instead of analyzing problems with specific structures, we utilize results from the Vapnik-Chervonenkis (VC) dimension and Probably Approximately Correct learning to provide a general framework that offers explicit feasibility bounds for SAA solutions under minimal structural or distributional assumption. We show that, as long as the hypothesis class formed by the feasbible region has a finite VC dimension, the infeasibility of SAA solutions decreases exponentially with computable rates and explicitly identifiable accompanying constants. We demonstrate how our bounds apply more generally and competitively compared to existing results.

It was shown recently by Su et al. (2016) that Nesterov's accelerated gradient method for minimizing a smooth convex function $f$ can be thought of as the time discretization of a second-order ODE, and that $f(x(t))$ converges to its optimal value at a rate of $\mathcal{O}(1/t^2)$ along any trajectory $x(t)$ of this ODE. A variational formulation was introduced in Wibisono et al. (2016) which allowed for accelerated convergence at a rate of $\mathcal{O}(1/t^p)$, for arbitrary $p>0$, in normed vector spaces. This framework was exploited in Duruisseaux et al. (2021) to design efficient explicit algorithms for symplectic accelerated optimization. In Alimisis et al. (2020), a second-order ODE was proposed as the continuous-time limit of a Riemannian accelerated algorithm, and it was shown that the objective function $f(x(t))$ converges to its optimal value at a rate of $\mathcal{O}(1/t^2)$ along solutions of this ODE. In this paper, we show that on Riemannian manifolds, the convergence rate of $f(x(t))$ to its optimal value can also be accelerated to an arbitrary convergence rate $\mathcal{O}(1/t^p)$, by considering a family of time-dependent Bregman Lagrangian and Hamiltonian systems on Riemannian manifolds. This generalizes the results of Wibisono et al. (2016) to Riemannian manifolds and also provides a variational framework for accelerated optimization on Riemannian manifolds. An approach based on the time-invariance property of the family of Bregman Lagrangians and Hamiltonians was used to construct very efficient optimization algorithms in Duruisseaux et al. (2021), and we establish a similar time-invariance property in the Riemannian setting. One expects that a geometric numerical integrator that is time-adaptive, symplectic, and Riemannian manifold preserving will yield a class of promising optimization algorithms on manifolds., 23 pages, 3 figures

We provide a proximal average with repect to a $1$-coercive Legendre function. In the sense of Bregman distance, the Bregman envelope of the proximal average is a convex combination of Bregman envelopes of individual functions. The Bregman proximal mapping of the average is a convex combination of convexified proximal mappings of individual functions. Techniques from variational analysis provide the keys for the Bregman proximal average.

We consider an elliptic linear-quadratic parameter estimation problem with a finite number of parameters. A novel a priori bound for the parameter error is proved and, based on this bound, an adaptive finite element method driven by an a posteriori error estimator is presented. Unlike prior results in the literature, our estimator, which is composed of standard energy error residual estimators for the state equation and suitable co-state problems, reflects the faster convergence of the parameter error compared to the (co)-state variables. We show optimal convergence rates of our method; in particular and unlike prior works, we prove that the estimator decreases with a rate that is the sum of the best approximation rates of the state and co-state variables. Experiments confirm that our method matches the convergence rate of the parameter error.

We study solution sensitivity for nonlinear programs (NLPs) whose structures are induced by graphs. These NLPs arise in many applications such as dynamic optimization, stochastic optimization, optimization with partial differential equations, and network optimization. We show that for a given pair of nodes, the sensitivity of the primal-dual solution at one node against a data perturbation at the other node decays exponentially with respect to the distance between these two nodes on the graph. In other words, the solution sensitivity decays as one moves away from the perturbation point. This result, which we call exponential decay of sensitivity, holds under the strong second-order sufficiency condition and the linear independence constraint qualification. We also present conditions under which the decay rate remains uniformly bounded; this allows us to characterize the sensitivity behavior of NLPs defined over subgraphs of infinite graphs. The theoretical developments are illustrated with numerical examples.

The focus of this paper is on stochastic variational inequalities (VI) under Markovian noise. A prominent application of our algorithmic developments is the stochastic policy evaluation problem in reinforcement learning. Prior investigations in the literature focused on temporal difference (TD) learning by employing nonsmooth finite time analysis motivated by stochastic subgradient descent leading to certain limitations. These encompass the requirement of analyzing a modified TD algorithm that involves projection to an a-priori defined Euclidean ball, achieving a non-optimal convergence rate and no clear way of deriving the beneficial effects of parallel implementation. Our approach remedies these shortcomings in the broader context of stochastic VIs and in particular when it comes to stochastic policy evaluation. We developed a variety of simple TD learning type algorithms motivated by its original version that maintain its simplicity, while offering distinct advantages from a non-asymptotic analysis point of view. We first provide an improved analysis of the standard TD algorithm that can benefit from parallel implementation. Then we present versions of a conditional TD algorithm (CTD), that involves periodic updates of the stochastic iterates, which reduce the bias and therefore exhibit improved iteration complexity. This brings us to the fast TD (FTD) algorithm which combines elements of CTD and the stochastic operator extrapolation method of the companion paper. For a novel index resetting policy FTD exhibits the best known convergence rate. We also devised a robust version of the algorithm that is particularly suitable for discounting factors close to 1., Comment: arXiv admin note: text overlap with arXiv:2011.02987

We propose new Riemannian preconditioned algorithms for low-rank tensor completion via the polyadic decomposition of a tensor. These algorithms exploit a non-Euclidean metric on the product space of the factor matrices of the low-rank tensor in the polyadic decomposition form. This new metric is designed using an approximation of the diagonal blocks of the Hessian of the tensor completion cost function, thus has a preconditioning effect on these algorithms. We prove that the proposed Riemannian gradient descent algorithm globally converges to a stationary point of the tensor completion problem, with convergence rate estimates using the $\L{}$ojasiewicz property. Numerical results on synthetic and real-world data suggest that the proposed algorithms are more efficient in memory and time compared to state-of-the-art algorithms. Moreover, the proposed algorithms display a greater tolerance for overestimated rank parameters in terms of the tensor recovery performance, thus enable a flexible choice of the rank parameter., Comment: 27 Pages, 10 figures, 4 tables

This paper is concerned with the column $\ell_{2,0}$-regularized factorization model of low-rank matrix recovery problems and its computation. The column $\ell_{2,0}$-norm of factor matrices is introduced to promote column sparsity of factors and low-rank solutions. For this nonconvex discontinuous optimization problem, we develop an alternating majorization-minimization (AMM) method with extrapolation, and a hybrid AMM in which a majorized alternating proximal method is proposed to seek an initial factor pair with less nonzero columns and the AMM with extrapolation is then employed to minimize of a smooth nonconvex loss. We provide the global convergence analysis for the proposed AMM methods and apply them to the matrix completion problem with non-uniform sampling schemes. Numerical experiments are conducted with synthetic and real data examples, and comparison results with the nuclear-norm regularized factorization model and the max-norm regularized convex model show that the column $\ell_{2,0}$-regularized factorization model has an advantage in offering solutions of lower error and rank within less time.

Chance constraints provide a principled framework to mitigate the risk of high-impact extreme events by modifying the controllable properties of a system. The low probability and rare occurrence of such events, however, impose severe sampling and computational requirements on classical solution methods that render them impractical. This work proposes a novel sampling-free method for solving rare chance constrained optimization problems affected by uncertainties that follow general Gaussian mixture distributions. By integrating modern developments in large deviation theory with tools from convex analysis and bilevel optimization, we propose tractable formulations that can be solved by off-the-shelf solvers. Our formulations enjoy several advantages compared to classical methods: their size and complexity is independent of event rarity, they do not require linearity or convexity assumptions on system constraints, and under easily verifiable conditions, serve as safe conservative approximations or asymptotically exact reformulations of the true problem. Computational experiments on linear, nonlinear and PDE-constrained problems from applications in portfolio management, structural engineering and fluid dynamics illustrate the broad applicability of our method and its advantages over classical sampling-based approaches in terms of both accuracy and efficiency.

We study the temperature control problem for Langevin diffusions in the context of non-convex optimization. The classical optimal control of such a problem is of the bang-bang type, which is overly sensitive to errors. A remedy is to allow the diffusions to explore other temperature values and hence smooth out the bang-bang control. We accomplish this by a stochastic relaxed control formulation incorporating randomization of the temperature control and regularizing its entropy. We derive a state-dependent, truncated exponential distribution, which can be used to sample temperatures in a Langevin algorithm, in terms of the solution to an HJB partial differential equation. We carry out a numerical experiment on a one-dimensional baseline example, in which the HJB equation can be easily solved, to compare the performance of the algorithm with three other available algorithms in search of a global optimum.

Random projection techniques based on Johnson-Lindenstrauss lemma are used for randomly aggregating the constraints or variables of optimization problems while approximately preserving their optimal values, that leads to smaller-scale optimization problems. D'Ambrosio et al. have applied random projection to a quadratic optimization problem so as to decrease the number of decision variables. Although the problem size becomes smaller, the projected problem will also almost surely be non-convex if the original problem is non-convex, and hence will be hard to solve. In this paper, by focusing on the fact that the level of the non-convexity of a non-convex quadratic optimization problem can be alleviated by random projection, we find an approximate global optimal value of the problem by attributing it to a convex problem with smaller size. To the best of our knowledge, our paper is the first to use random projection for convexification of non-convex optimization problems. We evaluate the approximation error between optimum values of a non-convex optimization problem and its convexified randomly projected problem., 25 pages, 4 figures