CP decomposition for tensors via alternating least squares with QR decomposition

The CP tensor decomposition is used in applications such as machine learning and signal processing to discover latent low-rank structure in multidimensional data. Computing a CP decomposition via an alternating least squares (ALS) method reduces the problem to several linear least squares problems. The standard way to solve these linear least squares subproblems is to use the normal equations, which inherit special tensor structure that can be exploited for computational efficiency. However, the normal equations are sensitive to numerical ill-conditioning, which can compromise the results of the decomposition. In this paper, we develop versions of the CP-ALS algorithm using the QR decomposition and the singular value decomposition, which are more numerically stable than the normal equations, to solve the linear least squares problems. Our algorithms utilize the tensor structure of the CP-ALS subproblems efficiently, have the same complexity as the standard CP-ALS algorithm when the input is dense and the rank is small, and are shown via examples to produce more stable results when ill-conditioning is present. Our MATLAB implementation achieves the same running time as the standard algorithm for small ranks, and we show that the new methods can obtain lower approximation error.     

A modified nearly exact method for solving low-rank trust region subproblem

In this paper, we first discuss how the nearly exact (NE) method proposed by Moré and Sorensen [14] for solving trust region (TR) subproblems can be modified to solve large-scale “low-rank” TR subproblems efficiently. Our modified algorithm completely avoids computation of Cholesky factorizations by instead relying primarily on the Sherman–Morrison–Woodbury formula for computing inverses of “diagonal plus low-rank” type matrices. We also implement a specific version of the modified log-barrier (MLB) algorithm proposed by Polyak [17] where the generated log-barrier subproblems are solved by a trust region method. The corresponding direction finding TR subproblems are of the low-rank type and are then solved by our modified NE method. We finally discuss the computational results of our implementation of the MLB method and its comparison with a version of LANCELOT [5] based on a collection extracted from CUTEr [12] of nonlinear programming problems with simple bound constraints.      

Overlapping Domain Decomposition Algorithms for General Sparse Matrices

Domain decomposition methods for finite element problems using a partition based on the underlying finite element mesh have been extensively studied. In this paper, we discuss algebraic extensions of the class of overlapping domain decomposition algorithms for general sparse matrices. The subproblems are created with an overlapping partition of the graph corresponding to the sparsity structure of the matrix. These algebraic domain decomposition methods are especially useful for unstructured mesh problems. We also discuss some difficulties encountered in the algebraic extension, particularly the issues related to the coarse solver.     

Nesterov acceleration of alternating least squares for canonical tensor decomposition: Momentum step size selection and restart mechanisms

We present Nesterov‐type acceleration techniques for alternating least squares (ALS) methods applied to canonical tensor decomposition. While Nesterov acceleration turns gradient descent into an optimal first‐order method for convex problems by adding a momentum term with a specific weight sequence, a direct application of this method and weight sequence to ALS results in erratic convergence behavior. This is so because ALS is accelerated instead of gradient descent for our nonconvex problem. Instead, we consider various restart mechanisms and suitable choices of momentum weights that enable effective acceleration. Our extensive empirical results show that the Nesterov‐accelerated ALS methods with restart can be dramatically more efficient than the stand‐alone ALS or Nesterov's accelerated gradient methods, when problems are ill‐conditioned or accurate solutions are desired. The resulting methods perform competitively with or superior to existing acceleration methods for ALS, including ALS acceleration by nonlinear conjugate gradient, nonlinear generalized minimal residual method, or limited‐memory Broyden‐Fletcher‐Goldfarb‐Shanno, and additionally enjoy the benefit of being much easier to implement. We also compare with Nesterov‐type updates where the momentum weight is determined by a line search (LS), which are equivalent or closely related to existing LS methods for ALS. On a large and ill‐conditioned 71×1,000×900 tensor consisting of readings from chemical sensors to track hazardous gases, the restarted Nesterov‐ALS method shows desirable robustness properties and outperforms any of the existing methods we compare with by a large factor. There is clear potential for extending our Nesterov‐type acceleration approach to accelerating other optimization algorithms than ALS applied to other nonconvex problems, such as Tucker tensor decomposition.     

A nonlinearly preconditioned conjugate gradient algorithm for rank‐R canonical tensor approximation

Alternating least squares (ALS) is often considered the workhorse algorithm for computing the rank‐R canonical tensor approximation, but for certain problems, its convergence can be very slow. The nonlinear conjugate gradient (NCG) method was recently proposed as an alternative to ALS, but the results indicated that NCG is usually not faster than ALS. To improve the convergence speed of NCG, we consider a nonlinearly preconditioned NCG (PNCG) algorithm for computing the rank‐R canonical tensor decomposition. Our approach uses ALS as a nonlinear preconditioner in the NCG algorithm. Alternatively, NCG can be viewed as an acceleration process for ALS. We demonstrate numerically that the convergence acceleration mechanism in PNCG often leads to important pay‐offs for difficult tensor decomposition problems, with convergence that is significantly faster and more robust than for the stand‐alone NCG or ALS algorithms. We consider several approaches for incorporating the nonlinear preconditioner into the NCG algorithm that have been described in the literature previously and have met with success in certain application areas. However, it appears that the nonlinearly PNCG approach has received relatively little attention in the broader community and remains underexplored both theoretically and experimentally. Thus, this paper serves several additional functions, by providing in one place a concise overview of several PNCG variants and their properties that have only been described in a few places scattered throughout the literature, by systematically comparing the performance of these PNCG variants for the tensor decomposition problem, and by drawing further attention to the usefulness of nonlinearly PNCG as a general tool. In addition, we briefly discuss the convergence of the PNCG algorithm. In particular, we obtain a new convergence result for one of the PNCG variants under suitable conditions, building on known convergence results for non‐preconditioned NCG. Copyright © 2014 John Wiley & Sons, Ltd.     

A sequential cutting plane algorithm for solving convex NLP problems

In this paper we look at a new algorithm for solving convex nonlinear programming optimization problems. The algorithm is a cutting plane-based method, where the sizes of the subproblems remain fixed, thus avoiding the issue with constantly growing subproblems we have for the classical Kelley’s cutting plane algorithm. Initial numerical experiments indicate that the algorithm is considerably faster than Kelley’s cutting plane algorithm and also competitive with existing nonlinear programming algorithms.     

Nonnegative non-redundant tensor decomposition

Nonnegative tensor decomposition allows us to analyze data in their ‘native’ form and to present results in the form of the sum of rank-1 tensors that does not nullify any parts of the factors. In this paper, we propose the geometrical structure of a basis vector frame for sum-of-rank-1 type decomposition of real-valued nonnegative tensors. The decomposition we propose reinterprets the orthogonality property of the singularvectors of matrices as a geometric constraint on the rank-1 matrix bases which leads to a geometrically constrained singularvector frame. Relaxing the orthogonality requirement, we developed a set of structured-bases that can be utilized to decompose any tensor into a similar constrained sum-of-rank-1 decomposition. The proposed approach is essentially a reparametrization and gives us an upper bound of the rank for tensors. At first, we describe the general case of tensor decomposition and then extend it to its nonnegative form. At the end of this paper, we show numerical results which conform to the proposed tensor model and utilize it for nonnegative data decomposition.     

Medium-term hydro-thermal coordination via hydro and thermal subsystem decomposition

Summary This paper addresses the medium-term hydro-thermal coordination problem in an electric energy system. That is, the problem of finding the energy production of every power plant (hydro or thermal) in every subperiod of a given planning period, so that the customer load is supplied at minimum cost. The planning horizon is typically one to two months and the first week of this planning period is modeled in detail. The solution method proposed decomposes the problem in two subproblems corresponding to the hydro and thermal subsystems. These two subproblems are coordinated using a coordinating function for every subperiod. The coordinating function of a given subperiod expresses total production cost in that subperiod as a function of the total hydro production in that subperiod. The decomposition proposed makes it possible to use specialized algorithms to solve the hydro and thermal subproblems. This results in a very efficient computational procedure. From an experimental point of view the coordinating mechanism is robust. A case study is provided. It considers 61 thermal plants, a hydro system including 8 cascaded hydro plants and a 48 subperiods planning period.     

Analysis of stochastic problem decomposition algorithms in computational grids

Stochastic programming usually represents uncertainty discretely by means of a scenario tree. This representation leads to an exponential growth of the size of stochastic mathematical problems when better accuracy is needed. Trying to solve the problem as a whole, considering all scenarios together, yields to huge memory requirements that surpass the capabilities of current computers. Thus, decomposition algorithms are employed to divide the problem into several smaller subproblems and to coordinate their solution in order to obtain the global optimum. This paper analyzes several decomposition strategies based on the classical Benders decomposition algorithm, and applies them in the emerging computational grid environments. Most decomposition algorithms are not able to take full advantage of all the computing power available in a grid system because of unavoidable dependencies inherent to the algorithms. However, a special decomposition method presented in this paper aims at reducing dependency among subproblems, to the point where all the subproblems can be sent simultaneously to the grid. All algorithms have been tested in a grid system, measuring execution times required to solve standard optimization problems and a real-size hydrothermal coordination problem. Numerical results are shown to confirm that this new method outperforms the classical ones when used in grid computing environments.     

The computation of the square roots of circulant matrices

In the first part of this paper, we investigate the reduced forms of circulant matrices and quasi-skew circulant matrices. By using their properties we present two efficient algorithms to compute the square roots of circulant matrices and quasi-skew circulant matrices, respectively. Those methods are faster than the traditional algorithm which is based on the Schur decomposition. In the second part, we further consider circulant H-matrices with positive diagonal entries and develop two algorithms for computing their principal square roots. Those two algorithms have the common advantage that is they only need matrix-matrix multiplications in their iterative sequences, an operation which can be done very efficiently on modern high performance computers.     

A heuristic lagrangean algorithm for the capacitated plant location problem

Lagrangean techniques have been widely applied to the uncapacitated plant location problem, and in some cases they have proven to be successfull even when capacitated problems with additional constraints are taken into account. In our paper we study the application of these techniques to the capacitated plant location problem when the model considered is a pure integer one. Several lagrangean decompositions are considered and for some of them heuristic algorithms have been designed to solve the resulting lagrangean subproblems, the heuristics consisting of a two phase procedure. The first (location phase) defines a set of multipliers from the analysis of the dual LP relaxation, and makes a choice of the plants considering the resulting subproblems as a particular case of the general assignment problems. Several heuristics have been studied for this second phase, based either on a decomposition of knapsack type subproblems through a definition of a set of penalties, or of looking into the duality gap and trying to reduce it. Computational experience is reported.     

On vector and matrix median computation

The aim of this paper is to gain more insight into vector and matrix medians and to investigate algorithms to compute them. We prove relations between vector and matrix means and medians, particularly regarding the classical structure tensor. Moreover, we examine matrix medians corresponding to different unitarily invariant matrix norms for the case of symmetric 2×2 matrices, which frequently arise in image processing. Our findings are explained and illustrated by numerical examples. To solve the corresponding minimization problems, we propose several algorithms. Existing approaches include Weiszfeld’s algorithm for the computation of ?₂ vector medians and semi-definite programming, in particular, second order cone programming, which has been used for matrix median computation. In this paper, we adapt Weiszfeld’s algorithm for our setting and show that also two splitting methods, namely the alternating direction method of multipliers and the parallel proximal algorithm, can be applied for generalized vector and matrix median computations. Besides, we compare the performance of these algorithms numerically and apply them within local median filters.     

An Improved Gradient Projection-based Decomposition Technique for Support Vector Machines

In this paper we propose some improvements to a recent decomposition technique for the large quadratic program arising in training support vector machines. As standard decomposition approaches, the technique we consider is based on the idea to optimize, at each iteration, a subset of the variables through the solution of a quadratic programming subproblem. The innovative features of this approach consist in using a very effective gradient projection method for the inner subproblems and a special rule for selecting the variables to be optimized at each step. These features allow to obtain promising performance by decomposing the problem into few large subproblems instead of many small subproblems as usually done by other decomposition schemes. We improve this technique by introducing a new inner solver and a simple strategy for reducing the computational cost of each iteration. We evaluate the effectiveness of these improvements by solving large-scale benchmark problems and by comparison with a widely used decomposition package.     

Prices stabilization for inexact unit-commitment problems

A widespread and successful approach to tackle unit-commitment problems is constraint decomposition: by dualizing the linking constraints, the large-scale nonconvex problem decomposes into smaller independent subproblems. The dual problem consists then in finding the best Lagrangian multiplier (the optimal “price”); it is solved by a convex nonsmooth optimization method. Realistic modeling of technical production constraints makes the subproblems themselves difficult to solve exactly. Nonsmooth optimization algorithms can cope with inexact solutions of the subproblems. In this case however, we observe that the computed dual solutions show a noisy and unstable behaviour, that could prevent their use as price indicators. In this paper, we present a simple and easy-to-implement way to stabilize dual optimal solutions, by penalizing the noisy behaviour of the prices in the dual objective. After studying the impact of a general stabilization term on the model and the resolution scheme, we focus on the penalization by discrete total variation, showing the consistency of the approach. We illustrate our stabilization on a synthetic example, and real-life problems from EDF (the French Electricity Board).     

Nonnegative singular control systems using the Drazin projector

In this work we study conditions for guaranteeing the nonnegativity of a discrete-time singular control system. A first approach can be found in the literature for general systems, using the whole coefficient matrices. Also, the particular case of matrices of index 1 has been treated by using a block decomposition and the group-projector of the matrix that gives the singularity to the system. In order to complete this study, an analysis of the nonnegativity of a singular control system for matrices having arbitrary index is done by means of the core–nilpotent decomposition. This technique allows us to reduce the size of the original matrices, improving the results where the whole coefficients are involved.     

A General Multiresolution Method for Fitting Functions on the Sphere

In [7], Lyche and Schumaker have described a method for fitting functions of class C ¹ on the sphere which is based on tensor products of quadratic polynomial splines and trigonometric splines of order three associated with uniform knots. In this paper, we present a multiresolution method leading to C ²-functions on the sphere, using tensor products of polynomial and trigonometric splines of odd order with arbitrary simple knot sequences. We determine the decomposition and reconstruction matrices corresponding to the polynomial and trigonometric spline spaces. We describe the general tensor product decomposition and reconstruction algorithms in matrix form which are convenient for the compression of surfaces. We give the different steps of the computer implementation of these algorithms and, finally, we present a test example.     

Nonsmooth Equation Based BFGS Method for Solving KKT Systems in Mathematical Programming

In this paper, we present a BFGS method for solving a KKT system in mathematical programming, based on a nonsmooth equation reformulation of the KKT system. We split successively the nonsmooth equation into equivalent equations with a particular structure. Based on the splitting, we develop a BFGS method in which the subproblems are systems of linear equations with symmetric and positive-definite coefficient matrices. A suitable line search is introduced under which the generated iterates exhibit an approximate norm descent property. The method is well defined and, under suitable conditions, converges to a KKT point globally and superlinearly without any convexity assumption on the problem.     

Parallel bundle-based decomposition for large-scale structured mathematical programming problems

In this paper, we present parallel bundle-based decomposition algorithms to solve a class of structured large-scale convex optimization problems. An example in this class of problems is the block-angular linear programming problem. By dualizing, we transform the original problem to an unconstrained nonsmooth concave optimization problem which is in turn solved by using a modified bundle method. Further, this dual problem consists of a collection of smaller independent subproblems which give rise to the parallel algorithms. We discuss the implementation on the CRYSTAL multi-computer. Finally, we present computational experience with block-angular linear programming problems and observe that more than 70% efficiency can be obtained using up to eleven processors for one group of test problems, and more than 60% efficiency can be obtained for relatively smaller problems using up to five processors for another group of problems.     

Decomposition with branch-and-cut approaches for two-stage stochastic mixed-integer programming

Decomposition has proved to be one of the more effective tools for the solution of large-scale problems, especially those arising in stochastic programming. A decomposition method with wide applicability is Benders' decomposition, which has been applied to both stochastic programming as well as integer programming problems. However, this method of decomposition relies on convexity of the value function of linear programming subproblems. This paper is devoted to a class of problems in which the second-stage subproblem(s) may impose integer restrictions on some variables. The value function of such integer subproblem(s) is not convex, and new approaches must be designed. In this paper, we discuss alternative decomposition methods in which the second-stage integer subproblems are solved using branch-and-cut methods. One of the main advantages of our decomposition scheme is that Stochastic Mixed-Integer Programming (SMIP) problems can be solved by dividing a large problem into smaller MIP subproblems that can be solved in parallel. This paper lays the foundation for such decomposition methods for two-stage stochastic mixed-integer programs.     

A comparison between the complex symmetric based and classical computation of the singular value decomposition of normal matrices

An algorithm for computing the singular value decomposition of normal matrices using intermediate complex symmetric matrices is proposed. This algorithm, as most eigenvalue and singular value algorithms, consists of two steps. It is based on combining the unitarily equivalence of normal matrices to complex symmetric tridiagonal form with the symmetric singular value decomposition of complex symmetric matrices. Numerical experiments are included comparing several algorithms, with respect to speed and accuracy, for computing the symmetric singular value decomposition (also known as the Takagi factorization). Next we compare the novel approach with the classical Golub-Kahan method for computing the singular value decomposition of normal matrices: it is faster, consumes less memory, but on the other hand the results are significantly less accurate.