Parallel solution of sparse linear least squares problems on distributed-memory multiprocessors

A wide range of problems in computer graphics and vision can be formulated as sparse least squares problems. For example, Laplacian mesh deformation, Least Squares Conformal Maps, Poisson image editing, and as-rigid-as-possible (ARAP) image warping involve solving a linear or non-linear sparse least squares problem. High performance is crucial in many of these applications for interactive user feedback. For these applications, we show that the matrices produced by factorization methods such as QR have a special structure: the off-diagonal blocks are low-rank. We leverage this property to produce a fast sparse approximate QR factorization, spaQR, for these matrices in near-linear time. In our benchmarks, spaQR shows up to 57% improvement over solving the normal equations using Cholesky and 63% improvement over a standard preconditioner with Conjugate Gradient Least Squares (CGLS).

In this paper, we describe scalable parallel algorithms for symmetric sparse matrix factorization, analyze their performance and scalability, and present experimental results for up to 1,024 processors on a Gray T3D parallel computer. Through our analysis and experimental results, we demonstrate that our algorithms substantially improve the state of the art in parallel direct solution of sparse linear systems-both in terms of scalability and overall performance. It is a well known fact that dense matrix factorization scales well and can be implemented efficiently on parallel computers. In this paper, we present the first algorithms to factor a wide class of sparse matrices (including those arising from two- and three-dimensional finite element problems) that are asymptotically as scalable as dense matrix factorization algorithms on a variety of parallel architectures. Our algorithms incur less communication overhead and are more scalable than any previously known parallel formulation of sparse matrix factorization. Although, in this paper, we discuss Cholesky factorization of symmetric positive definite matrices, the algorithms can be adapted for solving sparse linear least squares problems and for Gaussian elimination of diagonally dominant matrices that are almost symmetric in structure. An implementation of one of our sparse Cholesky factorization algorithms delivers up to 20 GFlops on a Gray T3D for medium-size structural engineering and linear programming problems. To the best of our knowledge, this is the highest performance ever obtained for sparse Cholesky factorization on any supercomputer.

There is often a trade-off between preserving sparsity and numerical stability in sparse matrix factorizations. In applications like the direct solution of Equality Constrained Least Squares problem, the accurate detection of the rank of a large and sparse constraint matrix is a key issue. Column pivoting is not suitable for distributed memory machines because it forces the program into a lock-step mode, preventing any overlapping of computations. So factorization algorithms on such machines need to use a reliable, yet inexpensive incremental condition estimator to decide on which columns to be included. We describe an incremental condition estimator that can be used during a sparse QR factorization. We show that it is quite reliable and is well suited for use on parallel machines. We supply experimental results to support its effectiveness as well as suitability for parallel architectures.

A wavelet network (WN) is a feed-forward neural network that uses wavelets as activation functions for the neurons in its hidden layer. By predetermining the wavelet positions and dilations, the WN can turn into a linear regression model. The common approach for the construction of these WN families is to use least-squares type algorithms. In this letter, we propose a novel approach by formulating a WN as a sparse linear regression problem, which we call a sparse wavelet network (SWN). In this WN, the problem of calculating the unknown inner parameters of the network becomes that of finding the sparse solution of an under-determined system of linear equations. Our sparse solution algorithm is a non-convex sparse relaxation approach inspired by smoothed L0 (SL0), a distinguished sparse recovery algorithm. The proposed SWN can be applied as a tool for the prediction and identification of dynamical systems.

It has been observed that it is possible to represent the solution to certain interior scattering problems, such as simple cavities, in terms of solutions to the exterior problem and the pseudo inverse of the plane wave transform (PWT) operator. This may provide a pathway to sparse direct representations for interior problems provided that sparse representations of the exterior solution operator and the pseudo inverse of the PWT are available. This paper considers the latter possibility. The approach used consists of first compressing the PWT via spatially localized current modes that radiate far fields only to specific angular regions. It is shown that the resulting operator has a sparse QR factorization, which permits the efficient direct implementation of the desired pseudo inverse operation.

In QR factorization of a sparse m{times}n matrix A (m {ge} n) the orthogonal factor Q is often stored implicitly as a lower trapezoidal matrix H known as the Householder matrix. This paper presents a simple characterization of the row structure of Q, which could be used as the basis for a sparse data structure that can store Q explicitly. The new characterization is a simple extension of a well known row-oriented characterization of the structure of H. Hare, Johnson, Olesky, and van den Driessche have recently provided a complete sparsity analysis of the QR factorization. Let U be the matrix consisting of the first n columns of Q. Using results from, we show that the data structures for H and U resulting from our characterizations are tight when A is a strong Hall matrix. We also show that H and the lower trapezoidal part of U have the same sparsity characterization when A is strong Hall. We then show that this characterization can be extended to any weak Hall matrix that has been permuted into block upper triangular form. Finally, we show that permuting to block triangular form never increases the fill incurred during the factorization.

In this paper, we propose a new class of iteratively re-weighted least squares (IRLS) for sparse recovery problems. The proposed methods are inspired by constrained maximum-likelihood estimation under a Gaussian scale mixture (GSM) distribution assumption. In the noise-free setting, we provide sufficient conditions ensuring the convergence of the sequences generated by these algorithms to the set of fixed points of the maps that rule their dynamics and derive conditions verifiable a posteriori for the convergence to a sparse solution. We further prove that these algorithms are quadratically fast in a neighborhood of a sparse solution. We show through numerical experiments that the proposed methods outperform classical IRLS for l_p-minimization with p\in(0,1] in terms of speed and of sparsity-undersampling tradeoff and are robust even in presence of noise. The simplicity and the theoretical guarantees provided in this paper make this class of algorithms an attractive solution for sparse recovery problems

As with many other linear algebra algorithms, devising a portable implementation of sparse Cholesky factorization that performs well on the broad range of computer architectures currently available is a formidable challenge. Even after limiting our attention to machines with only one processor, as we have done in this report, there are still several interesting issues to consider. For dense matrices, it is well known that block factorization algorithms are the best means of achieving this goal. We take this approach for sparse factorization as well. This paper has two primary goals. First, we examine two sparse Cholesky factorization algorithms, the multifrontal method and a blocked left-looking sparse Cholesky method, in a systematic and consistent fashion, both to illustrate the strengths of the blocking techniques in general and to obtain a fair evaluation of the two approaches. Second, we assess the impact of various implementation techniques on time and storage efficiency, paying particularly close attention to the work-storage requirement of the two methods and their variants.

We present algorithms to determine the number of nonzeros in each row and column of the factors of a sparse matrix, for both the QR factorization and the LU factorization with partial pivoting. The algorithms use only the nonzero structure of the input matrix, and run in time nearly linear in the number of nonzeros in that matrix. They may be used to set up data structures or schedule parallel operations in advance of the numerical factorization. The row and column counts we compute are upper bounds on the actual counts. If the input matrix is strong Hall and there is no coincidental numerical cancellation, the counts are exact for QR factorization and are the tightest bounds possible for LU factorization. These algorithms are based on our earlier work on computing row and column counts for sparse Cholesky factorization, plus an efficient method to compute the column elimination tree of a sparse matrix without explicitly forming the product of the matrix and its transpose.

Several algorithms for orthogonal factorization on distributed memory multiprocessors are designed and implemented. Two of the algorithms employ Householder transformations, a third is based on Givens rotations, and a fourth hybrid algorithm uses Householder transformations and Givens rotations in different phases.The arithmetic and communication complexities of the algorithms are analyzed. The analyses show that the sequential arithmetic terms play a more important role than the communication terms in determining the running times and efficiencies of these algorithms. The hybrid algorithm is the fastest algorithm overall, since its arithmetic cost is lower than the Householder algorithms and its communication cost does not increase with the column length of the matrix. The observed execution times of the implementations on an iPSC-286 agree quite well with the complexity analyses. It is also shown that the efficiencies can be approximated using only the arithmetic costs of the algorithms.

Mathematical model construction of complicate physical phenomenon often leads to the setting and solving problems of parameters optimal control in differential equations in partial derivatives. Chosen equation with boundary and initial conditions is usually mathematical model basis of the object, which is under analysis. Optimal control of right-hand side function in non-linear problem for inhomogeneous biharmonic has been investigated. With the help of various gradient methods the problems of parameters control in such equations are solved successfully. Herewith linear problem is solved with the potential method on every step. In this work biharmonic potentials with logarithmic singularity are under consideration. That is why parameters optimization in these problems are conducted together with elimination of their incorrectness. It should be noted that system optimal control problems, described with biharmonic equations in irregular shape region, are studied poorly. Therefore the object of investigation is one of boundary conditions for inhomogeneous biharmonic equation. The circumstance which complicated the problem was irregularity of function domain. Clearly that this problem can be solved with computational mathematics methods. Without precise solution of linear problem, it is impossible to use gradient method, build converging iterative process, and obtain precise solution of optimization problem. Algorithm for linear boundary value problem solution with boundary integral equations overcomes this problem successfully. Physical examples of numerical implementation have been presented, analysis of obtained solutions have been conducted. Their accuracy, algorithm simplicity and time spent evidence about this approach promising for practical results obtaining in plate theory and mathematical physics problems successful numerical solving

BackgroundA genetic network can be represented as a directed graph in which a node corresponds to a gene and a directed edge specifies the direction of influence of one gene on another. The reconstruction of such networks from transcript profiling data remains an important yet challenging endeavor. A transcript profile specifies the abundances of many genes in a biological sample of interest. Prevailing strategies for learning the structure of a genetic network from high-dimensional transcript profiling data assume sparsity and linearity. Many methods consider relatively small directed graphs, inferring graphs with up to a few hundred nodes. This work examines large undirected graphs representations of genetic networks, graphs with many thousands of nodes where an undirected edge between two nodes does not indicate the direction of influence, and the problem of estimating the structure of such a sparse linear genetic network (SLGN) from transcript profiling data.ResultsThe structure learning task is cast as a sparse linear regression problem which is then posed as a LASSO (l1-constrained fitting) problem and solved finally by formulating a Linear Program (LP). A bound on the Generalization Error of this approach is given in terms of the Leave-One-Out Error. The accuracy and utility of LP-SLGNs is assessed quantitatively and qualitatively using simulated and real data. The Dialogue for Reverse Engineering Assessments and Methods (DREAM) initiative provides gold standard data sets and evaluation metrics that enable and facilitate the comparison of algorithms for deducing the structure of networks. The structures of LP-SLGNs estimated from the INSILICO1, INSILICO2 and INSILICO3 simulated DREAM2 data sets are comparable to those proposed by the first and/or second ranked teams in the DREAM2 competition. The structures of LP-SLGNs estimated from two published Saccharomyces cerevisae cell cycle transcript profiling data sets capture known regulatory associations. In each S. cerevisiae LP-SLGN, the number of nodes with a particular degree follows an approximate power law suggesting that its degree distributions is similar to that observed in real-world networks. Inspection of these LP-SLGNs suggests biological hypotheses amenable to experimental verification.ConclusionA statistically robust and computationally efficient LP-based method for estimating the topology of a large sparse undirected graph from high-dimensional data yields representations of genetic networks that are biologically plausible and useful abstractions of the structures of real genetic networks. Analysis of the statistical and topological properties of learned LP-SLGNs may have practical value; for example, genes with high random walk betweenness, a measure of the centrality of a node in a graph, are good candidates for intervention studies and hence integrated computational – experimental investigations designed to infer more realistic and sophisticated probabilistic directed graphical model representations of genetic networks. The LP-based solutions of the sparse linear regression problem described here may provide a method for learning the structure of transcription factor networks from transcript profiling and transcription factor binding motif data.

The problem of computing sparse solutions to linear inverse problems arises in a large number of signal processing application areas. We address the problem of finding sparse solutions to linear inverse problems when there are multiple measurement vectors (MMV) and the solutions are assumed to have a common, but unknown, sparsity profile. This is an important extension to the single measurement sparse solution problem that has been extensively studied in the past. Of particular interest are methods based on minimizing diversity measures. A measure appropriate for the multiple measurement problem is developed, and an algorithm is derived based on its minimization. The algorithm developed, M-FOCUSS, generalizes the focal underdetermined system solver (FOCUSS) algorithm developed for the single measurement case. The convergence of the algorithm is established and a simulation study is conducted to evaluate its effectiveness. The results clearly show the ability of M-FOCUSS to utilize multiple measurement vectors for accurate identification of the sparsity structure and sparse solution computation.

Feedback stabilization of a second-order system: A nonmodal approach

Geophysical Inverse Theory and Regularization Problems