Search

International audience; We introduce a model of graph-constrained dynamic choice with reinforcement modeled by positively $\alpha$-homogeneous rewards. We show that its empirical process, which can be written as a stochastic approximation recursion with Markov noise, has the same probability law as a certain vertex reinforced random walk. We use this equivalence to show that for $\alpha > 0$, the asymptotic outcome concentrates around the optimum in a certain limiting sense when 'annealed' by letting $\alpha \to \infty$ slowly.

We define and analyze the notion of variational Wardrop equilibrium for nonatomic aggregative games with an infinity of players types. These equilibria are characterized through an infinite-dimensional variational inequality. We show, under monotonicity conditions, a convergence theorem enables to approximate such an equilibrium with arbitrary precision. To this end, we introduce a sequence of nonatomic games with a finite number of players types, which approximates the initial game. We show the existence of a symmetric Wardrop equilibrium in each of these games. We prove that those symmetric equilibria converge to an equilibrium of the infinite game, and that they can be computed as solutions of finite-dimensional variational inequalities. The model is illustrated through an example from smart grids: the description of a large population of electricity consumers by a parametric distribution gives a nonatomic game with an infinity of different players types, with actions subject to coupling constraints., arXiv admin note: substantial text overlap with arXiv:1806.06230

We propose a computationally efficient method to construct nonparametric, heteroscedastic prediction bands for uncertainty quantification, with or without any user-specified predictive model. Our approach provides an alternative to the now-standard conformal prediction for uncertainty quantification, with novel theoretical insights and computational advantages. The data-adaptive prediction band is universally applicable with minimal distributional assumptions, has strong non-asymptotic coverage properties, and is easy to implement using standard convex programs. Our approach can be viewed as a novel variance interpolation with confidence and further leverages techniques from semi-definite programming and sum-of-squares optimization. Theoretical and numerical performances for the proposed approach for uncertainty quantification are analyzed., 21 pages, 4 figures

Following the breakthrough work of Tardos (Oper Res 34:250–256, 1986) in the bit-complexity model, Vavasis and Ye (Math Program 74(1):79–120, 1996) gave the first exact algorithm for linear programming in the real model of computation with running time depending only on the constraint matrix. For solving a linear program (LP) $$\max \, c^\top x,\, Ax = b,\, x \ge 0,\, A \in \mathbb {R}^{m \times n}$$ max c ⊤ x , A x = b , x ≥ 0 , A ∈ R m × n , Vavasis and Ye developed a primal-dual interior point method using a ‘layered least squares’ (LLS) step, and showed that $$O(n^{3.5} \log (\bar{\chi }_A+n))$$ O ( n 3.5 log ( χ ¯ A + n ) ) iterations suffice to solve (LP) exactly, where $$\bar{\chi }_A$$ χ ¯ A is a condition measure controlling the size of solutions to linear systems related to A. Monteiro and Tsuchiya (SIAM J Optim 13(4):1054–1079, 2003), noting that the central path is invariant under rescalings of the columns of A and c, asked whether there exists an LP algorithm depending instead on the measure $$\bar{\chi }^*_A$$ χ ¯ A ∗ , defined as the minimum $$\bar{\chi }_{AD}$$ χ ¯ AD value achievable by a column rescaling AD of A, and gave strong evidence that this should be the case. We resolve this open question affirmatively. Our first main contribution is an $$O(m^2 n^2 + n^3)$$ O ( m 2 n 2 + n 3 ) time algorithm which works on the linear matroid of A to compute a nearly optimal diagonal rescaling D satisfying $$\bar{\chi }_{AD} \le n(\bar{\chi }_A^*)^3$$ χ ¯ AD ≤ n ( χ ¯ A ∗ ) 3 . This algorithm also allows us to approximate the value of $$\bar{\chi }_A$$ χ ¯ A up to a factor $$n (\bar{\chi }_A^*)^2$$ n ( χ ¯ A ∗ ) 2 . This result is in surprising contrast to that of Tunçel (Math Program 86(1):219–223, 1999), who showed NP-hardness for approximating $$\bar{\chi }_A$$ χ ¯ A to within $$2^{\textrm{poly}(\textrm{rank}(A))}$$ 2 poly ( rank ( A ) ) . The key insight for our algorithm is to work with ratios $$g_i/g_j$$ g i / g j of circuits of A—i.e., minimal linear dependencies $$Ag=0$$ A g = 0 —which allow us to approximate the value of $$\bar{\chi }_A^*$$ χ ¯ A ∗ by a maximum geometric mean cycle computation in what we call the ‘circuit ratio digraph’ of A. While this resolves Monteiro and Tsuchiya’s question by appropriate preprocessing, it falls short of providing either a truly scaling invariant algorithm or an improvement upon the base LLS analysis. In this vein, as our second main contribution we develop a scaling invariant LLS algorithm, which uses and dynamically maintains improving estimates of the circuit ratio digraph, together with a refined potential function based analysis for LLS algorithms in general. With this analysis, we derive an improved $$O(n^{2.5} \log (n)\log (\bar{\chi }^*_A+n))$$ O ( n 2.5 log ( n ) log ( χ ¯ A ∗ + n ) ) iteration bound for optimally solving (LP) using our algorithm. The same argument also yields a factor $$n/\log n$$ n / log n improvement on the iteration complexity bound of the original Vavasis–Ye algorithm.

We propose subsampling as a unified algorithmic technique for submodular maximization in centralized and online settings. The idea is simple: independently sample elements from the ground set, and use simple combinatorial techniques (such as greedy or local search) on these sampled elements. We show that this approach leads to optimal/state-of-the-art results despite being much simpler than existing methods. In the usual offline setting, we present SampleGreedy, which obtains a $(p + 2 + o(1))$-approximation for maximizing a submodular function subject to a $p$-extendible system using $O(n + nk/p)$ evaluation and feasibility queries, where $k$ is the size of the largest feasible set. The approximation ratio improves to $p+1$ and $p$ for monotone submodular and linear objectives, respectively. In the streaming setting, we present SampleStreaming, which obtains a $(4p +2 - o(1))$-approximation for maximizing a submodular function subject to a $p$-matchoid using $O(k)$ memory and $O(km/p)$ evaluation and feasibility queries per element, where $m$ is the number of matroids defining the $p$-matchoid. The approximation ratio improves to $4p$ for monotone submodular objectives. We empirically demonstrate the effectiveness of our algorithms on video summarization, location summarization, and movie recommendation tasks., Comment: arXiv admin note: text overlap with arXiv:1802.07098

We study a noisy tensor completion problem of broad practical interest, namely, the reconstruction of a low-rank tensor from highly incomplete and randomly corrupted observations of its entries. While a variety of prior work has been dedicated to this problem, prior algorithms either are computationally too expensive for large-scale applications, or come with sub-optimal statistical guarantees. Focusing on "incoherent" and well-conditioned tensors of a constant CP rank, we propose a two-stage nonconvex algorithm -- (vanilla) gradient descent following a rough initialization -- that achieves the best of both worlds. Specifically, the proposed nonconvex algorithm faithfully completes the tensor and retrieves all individual tensor factors within nearly linear time, while at the same time enjoying near-optimal statistical guarantees (i.e. minimal sample complexity and optimal estimation accuracy). The estimation errors are evenly spread out across all entries, thus achieving optimal $\ell_{\infty}$ statistical accuracy. We have also discussed how to extend our approach to accommodate asymmetric tensors. The insight conveyed through our analysis of nonconvex optimization might have implications for other tensor estimation problems., Comment: Accepted to Operations Research

In this paper, we present a new model and two mechanisms for auctions in two-sided markets of buyers and sellers, where budget constraints are imposed on buyers. Our model incorporates polymatroidal environments, and is applicable to a wide variety of models that include multiunit auctions, matching markets and reservation exchange markets. Our mechanisms are build on polymatroidal network flow model by Lawler and Martel, and enjoy various nice properties such as incentive compatibility of buyers, individual rationality, pareto optimality, strong budget balance. The first mechanism is a simple "reduce-to-recover" algorithm that reduces the market to be one-sided, applies the polyhedral clinching auction by Goel et al, and lifts the resulting allocation to the original two-sided market via polymatroidal network flow. The second mechanism is a two-sided generalization of the polyhedral clinching auction, which improves the first mechanism in terms of the fairness of revenue sharing on sellers. Both mechanisms are implemented by polymatroid algorithms. We demonstrate how our framework is applied to internet display ad auctions.

Differential-algebraic equations (DAEs) are widely used for modeling of dynamical systems. In numerical analysis of DAEs, consistent initialization and index reduction are important preprocessing prior to numerical integration. Existing DAE solvers commonly adopt structural preprocessing methods based on combinatorial optimization. Unfortunately, the structural methods fail if the DAE has numerical or symbolic cancellations. For such DAEs, methods have been proposed to modify them to other DAEs to which the structural methods are applicable, based on the combinatorial relaxation technique. Existing modification methods, however, work only for a class of DAEs that are linear or close to linear. This paper presents two new modification methods for nonlinear DAEs: the substitution method and the augmentation method. Both methods are based on the combinatorial relaxation approach and are applicable to a large class of nonlinear DAEs. The substitution method symbolically solves equations for some derivatives based on the implicit function theorem and substitutes the solution back into the system. Instead of solving equations, the augmentation method modifies DAEs by appending new variables and equations. The augmentation method has advantages that the equation solving is not needed and the sparsity of DAEs is retained. It is shown in numerical experiments that both methods, especially the augmentation method, successfully modify high-index DAEs that the DAE solver in MATLAB cannot handle., Comment: A preliminary version of this paper is to appear in Proceedings of the 44th International Symposium on Symbolic and Algebraic Computation (ISSAC 2019), Beijing, China, July 2019

Matrix completion is about recovering a matrix from its partial revealed entries, and it can often be achieved by exploiting the inherent simplicity or low dimensional structure of the target matrix. For instance, a typical notion of matrix simplicity is low rank. In this paper we study matrix completion based on another low dimensional structure, namely the low rank Hankel structure in the Fourier domain. It is shown that matrices with this structure can be exactly recovered by solving a convex optimization program provided the sampling complexity is nearly optimal. Empirical results are also presented to justify the effectiveness of the convex method.

Given a graph $G=(V,E)$ with a weight $w_v$ associated with each vertex $v\in V$, the maximum weighted induced forest problem (MWIF) consists of encountering a maximum weighted subset $V'\subseteq V$ of the vertices such that $V'$ induces a forest. This NP-hard problem is known to be equivalent to the minimum weighted feedback vertex set problem, which has applicability in a variety of domains. The closely related maximum weighted induced tree problem (MWIT), on the other hand, requires that the subset $V'\subseteq V$ induces a tree. We propose two new integer programming formulations with an exponential number of constraints and branch-and-cut procedures for MWIF. Computational experiments using benchmark instances are performed comparing several formulations, including the newly proposed approaches and those available in the literature, when solved by a standard commercial mixed integer programming solver. More specifically, five formulations are compared, two compact (i.e., with a polynomial number of variables and constraints) ones and the three others with an exponential number of constraints. The experiments show that a new formulation for the problem based on directed cutset inequalities for eliminating cycles (DCUT) offers stronger linear relaxation bounds earlier in the search process. The results also indicate that the other new formulation, denoted tree with cycle elimination (TCYC), outperforms those available in the literature when it comes to the average times for proving optimality for the small instances, especially the more challenging ones. Additionally, this formulation can achieve much lower average times for solving the larger random instances that can be optimally solved. Furthermore, we show how the formulations for MWIF can be easily extended for MWIT. Such extension allowed us to compare the optimal solution values of the two problems.

We provide a control-theoretic perspective on optimal tensor algorithms for minimizing a convex function in a finite-dimensional Euclidean space. Given a function $\Phi: \mathbb{R}^d \rightarrow \mathbb{R}$ that is convex and twice continuously differentiable, we study a closed-loop control system that is governed by the operators $\nabla \Phi$ and $\nabla^2 \Phi$ together with a feedback control law $\lambda(\cdot)$ satisfying the algebraic equation $(\lambda(t))^p\|\nabla\Phi(x(t))\|^{p-1} = \theta$ for some $\theta \in (0, 1)$. Our first contribution is to prove the existence and uniqueness of a local solution to this system via the Banach fixed-point theorem. We present a simple yet nontrivial Lyapunov function that allows us to establish the existence and uniqueness of a global solution under certain regularity conditions and analyze the convergence properties of trajectories. The rate of convergence is $O(1/t^{(3p+1)/2})$ in terms of objective function gap and $O(1/t^{3p})$ in terms of squared gradient norm. Our second contribution is to provide two algorithmic frameworks obtained from discretization of our continuous-time system, one of which generalizes the large-step A-HPE framework and the other of which leads to a new optimal $p$-th order tensor algorithm. While our discrete-time analysis can be seen as a simplification and generalization of~\citet{Monteiro-2013-Accelerated}, it is largely motivated by the aforementioned continuous-time analysis, demonstrating the fundamental role that the feedback control plays in optimal acceleration and the clear advantage that the continuous-time perspective brings to algorithmic design. A highlight of our analysis is that we show that all of the $p$-th order optimal tensor algorithms that we discuss minimize the squared gradient norm at a rate of $O(k^{-3p})$, which complements the recent analysis., Comment: Accepted by Mathematical Programming Series A; 45 pages

We present two sampled quasi-Newton methods (sampled LBFGS and sampled LSR1) for solving empirical risk minimization problems that arise in machine learning. Contrary to the classical variants of these methods that sequentially build Hessian or inverse Hessian approximations as the optimization progresses, our proposed methods sample points randomly around the current iterate at every iteration to produce these approximations. As a result, the approximations constructed make use of more reliable (recent and local) information, and do not depend on past iterate information that could be significantly stale. Our proposed algorithms are efficient in terms of accessed data points (epochs) and have enough concurrency to take advantage of parallel/distributed computing environments. We provide convergence guarantees for our proposed methods. Numerical tests on a toy classification problem as well as on popular benchmarking binary classification and neural network training tasks reveal that the methods outperform their classical variants., Comment: 50 pages, 33 figures

Advances in multi-spectral detectors are causing a paradigm shift in X-ray Computed Tomography (CT). Spectral information acquired from these detectors can be used to extract volumetric material composition maps of the object of interest. If the materials and their spectral responses are known a priori, the image reconstruction step is rather straightforward. If they are not known, however, the maps as well as the responses need to be estimated jointly. A conventional workflow in spectral CT involves performing volume reconstruction followed by material decomposition, or vice versa. However, these methods inherently suffer from the ill-posedness of the joint reconstruction problem. To resolve this issue, we propose 'A Dictionary-based Joint reconstruction and Unmixing method for Spectral Tomography' (ADJUST). Our formulation relies on forming a dictionary of spectral signatures of materials common in CT and prior knowledge of the number of materials present in an object. In particular, we decompose the spectral volume linearly in terms of spatial material maps, a spectral dictionary, and the indicator of materials for the dictionary elements. We propose a memory-efficient accelerated alternating proximal gradient method to find an approximate solution to the resulting bi-convex problem. From numerical demonstrations on several synthetic phantoms, we observe that ADJUST performs exceedingly well compared to other state-of-the-art methods. Additionally, we address the robustness of ADJUST against limited and noisy measurement patterns. The demonstration of the proposed approach on a spectral micro-CT dataset shows its potential for real-world applications. Code is available at https://github.com/mzeegers/ADJUST., Comment: This paper is under consideration at Inverse Problems with minor revisions. 33 pages, 24 figures

We study stochastic projection-free methods for constrained optimization of smooth functions on Riemannian manifolds, i.e., with additional constraints beyond the parameter domain being a manifold. Specifically, we introduce stochastic Riemannian Frank-Wolfe methods for nonconvex and geodesically convex problems. We present algorithms for both purely stochastic optimization and finite-sum problems. For the latter, we develop variance-reduced methods, including a Riemannian adaptation of the recently proposed Spider technique. For all settings, we recover convergence rates that are comparable to the best-known rates for their Euclidean counterparts. Finally, we discuss applications to two classic tasks: The computation of the Karcher mean of positive definite matrices and Wasserstein barycenters for multivariate normal distributions. For both tasks, stochastic Fw methods yield state-of-the-art empirical performance., Under Review

Let A ∈ Zm×n be an integral matrix and a, b, c ∈ Z satisfy a ≥ b ≥ c ≥ 0. The question is to recognize whether A is {a, b, c}-modular, i.e., whether the set of n×n subdeterminants of A in absolute value is {a, b, c}. We will succeed in solving this problem in polynomial time unless A possesses a duplicative relation, that is, A has nonzero n × n subdeterminants k1 and k2 satisfying 2·|k1|=|k2|. This is an extension of the well-known recognition algorithm for totally unimodular matrices. As a consequence of our analysis, we present a polynomial time algorithm to solve integer programs in standard form over {a, b, c}-modular constraint matrices for any constants a, b and c., Vietnam Journal of Mathematics, 50 (2), ISSN:2305-2228, ISSN:0866-7179, ISSN:2305-221X

This paper is concerned with the squared F(robenius)-norm regularized factorization form for noisy low-rank matrix recovery problems. Under a suitable assumption on the restricted condition number of the Hessian matrix of the loss function, we establish an error bound to the true matrix for the non-strict critical points with rank not more than that of the true matrix. Then, for the squared F-norm regularized factorized least squares loss function, we establish its KL property of exponent 1/2 on the global optimal solution set under the noisy and full sample setting, and achieve this property at its certain class of critical points under the noisy and partial sample setting. These theoretical findings are also confirmed by solving the squared F-norm regularized factorization problem with an accelerated alternating minimization method.

We introduce a novel method for the implementation of shape optimization for non-parameterized shapes in fluid dynamics applications, where we propose to use the shape derivative to determine deformation fields with the help of the $$p-$$ p - Laplacian for $$p > 2$$ p > 2 . This approach is closely related to the computation of steepest descent directions of the shape functional in the $$W^{1,\infty }-$$ W 1 , ∞ - topology and refers to the recent publication Deckelnick et al. (A novel $$W^{1,\infty}$$ W 1 , ∞ approach to shape optimisation with Lipschitz domains, 2021), where this idea is proposed. Our approach is demonstrated for shape optimization related to drag-minimal free floating bodies. The method is validated against existing approaches with respect to convergence of the optimization algorithm, the obtained shape, and regarding the quality of the computational grid after large deformations. Our numerical results strongly indicate that shape optimization related to the $$W^{1,\infty }$$ W 1 , ∞ -topology—though numerically more demanding—seems to be superior over the classical approaches invoking Hilbert space methods, concerning the convergence, the obtained shapes and the mesh quality after large deformations, in particular when the optimal shape features sharp corners.

Many geometric optimization problems can be reduced to finding points in space (centers) minimizing an objective function which continuously depends on the distances from the centers to given input points. Examples are k-Means, Geometric k-Median/Center, Continuous Facility Location, m-Variance, etc. We prove that, for any fixed $$\varepsilon >0$$ , every set of n input points in fixed-dimensional space with the metric induced by any vector norm admits a set of O(n) candidate centers which can be computed in almost linear time and which contains a $$(1+\varepsilon )$$ -approximation of each point of space with respect to the distances to all the input points. It gives a universal approximation-preserving reduction of geometric center-based problems with arbitrary continuity-type objective functions to their discrete versions where the centers are selected from a fairly small set of candidates. The existence of such a linear-size set of candidates is also shown for any metric space of fixed doubling dimension.

We present a new multilevel minimization framework for the training of deep residual networks (ResNets), which has the potential to significantly reduce training time and effort. Our framework is based on the dynamical system’s viewpoint, which formulates a ResNet as the discretization of an initial value problem. The training process is then formulated as a time-dependent optimal control problem, which we discretize using different time-discretization parameters, eventually generating multilevel-hierarchy of auxiliary networks with different resolutions. The training of the original ResNet is then enhanced by training the auxiliary networks with reduced resolutions. By design, our framework is conveniently independent of the choice of the training strategy chosen on each level of the multilevel hierarchy. By means of numerical examples, we analyze the convergence behavior of the proposed method and demonstrate its robustness. For our examples we employ a multilevel gradient-based methods. Comparisons with standard single level methods show a speedup of more than factor three while achieving the same validation accuracy.

Gradient-based optimization algorithms can be studied from the perspective of limiting ordinary differential equations (ODEs). Motivated by the fact that existing ODEs do not distinguish between two fundamentally different algorithms---Nesterov's accelerated gradient method for strongly convex functions (NAG-SC) and Polyak's heavy-ball method---we study an alternative limiting process that yields high-resolution ODEs. We show that these ODEs permit a general Lyapunov function framework for the analysis of convergence in both continuous and discrete time. We also show that these ODEs are more accurate surrogates for the underlying algorithms; in particular, they not only distinguish between NAG-SC and Polyak's heavy-ball method, but they allow the identification of a term that we refer to as "gradient correction" that is present in NAG-SC but not in the heavy-ball method and is responsible for the qualitative difference in convergence of the two methods. We also use the high-resolution ODE framework to study Nesterov's accelerated gradient method for (non-strongly) convex functions, uncovering a hitherto unknown result---that NAG-C minimizes the squared gradient norm at an inverse cubic rate. Finally, by modifying the high-resolution ODE of NAG-C, we obtain a family of new optimization methods that are shown to maintain the accelerated convergence rates of NAG-C for smooth convex functions., 82 pages, 11 figures

The line sum optimization problem asks for a ( 0 , 1 ) -matrix minimizing the sum of given functions evaluated at its row and column sums. We show that the uniform problem, with identical row functions and identical column functions, and the monotone problem, over matrices with nonincreasing row and column sums, are polynomial time solvable.

We show algorithms for computing representative families for matroid intersections and use them in fixed-parameter algorithms for set packing, set covering, and facility location problems with multiple matroid constraints. We complement our tractability results by hardness results., Comment: Restructuring (focus on representative families instead of facility location), slight running time improvements, algorithms factored out

We address the problem of prescribing an optimal decision in a framework where the cost function depends on uncertain problem parameters that need to be learned from data. Earlier work proposed prescriptive formulations based on supervised machine learning methods. These prescriptive methods can factor in contextual information on a potentially large number of covariates to take context specific actions which are superior to any static decision. When working with noisy or corrupt data, however, such nominal prescriptive methods can be prone to adverse overfitting phenomena and fail to generalize on out-of-sample data. In this paper we combine ideas from robust optimization and the statistical bootstrap to propose novel prescriptive methods which safeguard against overfitting. We show indeed that a particular entropic robust counterpart to such nominal formulations guarantees good performance on synthetic bootstrap data. As bootstrap data is often a sensible proxy to actual out-of-sample data, our robust counterpart can be interpreted to directly encourage good out-of-sample performance. The associated robust prescriptive methods furthermore reduce to convenient tractable convex optimization problems in the context of local learning methods such as nearest neighbors and Nadaraya-Watson learning. We illustrate our data-driven decision-making framework and our novel robustness notion on a small newsvendor problem.

Wasserstein balls, which contain all probability measures within a pre-specified Wasserstein distance to a reference measure, have recently enjoyed wide popularity in the distributionally robust optimization and machine learning communities to formulate and solve data-driven optimization problems with rigorous statistical guarantees. In this technical note we prove that the Wasserstein ball is weakly compact under mild conditions, and we offer necessary and sufficient conditions for the existence of optimal solutions. We also characterize the sparsity of solutions if the Wasserstein ball is centred at a discrete reference measure. In comparison with the existing literature, which has proved similar results under different conditions, our proofs are self-contained and shorter, yet mathematically rigorous, and our necessary and sufficient conditions for the existence of optimal solutions are easily verifiable in practice.

In this paper, we introduce an asynchronous decentralized accelerated stochastic gradient descent type of algorithm for decentralized stochastic optimization. Considering communication and synchronization costs are the major bottlenecks for decentralized optimization, we attempt to reduce these costs from an algorithmic design aspect, in particular, we are able to reduce the number of agents involved in one round of update via randomization. Our major contribution is to develop a class of accelerated randomized decentralized algorithms for solving general convex composite problems. We establish $\mathcal {O}(1/\epsilon)$ (resp., $\mathcal {O}(1/\sqrt {\epsilon })$ ) communication complexity and $\mathcal {O}(1/\epsilon ^{2})$ (resp., $\mathcal {O}(1/\epsilon)$ ) sampling complexity for solving general convex (resp., strongly convex) problems. It worths mentioning that our proposing algorithm only sublinear depends on the Lipschitz constant if there is a smooth component presented in the objective function. Moreover, we also conduct some preliminary numerical experiments to demonstrate the advantages of our proposing algorithms over the state-of-the-art synchronous decentralized algorithm.

In a recent paper by Iglesias, Rumpf and Scherzer (Found. Comput. Math. 18(4), 2018) a variational model for deformations matching a pair of shapes given as level set functions was proposed. Its main feature is the presence of anisotropic energies active only in a narrow band around the hypersurfaces that resemble the behavior of elastic shells. In this work we consider some extensions and further analysis of that model. First, we present a symmetric energy functional such that given two particular shapes, it assigns the same energy to any given deformation as to its inverse when the roles of the shapes are interchanged, and introduce the adequate parameter scaling to recover a surface problem when the width of the narrow band vanishes. Then, we obtain existence of minimizing deformations for the symmetric energy in classes of bi-Sobolev homeomorphisms for small enough widths, and prove a $\Gamma$-convergence result for the corresponding non-symmetric energies as the width tends to zero. Finally, numerical results on realistic shape matching applications demonstrating the effect of the symmetric energy are presented., Comment: 28 pages, 4 figures, 2 tables. V3: Final accepted version

This paper discusses an outer-approximation guided optimization method for constrained neural network inverse problems with rectified linear units. The constrained neural network inverse problems refer to an optimization problem to find the best set of input values of a given trained neural network in order to produce a predefined desired output in presence of constraints on input values. This paper analyzes the characteristics of optimal solutions of neural network inverse problems with rectified activation units and proposes an outer-approximation algorithm by exploiting their characteristics. The proposed outer-approximation guided optimization comprises primal and dual phases. The primal phase incorporates neighbor curvatures with neighbor outer-approximations to expedite the process. The dual phase identifies and utilizes the structure of local convex regions to improve the convergence to a local optimal solution. At last, computation experiments demonstrate the superiority of the proposed algorithm compared to a projected gradient method.

A variety of optimization problems takes the form of a minimum norm optimization. In this paper, we study the change of optimal values between two incrementally constructed least norm optimization problems, with new measurements included in the second one. We prove an exact equation to calculate the change of optimal values in the linear least norm optimization problem. With the result in this paper, the change of the optimal values can be pre-calculated as a metric to guide online decision makings, without solving the second optimization problem as long the solution and covariance of the first optimization problem are available. The result can be extended to linear least distance optimization problems, and nonlinear least distance optimization with (nonlinear) equality constraints through linearizations. This derivation in this paper provides a theoretically sound explanation to the empirical observations shown in RA-L 2018 bai et al. As an additional contribution, we propose another optimization problem, i.e. aligning two trajectories at given poses, to further demonstrate how to use the metric. The accuracy of the metric is validated with numerical examples, which is quite satisfactory in general (see the experiments in RA-L 2018 bai et al.} as well), unless in some extremely adverse scenarios. Last but not least, calculating the optimal value by the proposed metric is at least one magnitude faster than solving the corresponding optimization problems directly., 6 pages on IEEE International Conference on Robotics and Automation (ICRA), 2020

Nuclear Magnetic Resonance (NMR) Spectroscopy is the second most used technique (after X-ray crystallography) for structural determination of proteins. A computational challenge in this technique involves solving a discrete optimization problem that assigns the resonance frequency to each atom in the protein. This paper introduces LIAN (LInear programming Assignment for NMR), a novel linear programming formulation of the problem which yields state-of-the-art results in simulated and experimental datasets., 28 pages, 10 figures

I study the minimax-optimal design for a two-arm controlled experiment where conditional mean outcomes may vary in a given set. When this set is permutation symmetric, the optimal design is complete randomization, and using a single partition (i.e., the design that only randomizes the treatment labels for each side of the partition) has minimax risk larger by a factor of $n-1$. More generally, the optimal design is shown to be the mixed-strategy optimal design (MSOD) of Kallus (2018). Notably, even when the set of conditional mean outcomes has structure (i.e., is not permutation symmetric), being minimax-optimal for variance still requires randomization beyond a single partition. Nonetheless, since this targets precision, it may still not ensure sufficient uniformity in randomization to enable randomization (i.e., design-based) inference by Fisher's exact test to appropriately detect violations of null. I therefore propose the inference-constrained MSOD, which is minimax-optimal among all designs subject to such uniformity constraints. On the way, I discuss Johansson et al. (2020) who recently compared rerandomization of Morgan and Rubin (2012) and the pure-strategy optimal design (PSOD) of Kallus (2018). I point out some errors therein and set straight that randomization is minimax-optimal and that the "no free lunch" theorem and example in Kallus (2018) are correct.

We propose a novel differentiable reformulation of the linearly-constrained $$\ell _1$$ ℓ 1 minimization problem, also known as the basis pursuit problem. The reformulation is inspired by the Laplacian paradigm of network theory and leads to a new family of gradient-based methods for the solution of $$\ell _1$$ ℓ 1 minimization problems. We analyze the iteration complexity of a natural solution approach to the reformulation, based on a multiplicative weights update scheme, as well as the iteration complexity of an accelerated gradient scheme. The results can be seen as bounds on the complexity of iteratively reweighted least squares (IRLS) type methods of basis pursuit.

The goal of this paper is to design a simplex algorithm for linear programs on lattice polytopes that traces “short” simplex paths from any given vertex to an optimal one. We consider a lattice polytope P contained in $$[0,k]^n$$ and defined via m linear inequalities. Our first contribution is a simplex algorithm that reaches an optimal vertex by tracing a path along the edges of P of length in $$O(n^{4} k\, \hbox {log}\, k).$$ The length of this path is independent from m and it is the best possible up to a polynomial function. In fact, it is only polynomially far from the worst-case diameter, which roughly grows as nk. Motivated by the fact that most known lattice polytopes are defined via $$0,\pm 1$$ constraint matrices, our second contribution is a more sophisticated simplex algorithm which exploits the largest absolute value $$\alpha $$ of the entries in the constraint matrix. We show that the length of the simplex path generated by this algorithm is in $$O(n^2k\, \hbox {log}\, ({nk} \alpha ))$$ . In particular, if $$\alpha $$ is bounded by a polynomial in n, k, then the length of the simplex path is in $$O(n^2k\, \hbox {log}\, (nk))$$ . For both algorithms, if P is “well described”, then the number of arithmetic operations needed to compute the next vertex in the path is polynomial in n, m, and $$\hbox {log}\, k$$ . If k is polynomially bounded in n and m, the algorithm runs in strongly polynomial time.

We propose a new simple and natural algorithm for learning the optimal [Formula: see text]-value function of a discounted-cost Markov decision process (MDP) when the transition kernels are unknown. Unlike the classical learning algorithms for MDPs, such as [Formula: see text]-learning and actor-critic algorithms, this algorithm does not depend on a stochastic approximation-based method. We show that our algorithm, which we call the empirical [Formula: see text]-value iteration algorithm, converges to the optimal [Formula: see text]-value function. We also give a rate of convergence or a nonasymptotic sample complexity bound and show that an asynchronous (or online) version of the algorithm will also work. Preliminary experimental results suggest a faster rate of convergence to a ballpark estimate for our algorithm compared with stochastic approximation-based algorithms.

We present an integer programming model to compute the strong rainbow connection number, $src(G)$, of any simple graph $G$. We introduce several enhancements to the proposed model, including a fast heuristic, and a variable elimination scheme. Moreover, we present a novel lower bound for $src(G)$ which may be of independent research interest. We solve the integer program both directly and using an alternative method based on iterative lower bound improvement, the latter of which we show to be highly effective in practice. To our knowledge, these are the first computational methods for the strong rainbow connection problem. We demonstrate the efficacy of our methods by computing the strong rainbow connection numbers of graphs containing up to $379$ vertices., Comment: 27 pages, 7 figures

The network reconfiguration problem seeks to find a rooted tree $T$ such that the energy of the (unique) feasible electrical flow over $T$ is minimized. The tree requirement on the support of the flow is motivated by operational constraints in electricity distribution networks. The bulk of existing results on convex optimization over vertices of polytopes and on the structure of electrical flows do not easily give guarantees for this problem, while many heuristic methods have been developed in the power systems community as early as 1989. Our main contribution is to give the first provable approximation guarantees for the network reconfiguration problem. We provide novel lower bounds and corresponding approximation factors for various settings ranging from $\min\{O(m-n), O(n)\}$ for general graphs, to $O(\sqrt{n})$ over grids with uniform resistances on edges, and $O(1)$ for grids with uniform edge resistances and demands. To obtain the result for general graphs, we propose a new method for (approximate) spectral graph sparsification, which may be of independent interest. Using insights from our theoretical results, we propose a general heuristic for the network reconfiguration problem that is orders of magnitude faster than existing methods in the literature, while obtaining comparable performance., 37 pages, 11 figures

We revisit the problem of \textit{online linear optimization} in case the set of feasible actions is accessible through an approximated linear optimization oracle with a factor $\alpha$ multiplicative approximation guarantee. This setting is in particular interesting since it captures natural online extensions of well-studied \textit{offline} linear optimization problems which are NP-hard, yet admit efficient approximation algorithms. The goal here is to minimize the $\alpha$\textit{-regret} which is the natural extension of the standard \textit{regret} in \textit{online learning} to this setting. We present new algorithms with significantly improved oracle complexity for both the full information and bandit variants of the problem. Mainly, for both variants, we present $\alpha$-regret bounds of $O(T^{-1/3})$, were $T$ is the number of prediction rounds, using only $O(\log{T})$ calls to the approximation oracle per iteration, on average. These are the first results to obtain both average oracle complexity of $O(\log{T})$ (or even poly-logarithmic in $T$) and $\alpha$-regret bound $O(T^{-c})$ for a constant $c>0$, for both variants., Comment: Accepted to Conference on Neural Information Processing System (NIPS) 2017

Using polarity, we give an outer polyhedral approximation for the epigraph of set functions. For a submodular function, we prove that the corresponding polar relaxation is exact; hence, it is equivalent to the Lov\'asz extension. The polar approach provides an alternative proof for the convex hull description of the epigraph of a submodular function. Computational experiments show that the inequalities from outer approximations can be effective as cutting planes for solving submodular as well as non-submodular set function minimization problems.

The Nemhauser-Trotter theorem states that the standard linear programming (LP) formulation for the stable set problem has a remarkable property, also known as (weak) persistency: for every optimal LP solution that assigns integer values to some variables, there exists an optimal integer solution in which these variables retain the same values. While the standard LP is defined by only non-negativity and edge constraints, a variety of other LP formulations have been studied and one may wonder whether any of them has this property as well. We show that any other formulation that satisfies mild conditions cannot have the persistency property on all graphs, unless it is always equal to the stable set polytope., Comment: 17 pages, 6 figures

This paper investigates the variance reduction techniques Antithetic Variates (AV) and Latin Hypercube Sampling (LHS) when used for sequential sampling in stochastic programming and presents a comparative computational study. It shows conditions under which the sequential sampling with AV and LHS satisfy finite stopping guarantees and are asymptotically valid, discussing LHS in detail. It computationally compares their use in both the sequential and non-sequential settings through a collection of two-stage stochastic linear programs with different characteristics. The numerical results show that while both AV and LHS can be preferable to random sampling in either setting, LHS typically dominates in the non-sequential setting while performing well sequentially and AV gains some advantages in the sequential setting. These results imply that, given the ease of implementation of these variance reduction techniques, armed with the same theoretical properties and improved empirical performance relative to random sampling, AV and LHS sequential procedures present attractive alternatives in practice for a class of stochastic programs.

The geometric problem of estimating an unknown compact convex set from evaluations of its support function arises in a range of scientific and engineering applications. Traditional approaches typically rely on estimators that minimize the error over all possible compact convex sets; in particular, these methods do not allow for the incorporation of prior structural information about the underlying set and the resulting estimates become increasingly more complicated to describe as the number of measurements available grows. We address both of these shortcomings by describing a framework for estimating tractably specified convex sets from support function evaluations. Building on the literature in convex optimization, our approach is based on estimators that minimize the error over structured families of convex sets that are specified as linear images of concisely described sets -- such as the simplex or the spectraplex -- in a higher-dimensional space that is not much larger than the ambient space. Convex sets parametrized in this manner are significant from a computational perspective as one can optimize linear functionals over such sets efficiently; they serve a different purpose in the inferential context of the present paper, namely, that of incorporating regularization in the reconstruction while still offering considerable expressive power. We provide a geometric characterization of the asymptotic behavior of our estimators, and our analysis relies on the property that certain sets which admit semialgebraic descriptions are Vapnik-Chervonenkis (VC) classes. Our numerical experiments highlight the utility of our framework over previous approaches in settings in which the measurements available are noisy or small in number as well as those in which the underlying set to be reconstructed is non-polyhedral., Comment: 35 pages, 80 figures

We introduce iterative methods named TriCG and TriMR for solving symmetric quasi-definite systems based on the orthogonal tridiagonalization process proposed by Saunders, Simon and Yip in 1988. TriCG and TriMR are tantamount to preconditioned Block-CG and Block-MINRES with two right-hand sides in which the two approximate solutions are summed at each iteration, but require less storage and work per iteration. We evaluate the performance of TriCG and TriMR on linear systems generated from the SuiteSparse Matrix Collection and from discretized and stablized Stokes equations. We compare TriCG and TriMR with SYMMLQ and MINRES, the recommended Krylov methods for symmetric and indefinite systems. In all our experiments, TriCG and TriMR terminate earlier than SYMMLQ and MINRES on a residual-based stopping condition with an improvement of up to 50% in terms of number of iterations. They also terminate more reliably than Block-CG and Block-MINRES. Experiments in quadruple and octuple precision suggest that loss of orthogonality in the basis vectors is significantly less pronounced in TriCG and TriMR than in Block-CG and Block-MINRES., Comment: 24 pages, 12 figures

This paper studies the problem of decomposing a low-rank positive-semidefinite matrix into symmetric factors with binary entries, either $\{\pm 1\}$ or $\{0,1\}$. This research answers fundamental questions about the existence and uniqueness of these decompositions. It also leads to tractable factorization algorithms that succeed under a mild deterministic condition. A companion paper addresses the related problem of decomposing a low-rank rectangular matrix into a binary factor and an unconstrained factor., 21(+4) pages, 3 figures

We propose a unified framework to address a family of classical mixed-integer optimization problems with logically constrained decision variables, including network design, facility location, unit commitment, sparse portfolio selection, binary quadratic optimization, sparse principal analysis and sparse learning problems. These problems exhibit logical relationships between continuous and discrete variables, which are usually reformulated linearly using a big-M formulation. In this work, we challenge this longstanding modeling practice and express the logical constraints in a non-linear way. By imposing a regularization condition, we reformulate these problems as convex binary optimization problems, which are solvable using an outer-approximation procedure. In numerical experiments, we establish that a general-purpose numerical strategy, which combines cutting-plane, first-order and local search methods, solves these problems faster and at a larger scale than state-of-the-art mixed-integer linear or second-order cone methods. Our approach successfully solves network design problems with 100s of nodes and provides solutions up to 40\% better than the state-of-the-art; sparse portfolio selection problems with up to 3,200 securities compared with 400 securities for previous attempts; and sparse regression problems with up to 100,000 covariates., Revised version (including title change). The old title was "A unified approach to mixed-integer optimization: Nonlinear formulations and scalable algorithms"

Clustering is one of the fundamental tasks in data analytics and machine learning. In many situations, different clusterings of the same data set become relevant. For example, different algorithms for the same clustering task may return dramatically different solutions. We are interested in applications in which one clustering has to be transformed into another; e.g., when a gradual transition from an old solution to a new one is required. In this paper, we devise methods for constructing such a transition based on linear programming and network theory. We use a so-called clustering-difference graph to model the desired transformation and provide methods for decomposing the graph into a sequence of elementary moves that accomplishes the transformation. These moves are equivalent to the edge directions, or circuits, of the underlying partition polytopes. Therefore, in addition to a conceptually new metric for measuring the distance between clusterings, we provide new bounds on the circuit diameter of these partition polytopes.

In this manuscript, we study the properties of a family of second-order differential equations with damping, its discretizations and their connections with accelerated optimization algorithms for $m$-strongly convex and $L$-smooth functions. In particular, using the Linear Matrix Inequality LMI framework developed by \emph{Fazlyab et. al. $(2018)$}, we derive analytically a (discrete) Lyapunov function for a two-parameter family of Nesterov optimization methods, which allows for the complete characterization of their convergence rate. In the appropriate limit, this family of methods may be seen as a discretization of a family of second-order ordinary differential equations for which we construct(continuous) Lyapunov functions by means of the LMI framework. The continuous Lyapunov functions may alternatively, be obtained by studying the limiting behaviour of their discrete counterparts. Finally, we show that the majority of typical discretizations of the family of ODEs, such as the Heavy ball method, do not possess Lyapunov functions with properties similar to those of the Lyapunov function constructed here for the Nesterov method., 21 pages, 1 figure

Low rank matrix recovery problems appear widely in statistics, combinatorics, and imaging. One celebrated method for solving these problems is to formulate and solve a semidefinite program (SDP). It is often known that the exact solution to the SDP with perfect data recovers the solution to the original low rank matrix recovery problem. It is more challenging to show that an approximate solution to the SDP formulated with noisy problem data acceptably solves the original problem; arguments are usually ad hoc for each problem setting, and can be complex. In this note, we identify a set of conditions that we call simplicity that limit the error due to noisy problem data or incomplete convergence. In this sense, simple SDPs are robust: simple SDPs can be (approximately) solved efficiently at scale; and the resulting approximate solutions, even with noisy data, can be trusted. Moreover, we show that simplicity holds generically, and also for many structured low rank matrix recovery problems, including the stochastic block model, $\mathbb{Z}_2$ synchronization, and matrix completion. Formally, we call an SDP simple if it has a surjective constraint map, admits a unique primal and dual solution pair, and satisfies strong duality and strict complementarity. However, simplicity is not a panacea: we show the Burer-Monteiro formulation of the SDP may have spurious second-order critical points, even for a simple SDP with a rank 1 solution., Comment: 24 pages, 1 figure, and 1 table

This paper studies the problem of maximizing the sum of traces of matrix quadratic forms on a product of Stiefel manifolds. This orthogonal trace-sum maximization (OTSM) problem generalizes many interesting problems such as generalized canonical correlation analysis (CCA), Procrustes analysis, and cryo-electron microscopy of the Nobel prize fame. For these applications finding global solutions is highly desirable but it has been unclear how to find even a stationary point, let alone testing its global optimality. Through a close inspection of Ky Fan's classical result (1949) on the variational formulation of the sum of largest eigenvalues of a symmetric matrix, and a semidefinite programming (SDP) relaxation of the latter, we first provide a simple method to certify global optimality of a given stationary point of OTSM. This method only requires testing whether a symmetric matrix is positive semidefinite. A by-product of this analysis is an unexpected strong duality between Shapiro-Botha (1988) and Zhang-Singer (2017). After showing that a popular algorithm for generalized CCA and Procrustes analysis may generate oscillating iterates, we propose a simple fix that provably guarantees convergence to a stationary point. The combination of our algorithm and certificate reveals novel global optima of various instances of OTSM., 37 pages, 3 figures

Smooth convex minimization over the unit trace-norm ball is an important optimization problem in machine learning, signal processing, statistics and other fields, that underlies many tasks in which one wishes to recover a low-rank matrix given certain measurements. While first-order methods for convex optimization enjoy optimal convergence rates, they require in worst-case to compute a full-rank SVD on each iteration, in order to compute the projection onto the trace-norm ball. These full-rank SVD computations however prohibit the application of such methods to large problems. A simple and natural heuristic to reduce the computational cost is to approximate the projection using only a low-rank SVD. This raises the question if, and under what conditions, this simple heuristic can indeed result in provable convergence to the optimal solution. In this paper we show that any optimal solution is a center of a Euclid. ball inside-which the projected-gradient mapping admits rank that is at most the multiplicity of the largest singular value of the gradient vector. Moreover, the radius of the ball scales with the spectral gap of this gradient vector. We show how this readily implies the local convergence (i.e., from a "warm-start" initialization) of standard first-order methods, using only low-rank SVD computations. We also quantify the effect of "over-parameterization", i.e., using SVD computations with higher rank, on the radius of this ball, showing it can increase dramatically with moderately larger rank. We extend our results also to the setting of optimization with trace-norm regularization and optimization over bounded-trace positive semidefinite matrices. Our theoretical investigation is supported by concrete empirical evidence that demonstrates the \textit{correct} convergence of first-order methods with low-rank projections on real-world datasets., Accepted to SIAM Journal on Optimization (SIOPT)

Conditional gradient methods have attracted much attention in both machine learning and optimization communities recently. These simple methods can guarantee the generation of sparse solutions. In addition, without the computation of full gradients, they can handle huge-scale problems sometimes even with an exponentially increasing number of decision variables. This paper aims to significantly expand the application areas of these methods by presenting new conditional gradient methods for solving convex optimization problems with general affine and nonlinear constraints. More specifically, we first present a new constraint extrapolated condition gradient (CoexCG) method that can achieve an ${\cal O}(1/\epsilon^2)$ iteration complexity for both smooth and structured nonsmooth function constrained convex optimization. We further develop novel variants of CoexCG, namely constraint extrapolated and dual regularized conditional gradient (CoexDurCG) methods, that can achieve similar iteration complexity to CoexCG but allow adaptive selection for algorithmic parameters. We illustrate the effectiveness of these methods for solving an important class of radiation therapy treatment planning problems arising from healthcare industry. To the best of our knowledge, all the algorithmic schemes and their complexity results are new in the area of projection-free methods.

We develop algebraic tools for statistical inference from samples of rotation matrices. This rests on the theory of D-modules in algebraic analysis. Noncommutative Gr\"obner bases are used to design numerical algorithms for maximum likelihood estimation, building on the holonomic gradient method of Sei, Shibata, Takemura, Ohara, and Takayama. We study the Fisher model for sampling from rotation matrices, and we apply our algorithms for data from the applied sciences. On the theoretical side, we generalize the underlying equivariant D-modules from SO(3) to arbitrary Lie groups. For compact groups, our D-ideals encode the normalizing constant of the Fisher model., Comment: 21 pages