Search

Numerous machine learning and industrial problems can be modeled as the minimization of a sum of $N$ so-called clipped convex functions (SCC), i.e. each term of the sum stems as the pointwise minimum between a constant and a convex function. In this work, we extend this framework to capture more problems of interest. Specifically, we allow each term of the sum to be a pointwise minimum of an arbitrary number of convex functions, called components, turning the objective into a sum of pointwise minima of convex functions (SMC). Problem (SCC) is NP-hard, highlighting an appeal for scalable local heuristics. In this spirit, one can express (SMC) objectives as the difference between two convex functions to leverage the possibility to apply (DC) algorithms to compute critical points of the problem. Our approach relies on a bi-convex reformulation of the problem. From there, we derive a family of local methods, dubbed as relaxed alternating minimization (r-AM) methods, that include classical alternating minimization (AM) as a special case. We prove that every accumulation point of r-AM is critical. In addition, we show the empirical superiority of r-AM, compared to traditional AM and (DC) approaches, on piecewise-linear regression and restricted facility location problems. Under mild assumptions, (SCC) can be cast as a mixed-integer convex program (MICP) using perspective functions. This approach can be generalized to (SMC) but introduces many copies of the primal variable. In contrast, we suggest a compact big-M based (MICP) equivalent formulation of (SMC), free of these extra variables. Finally, we showcase practical examples where solving our (MICP), restricted to a neighbourhood of a given candidate (i.e. output iterate of a local method), will either certify the candidate's optimality on that neighbourhood or providing a new point, strictly better, to restart the local method.

We consider the global minimization of a particular type of minimum structured optimization problems wherein the variables must belong to some basic set, the feasible domain is described by the intersection of a large number of functional constraints and the objective stems as the pointwise minimum of a collection of functional pieces. Among others, this setting includes all non-convex piecewise linear problems. The global minimum of a minimum structured problem can be computed using a simple enumeration scheme by solving a sequence of individual problems involving each a single piece and then taking the smallest individual minimum. We propose a new algorithm, called Upper-Lower Optimization (ULO), tackling problems from the aforementioned class. Our method does not require the solution of every individual problem listed in the baseline enumeration scheme, yielding potential substantial computational gains. It alternates between (a) local-search steps, which minimize upper models, and (b) suggestion steps, which optimize lower models. We prove that ULO can return a solution of any prescribed global optimality accuracy. According to a computational complexity model based on the cost of solving the individual problems, we analyze in which practical scenarios ULO is expected to outperform the baseline enumeration scheme. Finally, we empirically validate our approach on a set of piecewise linear programs (PL) of incremental difficulty.

This paper studies the last iterate of subgradient method with Polyak step size when applied to the minimization of a nonsmooth convex function with bounded subgradients. We show that the subgradient method with Polyak step size achieves a convergence rate $\mathcal{O}\left(\tfrac{1}{\sqrt[4]{N}}\right)$ in terms of the final iterate. An example is provided to show that this rate is exact and cannot be improved. We introduce an adaptive Polyak step size for which the subgradient method enjoys a convergence rate $\mathcal{O}\left(\tfrac{1}{\sqrt{N}}\right)$ for the last iterate. Its convergence rate matches exactly the lower bound on the performance of any black-box method on the considered problem class. Additionally, we propose an adaptive Polyak method with a momentum term, where the step sizes are independent of the number of iterates. We establish that the algorithm also attains the optimal convergence rate. We investigate the alternating projection method. We derive a convergence rate $\left( \frac{2N }{ 2N+1 } \right)^N\tfrac{R}{\sqrt{2N+1}}$ for the last iterate, where $R$ is a bound on the distance between the initial iterate and a solution. An example is also provided to illustrate the exactness of the rate.

We consider gradient descent with constant stepsizes and derive exact worst-case convergence rates on the minimum gradient norm of the iterates. Our analysis covers all possible stepsizes and arbitrary upper/lower bounds on the curvature of the objective function, thus including convex, strongly convex and weakly convex (hypoconvex) objective functions. Unlike prior approaches, which rely solely on inequalities connecting consecutive iterations, our analysis employs inequalities involving an iterate and its two predecessors. While this complicates the proofs to some extent, it enables us to achieve, for the first time, an exact full-range analysis of gradient descent for any constant stepsizes (covering, in particular, normalized stepsizes greater than one), whereas the literature contained only conjectured rates of this type. In the nonconvex case, allowing arbitrary bounds on upper and lower curvatures extends existing partial results that are valid only for gradient Lipschitz functions (i.e., where lower and upper bounds on curvature are equal), leading to improved rates for weakly convex functions. From our exact rates, we deduce the optimal constant stepsize for gradient descent. Leveraging our analysis, we also introduce a new variant of gradient descent based on a unique, fixed sequence of variable stepsizes, demonstrating its superiority over any constant stepsize schedule., Comment: This paper replaces the arXiv submission of the authors: Tight convergence rates of the gradient method on smooth hypoconvex functions (arXiv:2203.00775)

Some early violins have been reduced during their history to fit imposed morphological standards, while more recent ones have been built directly to these standards. We can observe differences between reduced and unreduced instruments, particularly in their contour lines and channel of minima. In a recent preliminary work, we computed and highlighted those two features for two instruments using triangular 3D meshes acquired by photogrammetry, whose fidelity has been assessed and validated with sub-millimetre accuracy. We propose here an extension to a corpus of 38 violins, violas and cellos, and introduce improved procedures, leading to a stronger discussion of the geometric analysis. We first recall the material we are working with. We then discuss how to derive the best reference plane for the violin alignment, which is crucial for the computation of contour lines and channel of minima. Finally, we show how to compute efficiently both characteristics and we illustrate our results with a few examples., Comment: Paper accepted (before reviewing) for the Florence Heri-Tech 2024 Conference

We consider a difference-of-convex formulation where one of the terms is allowed to be hypoconvex (or weakly convex). We first examine the precise behavior of a single iteration of the Difference-of-Convex algorithm (DCA), giving a tight characterization of the objective function decrease. This requires distinguishing between eight distinct parameter regimes. Our proofs are inspired by the performance estimation framework, but are much simplified compared to similar previous work. We then derive sublinear DCA convergence rates towards critical points, distinguishing between cases where at least one of the functions is smooth and where both functions are nonsmooth. We conjecture the tightness of these rates for four parameter regimes, based on strong numerical evidence obtained via performance estimation, as well as the leading constant in the asymptotic sublinear rate for two more regimes.

Selecting the fastest algorithm for a specific signal/image processing task is a challenging question. We propose an approach based on the Performance Estimation Problem framework that numerically and automatically computes the worst-case performance of a given optimization method on a class of functions. We first propose a computer-assisted analysis and comparison of several first-order primal optimization methods, namely, the gradient method, the forward-backward, Peaceman-Rachford, and Douglas-Rachford splittings. We tighten the existing convergence results of these algorithms and extend them to new classes of functions. Our analysis is then extended and evaluated in the context of the primal-dual Chambolle-Pock and Condat-V\~u methods., Comment: 5 pages, 7 figures

Consider a sum of convex functions, where the only information known about each individual summand is the location of a minimizer. In this work, we give an exact characterization of the set of possible minimizers of the sum. Our results cover several types of assumptions on the summands, such as smoothness or strong convexity. Our main tool is the use of necessary and sufficient conditions for interpolating the considered function classes, which leads to shorter and more direct proofs in comparison with previous work. We also address the setting where each summand minimizer is assumed to lie in a unit ball, and prove a tight bound on the norm of any minimizer of the sum.

We introduce the concept of inexact first-order oracle of degree q for a possibly nonconvex and nonsmooth function, which naturally appears in the context of approximate gradient, weak level of smoothness and other situations. Our definition is less conservative than those found in the existing literature, and it can be viewed as an interpolation between fully exact and the existing inexact first-order oracle definitions. We analyze the convergence behavior of a (fast) inexact proximal gradient method using such an oracle for solving (non)convex composite minimization problems. We derive complexity estimates and study the dependence between the accuracy of the oracle and the desired accuracy of the gradient or of the objective function. Our results show that better rates can be obtained both theoretically and in numerical simulations when q is large., Comment: 21 pages

Several recent works address the impact of inexact oracles in the convergence analysis of modern first-order optimization techniques, e.g. Bregman Proximal Gradient and Prox-Linear methods as well as their accelerated variants, extending their field of applicability. In this paper, we consider situations where the oracle's inexactness can be chosen upon demand, more precision coming at a computational price counterpart. Our main motivations arise from oracles requiring the solving of auxiliary subproblems or the inexact computation of involved quantities, e.g. a mini-batch stochastic gradient as a full-gradient estimate. We propose optimal inexactness schedules according to presumed oracle cost models and patterns of worst-case guarantees, covering among others convergence results of the aforementioned methods under the presence of inexactness. Specifically, we detail how to choose the level of inexactness at each iteration to obtain the best trade-off between convergence and computational investments. Furthermore, we highlight the benefits one can expect by tuning those oracles' quality instead of keeping it constant throughout. Finally, we provide extensive numerical experiments that support the practical interest of our approach, both in offline and online settings, applied to the Fast Gradient algorithm., Comment: 36 pages, 17 figures

We study the convergence of the last iterate in subgradient methods applied to the minimization of a nonsmooth convex function with bounded subgradients. We first introduce a proof technique that generalizes the standard analysis of subgradient methods. It is based on tracking the distance between the current iterate and a different reference point at each iteration. Using this technique, we obtain the exact worst-case convergence rate for the objective accuracy of the last iterate of the projected subgradient method with either constant step sizes or constant step lengths. Tightness is shown with a worst-case instance matching the established convergence rate. We also derive the value of the optimal constant step size when performing $N$ iterations, for which we find that the last iterate accuracy is smaller than $B R \sqrt{1+\log(N)/4}/{\sqrt{N+1}}$ %$\frac{B R \log N}{\sqrt{N+1}}$ , where $B$ is a bound on the subgradient norm and $R$ is a bound on the distance between the initial iterate and a minimizer. Finally, we introduce a new optimal subgradient method that achieves the best possible last-iterate accuracy after a given number $N$ of iterations. Its convergence rate ${B R}/{\sqrt{N+1}}$ matches exactly the lower bound on the performance of any black-box method on the considered problem class. We also show that there is no universal sequence of step sizes that simultaneously achieves this optimal rate at each iteration, meaning that the dependence of the step size sequence in $N$ is unavoidable.

Kernel methods provide a powerful framework for non parametric learning. They are based on kernel functions and allow learning in a rich functional space while applying linear statistical learning tools, such as Ridge Regression or Support Vector Machines. However, standard kernel methods suffer from a quadratic time and memory complexity in the number of data points and thus have limited applications in large-scale learning. In this paper, we propose Snacks, a new large-scale solver for Kernel Support Vector Machines. Specifically, Snacks relies on a Nystr\"om approximation of the kernel matrix and an accelerated variant of the stochastic subgradient method. We demonstrate formally through a detailed empirical evaluation, that it competes with other SVM solvers on a variety of benchmark datasets., Comment: 6 pages

The Performance Estimation Problem methodology makes it possible to determine the exact worst-case performance of an optimization method. In this work, we generalize this framework to first-order methods involving linear operators. This extension requires an explicit formulation of interpolation conditions for those linear operators. We consider the class of linear operators $\mathcal{M}:x \mapsto Mx$ where matrix $M$ has bounded singular values, and the class of linear operators where $M$ is symmetric and has bounded eigenvalues. We describe interpolation conditions for these classes, i.e. necessary and sufficient conditions that, given a list of pairs $\{(x_i,y_i)\}$, characterize the existence of a linear operator mapping $x_i$ to $y_i$ for all $i$. Using these conditions, we first identify the exact worst-case behavior of the gradient method applied to the composed objective $h\circ \mathcal{M}$, and observe that it always corresponds to $\mathcal{M}$ being a scaling operator. We then investigate the Chambolle-Pock method applied to $f+g\circ \mathcal{M}$, and improve the existing analysis to obtain a proof of the exact convergence rate of the primal-dual gap. In addition, we study how this method behaves on Lipschitz convex functions, and obtain a numerical convergence rate for the primal accuracy of the last iterate. We also show numerically that averaging iterates is beneficial in this setting., Comment: Proof of the main interpolation theorem is streamlined. Added results about the Chambolle-Pock algorithm including the exact convergence rate of the primal-dual gap in a standard setting

We propose a unifying framework for the automated computer-assisted worst-case analysis of cyclic block coordinate algorithms in the unconstrained smooth convex optimization setup. We compute exact worst-case bounds for the cyclic coordinate descent and the alternating minimization algorithms over the class of smooth convex functions, and provide sublinear upper and lower bounds on the worst-case rate for the standard class of functions with coordinate-wise Lipschitz gradients. We obtain in particular a new upper bound for cyclic coordinate descent that outperforms the best available ones by an order of magnitude. We also demonstrate the flexibility of our approach by providing new numerical bounds using simpler and more natural assumptions than those normally made for the analysis of block coordinate algorithms. Finally, we provide numerical evidence for the fact that a standard scheme that provably accelerates random coordinate descent to a $O(1/k^2)$ complexity is actually inefficient when used in a (deterministic) cyclic algorithm.

Purpose: To improve target coverage and reduce the dose in the surrounding organs-at-risks (OARs), we developed an image-guided treatment method based on a precomputed library of treatment plans controlled and delivered in real-time. Methods: A library of treatment plans is constructed by optimizing a plan for each breathing phase of a 4DCT. Treatments are delivered by simulation on a continuous sequence of synthetic CTs generated from real MRI sequences. During treatment, the plans for which the tumor is at a close distance to the current tumor position are selected to deliver their spots. The study is conducted on five liver cases. Results: We tested our approach under imperfect knowledge of the tumor positions with a 2 mm distance error. On average, compared to a 4D robustly optimized treatment plan, our approach led to a dose homogeneity increase of 5% (defined as 1-(D_5-D_95)/prescription) in the target and a mean liver dose decrease of 23. The treatment time was roughly increased by a factor of 2 but remained below 4 minutes on average. Conclusions: Our image-guided treatment framework outperforms state-of-the-art 4D-robust plans for all patients in this study on both target coverage and OARs sparing, with an acceptable increase in treatment time under the current accuracy of tumor tracking technology.

Nonnegative Matrix Factorization (NMF) models are widely used to recover linearly mixed nonnegative data. When the data is made of samplings of continuous signals, the factors in NMF can be constrained to be samples of nonnegative rational functions, which allow fairly general models; this is referred to as NMF using rational functions (R-NMF). We first show that, under mild assumptions, R-NMF has an essentially unique factorization unlike NMF, which is crucial in applications where ground-truth factors need to be recovered such as blind source separation problems. Then we present different approaches to solve R-NMF: the R-HANLS, R-ANLS and R-NLS methods. From our tests, no method significantly outperforms the others, and a trade-off should be done between time and accuracy. Indeed, R-HANLS is fast and accurate for large problems, while R-ANLS is more accurate, but also more resources demanding, both in time and memory. R-NLS is very accurate but only for small problems. Moreover, we show that R-NMF outperforms NMF in various tasks including the recovery of semi-synthetic continuous signals, and a classification problem of real hyperspectral signals., Comment: 13 pages

Some early violins have been reduced during their history to fit imposed morphological standards, while more recent ones have been built directly to these standards. We propose an objective photogrammetric approach to differentiate between a reduced and an unreduced instrument, whereby a three-dimensional mesh is studied geometrically by examining 2D slices. Our contribution is twofold. First, we validate the quality of the photogrammetric mesh through a comparison with reference images obtained by medical imaging, and conclude that a sub-millimetre accuracy is achieved. Then, we show how quantitative and qualitative features such as contour lines, channel of minima and a measure of asymmetry between the upper and lower surfaces of a violin can be automatically extracted from the validated photogrammetric meshes, allowing to successfully highlight differences between instruments.

Predicting incoming failures and scheduling maintenance based on sensors information in industrial machines is increasingly important to avoid downtime and machine failure. Different machine learning formulations can be used to solve the predictive maintenance problem. However, many of the approaches studied in the literature are not directly applicable to real-life scenarios. Indeed, many of those approaches usually either rely on labelled machine malfunctions in the case of classification and fault detection, or rely on finding a monotonic health indicator on which a prediction can be made in the case of regression and remaining useful life estimation, which is not always feasible. Moreover, the decision-making part of the problem is not always studied in conjunction with the prediction phase. This paper aims to design and compare different formulations for predictive maintenance in a two-level framework and design metrics that quantify both the failure detection performance as well as the timing of the maintenance decision. The first level is responsible for building a health indicator by aggregating features using a learning algorithm. The second level consists of a decision-making system that can trigger an alarm based on this health indicator. Three degrees of refinements are compared in the first level of the framework, from simple threshold-based univariate predictive technique to supervised learning methods based on the remaining time before failure. We choose to use the Support Vector Machine (SVM) and its variations as the common algorithm used in all the formulations. We apply and compare the different strategies on a real-world rotating machine case study and observe that while a simple model can already perform well, more sophisticated refinements enhance the predictions for well-chosen parameters.

We perform the first tight convergence analysis of the gradient method with varying step sizes when applied to smooth hypoconvex (weakly convex) functions. Hypoconvex functions are smooth nonconvex functions whose curvature is bounded and assumed to belong to the interval $[\mu, L]$, with $\mu<0$. Our convergence rates improve and extend the existing analysis for smooth nonconvex functions with $L$-Lipschitz gradient (which corresponds to the case $\mu=-L$), and smoothly interpolates between that class and the class of smooth convex functions. We obtain our results using the performance estimation framework adapted to hypoconvex functions, for which new interpolation conditions are derived. We derive explicit upper bounds on the minimum gradient norm of the iterates for a large range of step sizes, explain why all such rates share a common structure, and prove that these rates are tight when step sizes are smaller or equal to $1/L$. Finally, we identify the optimal constant step size that minimizes the worst-case of the gradient method applied to hypoconvex functions.

PEPit is a Python package aiming at simplifying the access to worst-case analyses of a large family of first-order optimization methods possibly involving gradient, projection, proximal, or linear optimization oracles, along with their approximate, or Bregman variants. In short, PEPit is a package enabling computer-assisted worst-case analyses of first-order optimization methods. The key underlying idea is to cast the problem of performing a worst-case analysis, often referred to as a performance estimation problem (PEP), as a semidefinite program (SDP) which can be solved numerically. To do that, the package users are only required to write first-order methods nearly as they would have implemented them. The package then takes care of the SDP modeling parts, and the worst-case analysis is performed numerically via a standard solver., Comment: Reference work for the PEPit package (available at https://github.com/bgoujaud/PEPit)

Nonnegative matrix factorization is the following problem: given a nonnegative input matrix $V$ and a factorization rank $K$, compute two nonnegative matrices, $W$ with $K$ columns and $H$ with $K$ rows, such that $WH$ approximates $V$ as well as possible. In this paper, we propose two new approaches for computing high-quality NMF solutions using conic optimization. These approaches rely on the same two steps. First, we reformulate NMF as minimizing a concave function over a product of convex cones--one approach is based on the exponential cone, and the other on the second-order cone. Then, we solve these reformulations iteratively: at each step, we minimize exactly, over the feasible set, a majorization of the objective functions obtained via linearization at the current iterate. Hence these subproblems are convex conic programs and can be solved efficiently using dedicated algorithms. We prove that our approaches reach a stationary point with an accuracy decreasing as $\mathcal{O}(\frac{1}{i})$, where $i$ denotes the iteration number. To the best of our knowledge, our analysis is the first to provide a convergence rate to stationary points for NMF. Furthermore, in the particular cases of rank-one factorizations (that is, $K=1$), we show that one of our formulations can be expressed as a convex optimization problem implying that optimal rank-one approximations can be computed efficiently. Finally, we show on several numerical examples that our approaches are able to frequently compute exact NMFs (that is, with $V = WH$), and compete favorably with the state of the art., Comment: Accepted in Optimization Methods and Software

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

We consider a Canonical Polyadic (CP) decomposition approach to low-rank tensor completion (LRTC) by incorporating external pairwise similarity relations through graph Laplacian regularization on the CP factor matrices. The usage of graph regularization entails benefits in the learning accuracy of LRTC, but at the same time, induces coupling graph Laplacian terms that hinder the optimization of the tensor completion model. In order to solve graph-regularized LRTC, we propose efficient alternating minimization algorithms by leveraging the block structure of the underlying CP decomposition-based model. For the subproblems of alternating minimization, a linear conjugate gradient subroutine is specifically adapted to graph-regularized LRTC. Alternatively, we circumvent the complicating coupling effects of graph Laplacian terms by using an alternating directions method of multipliers. Based on the Kurdyka-{\L}ojasiewicz property, we show that the sequence generated by the proposed algorithms globally converges to a critical point of the objective function. Moreover, the complexity and convergence rate are also derived. In addition, numerical experiments including synthetic data and real data show that the graph regularized tensor completion model has improved recovery results compared to those without graph regularization, and that the proposed algorithms achieve gains in time efficiency over existing algorithms.

Mathematical psychology has a long tradition of modeling probabilistic choice via distribution-free random utility models and associated random preference models. For such models, the predicted choice probabilities often form a bounded and convex polyhedral set, or polytope. Polyhedral combinatorics have thus played a key role in studying the mathematical structure of these models. However, standard methods for characterizing the polytopes of such models are subject to a combinatorial explosion in complexity as the number of choice alternatives increases. Specifically, this is the case for random preference models based on linear, weak, semi- and interval orders. For these, a complete, linear description of the polytope is currently known only for, at most, 5--8 choice alternatives. We leverage the method of extended formulations to break through those boundaries. For each of the four types of preferences, we build an appropriate network, and show that the associated network flow polytope provides an extended formulation of the polytope of the choice model. This extended formulation has a simple linear description that is more parsimonious than descriptions obtained by standard methods for large numbers of choice alternatives. The result is a computationally less demanding way of testing the probabilistic choice model on data. We sketch how the latter interfaces with recent developments in contemporary statistics.

We provide new tools for worst-case performance analysis of the gradient (or steepest descent) method of Cauchy for smooth strongly convex functions, and Newton's method for self-concordant functions, including the case of inexact search directions. The analysis uses semidefinite programming performance estimation, as pioneered by Drori and Teboulle [Mathematical Programming, 145(1-2):451-482, 2014], and extends recent performance estimation results for the method of Cauchy by the authors [Optimization Letters, 11(7), 1185-1199, 2017]. To illustrate the applicability of the tools, we demonstrate a novel complexity analysis of short step interior point methods using inexact search directions. As an example in this framework, we sketch how to give a rigorous worst-case complexity analysis of a recent interior point method by Abernethy and Hazan [PMLR, 48:2520-2528, 2016]., Comment: 22 pages, 1 figure. Title of earlier version was "Worst-case convergence analysis of gradient and Newton methods through semidefinite programming performance estimation"

This paper considers the problem of positive semidefinite factorization (PSD factorization), a generalization of exact nonnegative matrix factorization. Given an $m$-by-$n$ nonnegative matrix $X$ and an integer $k$, the PSD factorization problem consists in finding, if possible, symmetric $k$-by-$k$ positive semidefinite matrices $\{A^1,...,A^m\}$ and $\{B^1,...,B^n\}$ such that $X_{i,j}=\text{trace}(A^iB^j)$ for $i=1,...,m$, and $j=1,...n$. PSD factorization is NP-hard. In this work, we introduce several local optimization schemes to tackle this problem: a fast projected gradient method and two algorithms based on the coordinate descent framework. The main application of PSD factorization is the computation of semidefinite extensions, that is, the representations of polyhedrons as projections of spectrahedra, for which the matrix to be factorized is the slack matrix of the polyhedron. We compare the performance of our algorithms on this class of problems. In particular, we compute the PSD extensions of size $k=1+ \lceil \log_2(n) \rceil$ for the regular $n$-gons when $n=5$, $8$ and $10$. We also show how to generalize our algorithms to compute the square root rank (which is the size of the factors in a PSD factorization where all factor matrices $A^i$ and $B^j$ have rank one) and completely PSD factorizations (which is the special case where the input matrix is symmetric and equality $A^i=B^i$ is required for all $i$)., Comment: 21 pages, 3 figures, 3 tables

We study the worst-case convergence rates of the proximal gradient method for minimizing the sum of a smooth strongly convex function and a non-smooth convex function whose proximal operator is available. We establish the exact worst-case convergence rates of the proximal gradient method in this setting for any step size and for different standard performance measures: objective function accuracy, distance to optimality and residual gradient norm. The proof methodology relies on recent developments in performance estimation of first-order methods based on semidefinite programming. In the case of the proximal gradient method, this methodology allows obtaining exact and non-asymptotic worst-case guarantees that are conceptually very simple, although apparently new. On the way, we discuss how strong convexity can be replaced by weaker assumptions, while preserving the corresponding convergence rates. We also establish that the same fixed step size policy is optimal for all three performance measures. Finally, we extend recent results on the worst-case behavior of gradient descent with exact line search to the proximal case., Comment: Journal of Optimization Theory and Algorithms (DOI:10.1007/s10957-018-1298-1) [V4: correction of a minor typo and one updated reference]

Autonomous microgrid planning is a Mixed-Integer Non Convex decision problem that requires to consider investments in both distribution and generation capacity and represents significant computation challenges. We proposed in a previous publication a deterministic Second-Order Cone (SOC) relaxation of this problem that made it computationally tractable for realsize cases. However, this problem is subject to considerable uncertainty emanating from load consumption, RES-based generation and contingencies. In this paper, we thus present a robust optimization approach that extends our previous work by including load related uncertainty at the cost of a substantial increase of the computational burden. The results show that significantly higher investment and operational costs are incurred to account for the load related uncertainty.

We consider the gradient (or steepest) descent method with exact line search applied to a strongly convex function with Lipschitz continuous gradient. We establish the exact worst-case rate of convergence of this scheme, and show that this worst-case behavior is exhibited by a certain convex quadratic function. We also give the tight worst-case complexity bound for a noisy variant of gradient descent method, where exact line-search is performed in a search direction that differs from negative gradient by at most a prescribed relative tolerance. The proofs are computer-assisted, and rely on the resolutions of semidefinite programming performance estimation problems as introduced in the paper [Y. Drori and M. Teboulle. Performance of first-order methods for smooth convex minimization: a novel approach. Mathematical Programming, 145(1-2):451-482, 2014]., Comment: 11pages

We provide a framework for computing the exact worst-case performance of any algorithm belonging to a broad class of oracle-based first-order methods for composite convex optimization, including those performing explicit, projected, proximal, conditional and inexact (sub)gradient steps. We simultaneously obtain tight worst-case guarantees and explicit instances of optimization problems on which the algorithm reaches this worst-case. We achieve this by reducing the computation of the worst-case to solving a convex semidefinite program, generalizing previous works on performance estimation by Drori and Teboulle [13] and the authors [43]. We use these developments to obtain a tighter analysis of the proximal point algorithm and of several variants of fast proximal gradient, conditional gradient, subgradient and alternating projection methods. In particular, we present a new analytical worst-case guarantee for the proximal point algorithm that is twice better than previously known, and improve the standard worst-case guarantee for the conditional gradient method by more than a factor of two. We also show how the optimized gradient method proposed by Kim and Fessler in [22] can be extended by incorporating a projection or a proximal operator, which leads to an algorithm that converges in the worst-case twice as fast as the standard accelerated proximal gradient method [2]., Comment: Published in SIOPT (updated version with corrected typo) Code available at https://github.com/AdrienTaylor/Performance-Estimation-Toolbox

In this paper we present a complete iteration complexity analysis of inexact first order Lagrangian and penalty methods for solving cone constrained convex problems that have or may not have optimal Lagrange multipliers that close the duality gap. We first assume the existence of optimal Lagrange multipliers and study primal-dual first order methods based on inexact information and augmented Lagrangian smoothing or Nesterov type smoothing. For inexact (fast) gradient augmented Lagrangian methods we derive a total computational complexity of $\mathcal{O}\left(\frac{1}{\epsilon}\right)$ projections onto a simple primal set in order to attain an $\epsilon-$optimal solution of the conic convex problem. For the inexact fast gradient method combined with Nesterov type smoothing we derive computational complexity $\mathcal{O}\left(\frac{1}{\epsilon^{3/2}}\right)$ projections onto the same set. Then, we assume that optimal Lagrange multipliers for the cone constrained convex problem might not exist, and analyze the fast gradient method for solving penalty reformulations of the problem. For the fast gradient method combined with penalty framework we also derive a total computational complexity of $\mathcal{O}\left(\frac{1}{\epsilon^{3/2}}\right)$ projections onto a simple primal set to attain an $\epsilon-$optimal solution for the original problem., Comment: 35 pages, Technical Report, UPB, March 2017

In this paper, we propose new lower and upper bounds on the linear extension complexity of regular $n$-gons. Our bounds are based on the equivalence between the computation of (i) an extended formulation of size $r$ of a polytope $P$, and (ii) a rank-$r$ nonnegative factorization of a slack matrix of the polytope $P$. The lower bound is based on an improved bound for the rectangle covering number (also known as the boolean rank) of the slack matrix of the $n$-gons. The upper bound is a slight improvement of the result of Fiorini, Rothvoss and Tiwary [Extended Formulations for Polygons, Discrete Comput. Geom. 48(3), pp. 658-668, 2012]. The difference with their result is twofold: (i) our proof uses a purely algebraic argument while Fiorini et al. used a geometric argument, and (ii) we improve the base case allowing us to reduce their upper bound $2 \left\lceil \log_2(n) \right\rceil$ by one when $2^{k-1} < n \leq 2^{k-1}+2^{k-2}$ for some integer $k$. We conjecture that this new upper bound is tight, which is suggested by numerical experiments for small $n$. Moreover, this improved upper bound allows us to close the gap with the best known lower bound for certain regular $n$-gons (namely, $9 \leq n \leq 13$ and $21 \leq n \leq 24$) hence allowing for the first time to determine their extension complexity., Comment: 20 pages, 3 figures. New contribution: improved lower bound for the boolean rank of the slack matrices of n-gons

We show that the exact worst-case performance of fixed-step first-order methods for unconstrained optimization of smooth (possibly strongly) convex functions can be obtained by solving convex programs. Finding the worst-case performance of a black-box first-order method is formulated as an optimization problem over a set of smooth (strongly) convex functions and initial conditions. We develop closed-form necessary and sufficient conditions for smooth (strongly) convex interpolation, which provide a finite representation for those functions. This allows us to reformulate the worst-case performance estimation problem as an equivalent finite dimension-independent semidefinite optimization problem, whose exact solution can be recovered up to numerical precision. Optimal solutions to this performance estimation problem provide both worst-case performance bounds and explicit functions matching them, as our smooth (strongly) convex interpolation procedure is constructive. Our works build on those of Drori and Teboulle in [Math. Prog. 145 (1-2), 2014] who introduced and solved relaxations of the performance estimation problem for smooth convex functions. We apply our approach to different fixed-step first-order methods with several performance criteria, including objective function accuracy and gradient norm. We conjecture several numerically supported worst-case bounds on the performance of the fixed-step gradient, fast gradient and optimized gradient methods, both in the smooth convex and the smooth strongly convex cases, and deduce tight estimates of the optimal step size for the gradient method., Comment: Accepted for publication in Mathematical Programming, DOI 10.1007/s10107-016-1009-3. Performance estimation code: http://perso.uclouvain.be/adrien.taylor/#projects

The exact nonnegative matrix factorization (exact NMF) problem is the following: given an $m$-by-$n$ nonnegative matrix $X$ and a factorization rank $r$, find, if possible, an $m$-by-$r$ nonnegative matrix $W$ and an $r$-by-$n$ nonnegative matrix $H$ such that $X = WH$. In this paper, we propose two heuristics for exact NMF, one inspired from simulated annealing and the other from the greedy randomized adaptive search procedure. We show that these two heuristics are able to compute exact nonnegative factorizations for several classes of nonnegative matrices (namely, linear Euclidean distance matrices, slack matrices, unique-disjointness matrices, and randomly generated matrices) and as such demonstrate their superiority over standard multi-start strategies. We also consider a hybridization between these two heuristics that allows us to combine the advantages of both methods. Finally, we discuss the use of these heuristics to gain insight on the behavior of the nonnegative rank, i.e., the minimum factorization rank such that an exact NMF exists. In particular, we disprove a conjecture on the nonnegative rank of a Kronecker product, propose a new upper bound on the extension complexity of generic $n$-gons and conjecture the exact value of (i) the extension complexity of regular $n$-gons and (ii) the nonnegative rank of a submatrix of the slack matrix of the correlation polytope., Comment: 32 pages, 2 figures, 16 tables

Several practical multi-user multi-carrier communication systems are characterized by a multi-carrier interference channel system model where the interference is treated as noise. For these systems, spectrum optimization is a promising means to mitigate interference. This however corresponds to a challenging nonconvex optimization problem. Existing iterative convex approximation (ICA) methods consist in solving a series of improving convex approximations and are typically implemented in a per-user iterative approach. However they do not take this typical iterative implementation into account in their design. This paper proposes a novel class of iterative approximation methods that focuses explicitly on the per-user iterative implementation, which allows to relax the problem significantly, dropping joint convexity and even convexity requirements for the approximations. A systematic design framework is proposed to construct instances of this novel class, where several new iterative approximation methods are developed with improved per-user convex and nonconvex approximations that are both tighter and simpler to solve (in closed-form). As a result, these novel methods display a much faster convergence speed and require a significantly lower computational cost. Furthermore, a majority of the proposed methods can tackle the issue of getting stuck in bad locally optimal solutions, and hence improve solution quality compared to existing ICA methods., Comment: 33 pages, 7 figures. This work has been submitted for possible publication

Dynamic spectrum management is recognized as a key technique to tackle interference in multi-user multi-carrier communication systems and networks. However existing dynamic spectrum management algorithms may not be suitable when the available computation time and compute power are limited, i.e., when a very fast responsiveness is required. In this paper, we present a new paradigm, theory and algorithm for real-time dynamic spectrum management (RT-DSM) under tight real-time constraints. Specifically, a RT-DSM algorithm can be stopped at any point in time while guaranteeing a feasible and improved solution. This is enabled by the introduction of a novel difference-of-variables (DoV) transformation and problem reformulation, for which a primal coordinate ascent approach is proposed with exact line search via a logarithmicly scaled grid search. The concrete proposed algorithm is referred to as iterative power difference balancing (IPDB). Simulations for different realistic wireline and wireless interference limited systems demonstrate its good performance, low complexity and wide applicability under different configurations., Comment: 14 pages, 9 figures. This work has been submitted to the IEEE for possible publication

This paper addresses a weighted sum rate (WSR) maximization problem for downlink OFDMA aided by a decode-and-forward (DF) relay under a total power constraint. A novel subcarrier-pair based opportunistic DF relaying protocol is proposed. Specifically, user message bits are transmitted in two time slots. A subcarrier in the first slot can be paired with a subcarrier in the second slot for the DF relay-aided transmission to a user. In particular, the source and the relay can transmit simultaneously to implement beamforming at the subcarrier in the second slot. Each unpaired subcarrier in either the first or second slot is used for the source's direct transmission to a user. A benchmark protocol, same as the proposed one except that the transmit beamforming is not used for the relay-aided transmission, is also considered. For each protocol, a polynomial-complexity algorithm is developed to find at least an approximately optimum resource allocation (RA), by using continuous relaxation, the dual method, and Hungarian algorithm. Instrumental to the algorithm design is an elegant definition of optimization variables, motivated by the idea of regarding the unpaired subcarriers as virtual subcarrier pairs in the direct transmission mode. The effectiveness of the RA algorithm and the impact of relay position and total power on the protocols' performance are illustrated by numerical experiments. The proposed protocol always leads to a maximum WSR equal to or greater than that for the benchmark one, and the performance gain of using the proposed one is significant especially when the relay is in close proximity to the source and the total power is low. Theoretical analysis is presented to interpret these observations., Comment: 8 figures, accepted and to be published in IEEE Transactions on Signal Processing. arXiv admin note: text overlap with arXiv:1301.2935

We study the problem of separating audio sources from a single linear mixture. The goal is to find a decomposition of the single channel spectrogram into a sum of individual contributions associated to a certain number of sources. In this paper, we consider an informed source separation problem in which the input spectrogram is partly annotated. We propose a convex formulation that relies on a nuclear norm penalty to induce low rank for the contributions. We show experimentally that solving this model with a simple subgradient method outperforms a previously introduced nonnegative matrix factorization (NMF) technique, both in terms of source separation quality and computation time., Comment: Submitted to ESANN 2013 conference

Approximate matrix factorization techniques with both nonnegativity and orthogonality constraints, referred to as orthogonal nonnegative matrix factorization (ONMF), have been recently introduced and shown to work remarkably well for clustering tasks such as document classification. In this paper, we introduce two new methods to solve ONMF. First, we show athematical equivalence between ONMF and a weighted variant of spherical k-means, from which we derive our first method, a simple EM-like algorithm. This also allows us to determine when ONMF should be preferred to k-means and spherical k-means. Our second method is based on an augmented Lagrangian approach. Standard ONMF algorithms typically enforce nonnegativity for their iterates while trying to achieve orthogonality at the limit (e.g., using a proper penalization term or a suitably chosen search direction). Our method works the opposite way: orthogonality is strictly imposed at each step while nonnegativity is asymptotically obtained, using a quadratic penalty. Finally, we show that the two proposed approaches compare favorably with standard ONMF algorithms on synthetic, text and image data sets., Comment: 17 pages, 8 figures. New numerical experiments (document and synthetic data sets)

Nonnegative matrix factorization (NMF) is a data analysis technique used in a great variety of applications such as text mining, image processing, hyperspectral data analysis, computational biology, and clustering. In this paper, we consider two well-known algorithms designed to solve NMF problems, namely the multiplicative updates of Lee and Seung and the hierarchical alternating least squares of Cichocki et al. We propose a simple way to significantly accelerate these schemes, based on a careful analysis of the computational cost needed at each iteration, while preserving their convergence properties. This acceleration technique can also be applied to other algorithms, which we illustrate on the projected gradient method of Lin. The efficiency of the accelerated algorithms is empirically demonstrated on image and text datasets, and compares favorably with a state-of-the-art alternating nonnegative least squares algorithm., Comment: 17 pages, 10 figures. New Section 4 about the convergence of the accelerated algorithms; Removed Section 5 about efficiency of HALS. Accepted in Neural Computation

Weighted low-rank approximation (WLRA), a dimensionality reduction technique for data analysis, has been successfully used in several applications, such as in collaborative filtering to design recommender systems or in computer vision to recover structure from motion. In this paper, we study the computational complexity of WLRA and prove that it is NP-hard to find an approximate solution, even when a rank-one approximation is sought. Our proofs are based on a reduction from the maximum-edge biclique problem, and apply to strictly positive weights as well as binary weights (the latter corresponding to low-rank matrix approximation with missing data)., Comment: Proof of Lemma 4 (Lemma 3 in v1) has been corrected. Some remarks and comments have been added. Accepted in SIAM Journal on Matrix Analysis and Applications

The nonnegative rank of a nonnegative matrix is the minimum number of nonnegative rank-one factors needed to reconstruct it exactly. The problem of determining this rank and computing the corresponding nonnegative factors is difficult; however it has many potential applications, e.g., in data mining, graph theory and computational geometry. In particular, it can be used to characterize the minimal size of any extended reformulation of a given combinatorial optimization program. In this paper, we introduce and study a related quantity, called the restricted nonnegative rank. We show that computing this quantity is equivalent to a problem in polyhedral combinatorics, and fully characterize its computational complexity. This in turn sheds new light on the nonnegative rank problem, and in particular allows us to provide new improved lower bounds based on its geometric interpretation. We apply these results to slack matrices and linear Euclidean distance matrices and obtain counter-examples to two conjectures of Beasly and Laffey, namely we show that the nonnegative rank of linear Euclidean distance matrices is not necessarily equal to their dimension, and that the rank of a matrix is not always greater than the nonnegative rank of its square.

Nonnegative Matrix Factorization (NMF) is the problem of approximating a nonnegative matrix with the product of two low-rank nonnegative matrices and has been shown to be particularly useful in many applications, e.g., in text mining, image processing, computational biology, etc. In this paper, we explain how algorithms for NMF can be embedded into the framework of multilevel methods in order to accelerate their convergence. This technique can be applied in situations where data admit a good approximate representation in a lower dimensional space through linear transformations preserving nonnegativity. A simple multilevel strategy is described and is experimentally shown to speed up significantly three popular NMF algorithms (alternating nonnegative least squares, multiplicative updates and hierarchical alternating least squares) on several standard image datasets., Comment: 23 pages, 10 figures. Section 6 added discussing limitations of the method. Accepted in Journal of Computational and Applied Mathematics

Nonnegative Matrix Factorization consists in (approximately) factorizing a nonnegative data matrix by the product of two low-rank nonnegative matrices. It has been successfully applied as a data analysis technique in numerous domains, e.g., text mining, image processing, microarray data analysis, collaborative filtering, etc. We introduce a novel approach to solve NMF problems, based on the use of an underapproximation technique, and show its effectiveness to obtain sparse solutions. This approach, based on Lagrangian relaxation, allows the resolution of NMF problems in a recursive fashion. We also prove that the underapproximation problem is NP-hard for any fixed factorization rank, using a reduction of the maximum edge biclique problem in bipartite graphs. We test two variants of our underapproximation approach on several standard image datasets and show that they provide sparse part-based representations with low reconstruction error. Our results are comparable and sometimes superior to those obtained by two standard Sparse Nonnegative Matrix Factorization techniques., Comment: Version 2 removed the section about convex reformulations, which was not central to the development of our main results; added material to the introduction; added a review of previous related work (section 2.3); completely rewritten the last part (section 4) to provide extensive numerical results supporting our claims. Accepted in J. of Pattern Recognition

Nonnegative Matrix Factorization (NMF) is a data analysis technique which allows compression and interpretation of nonnegative data. NMF became widely studied after the publication of the seminal paper by Lee and Seung (Learning the Parts of Objects by Nonnegative Matrix Factorization, Nature, 1999, vol. 401, pp. 788--791), which introduced an algorithm based on Multiplicative Updates (MU). More recently, another class of methods called Hierarchical Alternating Least Squares (HALS) was introduced that seems to be much more efficient in practice. In this paper, we consider the problem of approximating a not necessarily nonnegative matrix with the product of two nonnegative matrices, which we refer to as Nonnegative Factorization (NF); this is the subproblem that HALS methods implicitly try to solve at each iteration. We prove that NF is NP-hard for any fixed factorization rank, using a reduction to the maximum edge biclique problem. We also generalize the multiplicative updates to NF, which allows us to shed some light on the differences between the MU and HALS algorithms for NMF and give an explanation for the better performance of HALS. Finally, we link stationary points of NF with feasible solutions of the biclique problem to obtain a new type of biclique finding algorithm (based on MU) whose iterations have an algorithmic complexity proportional to the number of edges in the graph, and show that it performs better than comparable existing methods.