Efficient Computation of ℓ1 Regularized Estimates in Gaussian Graphical Models

In this article, I propose an effcient algorithm to compute ?₁ regularized maximum likelihood estimates in the Gaussian graphical model. These estimators, recently proposed in an earlier article by Yuan and Lin, conduct parameter estimation and model selection simultaneously and have been shown to enjoy nice properties in both large and finite samples. To compute the estimates, however, can be very challenging in practice because of the high dimensionality and positive definiteness constraint on the covariance matrix. Taking advantage of the recent advance in semidefinite programming, Yuan and Lin suggested a sophisticated interior-point algorithm to solve the optimization problem. Although a polynomial time algorithm, the optimization technique is known not to be scalable for high-dimensional problems. Alternatively, this article shows that the estimates can be computed by iteratively solving a sequence of ?₁ regularized quadratic programs. By effectively exploiting the sparsity of the graphical structure, I propose a new algorithm that can be applied to problems of larger scale. When combined with a path-following strategy, the new algorithm can be used to efficiently approximate the entire solution path of the ?₁ regularized maximum likelihood estimates, which also facilitates the choice of tuning parameter. I demonstrate the efficacy and usefulness of the proposed algorithm on a few simulations and real datasets.     

Regularizers for structured sparsity

We study the problem of learning a sparse linear regression vector under additional conditions on the structure of its sparsity pattern. This problem is relevant in machine learning, statistics and signal processing. It is well known that a linear regression can benefit from knowledge that the underlying regression vector is sparse. The combinatorial problem of selecting the nonzero components of this vector can be “relaxed” by regularizing the squared error with a convex penalty function like the ?₁ norm. However, in many applications, additional conditions on the structure of the regression vector and its sparsity pattern are available. Incorporating this information into the learning method may lead to a significant decrease of the estimation error. In this paper, we present a family of convex penalty functions, which encode prior knowledge on the structure of the vector formed by the absolute values of the regression coefficients. This family subsumes the ?₁ norm and is flexible enough to include different models of sparsity patterns, which are of practical and theoretical importance. We establish the basic properties of these penalty functions and discuss some examples where they can be computed explicitly. Moreover, we present a convergent optimization algorithm for solving regularized least squares with these penalty functions. Numerical simulations highlight the benefit of structured sparsity and the advantage offered by our approach over the Lasso method and other related methods.     

Computing Syzygies by Faug��re��s F5 algorithm

In this paper, we introduce a new algorithm for computing a set of generators for the syzygies on a sequence of polynomials. For this, we extend a given sequence of polynomials to a Gr?bner basis using Faugère??s F₅ algorithm (A new efficient algorithm for computing Gr?bner bases without reduction to zero (F ₅). ISSAC, ACM Press, pp 75?C83, 2002). We show then that if we keep all the reductions to zero during this computation, then at termination (by adding principal syzygies) we obtain a basis for the module of syzygies on the input polynomials. We have implemented our algorithm in the computer algebra system Magma, and we evaluate its performance via some examples.     

Matrix stretching for sparse least squares problems

For linear least squares problems min _x ‖Ax − b‖₂, where A is sparse except for a few dense rows, a straightforward application of Cholesky or QR factorization will lead to catastrophic fill in the factor R. We consider handling such problems by a matrix stretching technique, where the dense rows are split into several more sparse rows. We develop both a recursive binary splitting algorithm and a more general splitting method. We show that for both schemes the stretched problem has the same set of solutions as the original least squares problem. Further, the condition number of the stretched problem differs from that of the original by only a modest factor, and hence the approach is numerically stable. Experimental results from applying the recursive binary scheme to a set of modified matrices from the Harwell‐Boeing collection are given. We conclude that when A has a small number of dense rows relative to its dimension, there is a significant gain in sparsity of the factor R. A crude estimate of the optimal number of splits is obtained by analysing a simple model problem. Copyright © 2000 John Wiley & Sons, Ltd.     

A sparsity preserving stochastic gradient methods for sparse regression

We propose a new stochastic first-order algorithm for solving sparse regression problems. In each iteration, our algorithm utilizes a stochastic oracle of the subgradient of the objective function. Our algorithm is based on a stochastic version of the estimate sequence technique introduced by Nesterov (Introductory lectures on convex optimization: a basic course, Kluwer, Amsterdam, 2003). The convergence rate of our algorithm depends continuously on the noise level of the gradient. In particular, in the limiting case of noiseless gradient, the convergence rate of our algorithm is the same as that of optimal deterministic gradient algorithms. We also establish some large deviation properties of our algorithm. Unlike existing stochastic gradient methods with optimal convergence rates, our algorithm has the advantage of readily enforcing sparsity at all iterations, which is a critical property for applications of sparse regressions.     

An efficient algorithm for sparse null space basis problem using ABS methods

We propose a way to use the Markowitz pivot selection criterion for choosing the parameters of the extended ABS class of algorithms to present an effective algorithm for generating sparse null space bases. We explain in detail an efficient implementation of the algorithm, making use of the special MATLAB 7.0 functions for sparse matrix operations and the inherent efficiency of the ANSI C programming language. We then compare our proposed algorithm with an implementation of an efficient algorithm proposed by Coleman and Pothen with respect to the computing time and the accuracy and the sparsity of the generated null space bases. Our extensive numerical results, using coefficient matrices of linear programming problems from the NETLIB set of test problems show the competitiveness of our implemented algorithm.     

One-Step Estimator Paths for Concave Regularization

The statistics literature of the past 15 years has established many favorable properties for sparse diminishing-bias regularization: techniques that can roughly be understood as providing estimation under penalty functions spanning the range of concavity between ?₀ and ?₁ norms. However, lasso ?₁-regularized estimation remains the standard tool for industrial Big Data applications because of its minimal computational cost and the presence of easy-to-apply rules for penalty selection. In response, this article proposes a simple new algorithm framework that requires no more computation than a lasso path: the path of one-step estimators (POSE) does ?₁ penalized regression estimation on a grid of decreasing penalties, but adapts coefficient-specific weights to decrease as a function of the coefficient estimated in the previous path step. This provides sparse diminishing-bias regularization at no extra cost over the fastest lasso algorithms. Moreover, our gamma lasso implementation of POSE is accompanied by a reliable heuristic for the fit degrees of freedom, so that standard information criteria can be applied in penalty selection. We also provide novel results on the distance between weighted-?₁ and ?₀ penalized predictors; this allows us to build intuition about POSE and other diminishing-bias regularization schemes. The methods and results are illustrated in extensive simulations and in application of logistic regression to evaluating the performance of hockey players. Supplementary materials for this article are available online.     

Asymptotically Optimal Algorithms for Job Shop Scheduling and Packet Routing

We propose asymptotically optimal algorithms for the job shop scheduling and packet routing problems. We propose a fluid relaxation for the job shop scheduling problem in which we replace discrete jobs with the flow of a continuous fluid. We compute an optimal solution of the fluid relaxation in closed form, obtain a lower bound C_max to the job shop scheduling problem, and construct a feasible schedule from the fluid relaxation with objective value at most where the constant in the O( · ) notation is independent of the number of jobs, but it depends on the processing time of the jobs, thus producing an asymptotically optimal schedule as the total number of jobs tends to infinity. If the initially present jobs increase proportionally, then our algorithm produces a schedule with value at most C_max + O(1). For the packet routing problem with fixed paths the previous algorithm applies directly. For the general packet routing problem we propose a linear programming relaxation that provides a lower bound C_max and an asymptotically optimal algorithm that uses the optimal solution of the relaxation with objective value at most Unlike asymptotically optimal algorithms that rely on probabilistic assumptions, our proposed algorithms make no probabilistic assumptions and they are asymptotically optimal for all instances with a large number of jobs (packets). In computational experiments our algorithms produce schedules which are within 1% of optimality even for moderately sized problems.     

Minimizing Sum of Truncated Convex Functions and Its Applications

In this article, we study a class of problems where the sum of truncated convex functions is minimized. In statistical applications, they are commonly encountered when ?₀-penalized models are fitted and usually lead to NP-Hard non-convex optimization problems. In this article, we propose a general algorithm for the global minimizer in low-dimensional settings. We also extend the algorithm to high-dimensional settings, where an approximate solution can be found efficiently. We introduce several applications where the sum of truncated convex functions is used, compare our proposed algorithm with other existing algorithms in simulation studies, and show its utility in edge-preserving image restoration on real data.     

Simple algorithms for approximating all roots of a polynomial with real roots

We present a simple algorithm for approximating all roots of a polynomial p(x) when it has only real roots. The algorithm is based on some interesting properties of the polynomials appearing in the Extended Euclidean Scheme for p(x) and p′(x). For example, it turns out that these polynomials are orthogonal; as a consequence, we are able to limit the precision required by our algorithm in intermediate steps. A parallel implementation of this algorithm yields a P-uniform NC₂ circuit, and the bit complexity of its sequential implementation is within a polylog factor of the bit complexity of the best known algorithm for the problem.     

SPAI Preconditioners for HPC Applications

Iterative methods to solve systems of linear equations Ax = b usually require preconditioners M to speed convergence. But the calculation of many preconditioners can be notoriously sequential. The sparse approximate inverse preconditioner (SPAI) has particular potential for high performance computing [1]. We have ported the SPAI algorithm to graphical processing units (GPUs) within NVIDIA's CUSP library [2] for sparse linear algebra. GPUs perform well on dense linear algebra problems where data resides for long periods on the device. Since the underlying minimization problems are independent, they are mapped to separate thread-blocks, and an optimized QR algorithm implemented using NVIDIA's CUBLAS library is employed on each. Traditionally the challenge has been to determine a sparsity pattern Sp( M ) of the preconditioner dynamically [3], which reduces the condition number of MA to a degree where a preconditioned iterative solver such as GMRES becomes computationally viable. Due to the extremely high performance of the GPU, it is possible to consider initial sparsity patterns much denser than have been previously considered. We therefore consider a fixed sparsity patterns to simplify the GPU implementation. We evaluate the performance of the resulting preconditioner on a standard set of sparse matrices, and compare SPAI to other preconditioners. (© 2012 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim)     

Compressed sensing with local structure: Uniform recovery guarantees for the sparsity in levels class

In compressed sensing, it is often desirable to consider signals possessing additional structure beyond sparsity. One such structured signal model – which forms the focus of this paper – is the local sparsity in levels class. This class has recently found applications in problems such as compressive imaging, multi-sensor acquisition systems and sparse regularization in inverse problems. In this paper we present uniform recovery guarantees for this class when the measurement matrix corresponds to a subsampled isometry. We do this by establishing a variant of the standard restricted isometry property for sparse in levels vectors, known as the restricted isometry property in levels. Interestingly, besides the usual log factors, our uniform recovery guarantees are simpler and less stringent than existing nonuniform recovery guarantees. For the particular case of discrete Fourier sampling with Haar wavelet sparsity, a corollary of our main theorem yields a new recovery guarantee which improves over the current state-of-the-art.     

An iterative algorithm for sparse and constrained recovery with applications to divergence-free current reconstructions in magneto-encephalography

We propose an iterative algorithm for the minimization of a ? ₁-norm penalized least squares functional, under additional linear constraints. The algorithm is fully explicit: it uses only matrix multiplications with the three matrices present in the problem (in the linear constraint, in the data misfit part and in the penalty term of the functional). None of the three matrices must be invertible. Convergence is proven in a finite-dimensional setting. We apply the algorithm to a synthetic problem in magneto-encephalography where it is used for the reconstruction of divergence-free current densities subject to a sparsity promoting penalty on the wavelet coefficients of the current densities. We discuss the effects of imposing zero divergence and of imposing joint sparsity (of the vector components of the current density) on the current density reconstruction.     

An Iterative Sparse-Group Lasso

In high-dimensional supervised learning problems, sparsity constraints in the solution often lead to better performance and interpretability of the results. For problems in which covariates are grouped and sparse structure are desired, both on group and within group levels, the sparse-group lasso (SGL) regularization method has proved to be very efficient. Under its simplest formulation, the solution provided by this method depends on two weight parameters that control the penalization on the coefficients. Selecting these weight parameters represents a major challenge. In most of the applications of the SGL, this problem is left aside, and the parameters are either fixed based on a prior information about the data, or chosen to minimize some error function in a grid of possible values. However, an appropriate choice of the parameters deserves more attention, considering that it plays a key role in the structure and interpretation of the solution. In this sense, we present a gradient-free coordinate descent algorithm that automatically selects the regularization parameters of the SGL. We focus on a more general formulation of this problem, which also includes individual penalizations for each group. The advantages of our approach are illustrated using both real and synthetic datasets. Supplementary materials for this article are available online.     

Linear time-periodic dynamical systems: an H2 analysis and a model reduction framework

Linear time-periodic (LTP) dynamical systems frequently appear in the modelling of phenomena related to fluid dynamics, electronic circuits and structural mechanics via linearization centred around known periodic orbits of nonlinear models. Such LTP systems can reach orders that make repeated simulation or other necessary analysis prohibitive, motivating the need for model reduction. We develop here an algorithmic framework for constructing reduced models that retains the LTP structure of the original LTP system. Our approach generalizes optimal approaches that have been established previously for linear time-invariant (LTI) model reduction problems. We employ an extension of the usual H₂ Hardy space defined for the LTI setting to time-periodic systems and within this broader framework develop an a posteriori error bound expressible in terms of related LTI systems. Optimization of this bound motivates our algorithm. We illustrate the success of our method on three numerical examples.     

Consistent learning by composite proximal thresholding

We investigate the modeling and the numerical solution of machine learning problems with prediction functions which are linear combinations of elements of a possibly infinite dictionary of functions. We propose a novel flexible composite regularization model, which makes it possible to incorporate various priors on the coefficients of the prediction function, including sparsity and hard constraints. We show that the estimators obtained by minimizing the regularized empirical risk are consistent in a statistical sense, and we design an error-tolerant composite proximal thresholding algorithm for computing such estimators. New results on the asymptotic behavior of the proximal forward–backward splitting method are derived and exploited to establish the convergence properties of the proposed algorithm. In particular, our method features a o(1 / m) convergence rate in objective values.     

Low‐rank approximation of tensors via sparse optimization

The goal of this paper is to find a low‐rank approximation for a given nth tensor. Specifically, we give a computable strategy on calculating the rank of a given tensor, based on approximating the solution to an NP‐hard problem. In this paper, we formulate a sparse optimization problem via an l₁‐regularization to find a low‐rank approximation of tensors. To solve this sparse optimization problem, we propose a rescaling algorithm of the proximal alternating minimization and study the theoretical convergence of this algorithm. Furthermore, we discuss the probabilistic consistency of the sparsity result and suggest a way to choose the regularization parameter for practical computation. In the simulation experiments, the performance of our algorithm supports that our method provides an efficient estimate on the number of rank‐one tensor components in a given tensor. Moreover, this algorithm is also applied to surveillance videos for low‐rank approximation.     

Approximation schemes for parallel machine scheduling with availability constraints

We investigate the problems of scheduling n weighted jobs to m parallel machines with availability constraints. We consider two different models of availability constraints: the preventive model, in which the unavailability is due to preventive machine maintenance, and the fixed job model, in which the unavailability is due to a priori assignment of some of the n jobs to certain machines at certain times. Both models have applications such as turnaround scheduling or overlay computing. In both models, the objective is to minimize the total weighted completion time. We assume that m is a constant, and that the jobs are non-resumable.For the preventive model, it has been shown that there is no approximation algorithm if all machines have unavailable intervals even if w_i=p_i for all jobs. In this paper, we assume that there is one machine that is permanently available and that the processing time of each job is equal to its weight for all jobs. We develop the first polynomial-time approximation scheme (PTAS) when there is a constant number of unavailable intervals. One main feature of our algorithm is that the classification of large and small jobs is with respect to each individual interval, and thus not fixed. This classification allows us (1) to enumerate the assignments of large jobs efficiently; and (2) to move small jobs around without increasing the objective value too much, and thus derive our PTAS. Next, we show that there is no fully polynomial-time approximation scheme (FPTAS) in this case unless P=NP.For the fixed job model, it has been shown that if job weights are arbitrary then there is no constant approximation for a single machine with 2 fixed jobs or for two machines with one fixed job on each machine, unless P=NP. In this paper, we assume that the weight of a job is the same as its processing time for all jobs. We show that the PTAS for the preventive model can be extended to solve this problem when the number of fixed jobs and the number of machines are both constants.     

Properties of a Class of Nonlinear Transformations Over Euclidean Jordan Algebras with Applications to Complementarity Problems

For any function φ from ?^r to ?^r, Tao and Gowda [Math. Oper. Res., 30 (2005), pp. 985–1004] introduced a corresponding nonlinear transformation R_φ over a Euclidean Jordan algebra (which is called a relaxation transformation) and established some useful relations between φ and R_φ. In this paper, we further investigate some interconnections between properties of φ and properties of R_φ, including the properties of continuity, (local) Lipschitz continuity, directional differentiability, (continuous) differentiability, semismoothness, monotonicity, the P₀-property, and the uniform P-property. As an application, we investigate the symmetric cone complementarity problem with a relaxation transformation. A property of the solution set of this class of problems is given. We also investigate a smoothing algorithm for solving this class of problems and show that the algorithm is globally convergent under an assumption that the solution set of the problem concerned is nonempty.     

Nonlinear Least Squares in ℝ N

Recent research and new paradigms in mathematics, engineering, and science assume nonlinear signal models of the form ?=∪ _i∈I V _i consisting of a union of subspaces V _i instead of a single subspace ?=V. These models have been used in sampling and reconstruction of signals with finite rate of innovation, the Generalized Principle Component Analysis and the subspace segmentation problem in computer vision, and problems related to sparsity, compressed sensing, and dictionary design. In this paper, we develop an algorithm that searches for the best nonlinear model of the form ?=∪ _i=1 ^l V _i ?? ^N that is optimally compatible with a set of observations ?={f ₁,…,f _m }?? ^N . When l=1 this becomes the classical least squares optimization. Thus, this problem is a nonlinear version of the least squares problem. We test our algorithm on synthetic data as well as images.