A Schur complement approach to preconditioning sparse linear least-squares problems with some dense rows

The effectiveness of sparse matrix techniques for directly solving large-scale linear least-squares problems is severely limited if the system matrix A has one or more nearly dense rows. In this paper, we partition the rows of A into sparse rows and dense rows (A_s and A_d) and apply the Schur complement approach. A potential difficulty is that the reduced normal matrix AsTA_s is often rank-deficient, even if A is of full rank. To overcome this, we propose explicitly removing null columns of A_s and then employing a regularization parameter and using the resulting Cholesky factors as a preconditioner for an iterative solver applied to the symmetric indefinite reduced augmented system. We consider complete factorizations as well as incomplete Cholesky factorizations of the shifted reduced normal matrix. Numerical experiments are performed on a range of large least-squares problems arising from practical applications. These demonstrate the effectiveness of the proposed approach when combined with either a sparse parallel direct solver or a robust incomplete Cholesky factorization algorithm.     

On Andrunakievich’s chain and Koethe’s problem

In 1969 Andrunakievich asked whether one gets a ring without nonzero nil left ideals from an arbitrary ring R by factoring out the ideal A(R) which is the sum of all nil left ideals of R. Recently, it was shown that this problem is equivalent to Koethe’s problem. In this context one may consider the chain of ideals, which starts with A ₁(R) = A(R) ⊆ A ₂(R), where A ₂(R)/A ₁(R) = A(R/A ₁(R)), and extends by repeating this process. We study the properties of this chain and show that, assuming a negative solution of Koethe’s problem, this chain can terminate at any given ordinal number.     

Stability analysis of a general toeplitz systems solver

We show that a fast algorithm for theQR factorization of a Toeplitz or Hankel matrixA is weakly stable in the sense thatR ^T R is close toA ^T A. Thus, when the algorithm is used to solve the semi-normal equationsR ^TRx=A^T b, we obtain a weakly stable method for the solution of a nonsingular Toeplitz or Hankel linear systemAx=b. The algorithm also applies to the solution of the full-rank Toeplitz or Hankel least squares problem min ||Ax-b||₂.     

Matrix enlarging methods and their application

This paper explores several methods for matrix enlarging, where an enlarged matrixà is constructed from a given matrixA. The methods explored include matrix primitization, stretching and node splitting. Graph interpretations of these methods are provided. Solving linear problems using enlarged matrices yields the answer to the originalAx=b problem.à can exhibit several desirable properties. For example,à can be constructed so that the valence of any row and/or column is smaller than some desired number (≥4). This is beneficial for algorithms that depend on the square of the number of entries of a row or column. Most particularly, matrix enlarging can results in a reduction of the fill-in in theR matrix which occurs during orthogonal factorization as a result of dense rows. Numerical experiments support these conjectures.     

Solution of sparse quasi-square rectangular systems by Gaussian elimination

We present a general method for the linear least-squares solutionof overdetermined and underdetermined systems. The method isparticularly efficient when the coefficient matrix is quasi-square,that is when the number of rows and number of columns is almostthe same. The numerical methods for linear least-squares problemsand minimum-norm solutions do not generally take account ofthis special characteristic. The proposed method is based onLU factorization of the original quasi-square matrix A, assumingthat A has full rank. In the overdetermined case, the LU factorsare used to compute a basis for the null space of A^T. The right-handside vector b is then projected onto this subspace and the least-squaressolution is obtained from the solution of this reduced problem.In the case of underdetermined systems, the desired solutionis again obtained through the solution of a reduced system.The use of this method may lead to important savings in computationaltime for both dense and sparse matrices. It is also shown inthe paper that, even in cases where the matrices are quite small,sparse solvers perform better than dense solvers. Some practicalexamples that illustrate the use of the method are included.     

Splitting recursively enumerable subalgebras in recursive Boolean algebras

A Boolean algebraB= is recursive ifB is a recursive subset of ω and the operations Λ, v and ┌ are partial recursive. A subalgebraC ofB is recursive an (r.e.) ifC is a recursive (r.e.) subset of B. Given an r.e. subalgebraA, we sayA can be split into two r.e. subalgebrasA ₁ andA ₂ if (A ₁ ∪A ₂)^*=A andA ₁ ∩A ₂={0, 1}. In this paper we show that any nonrecursive r.e. subalgebra ofB can be split into two nonrecursive r.e. subalgebras ofB. This is a natural analogue of the Friedberg's splitting theorem in ω recursion theory.     

Commuting extensions and cubature formulae

Based on a novel point of view on 1-dimensional Gaussian quadrature, we present a new approach to d-dimensional cubature formulae. It is well known that the nodes of 1-dimensional Gaussian quadrature can be computed as eigenvalues of the so-called Jacobi matrix. The d-dimensional analog is that cubature nodes can be obtained from the eigenvalues of certain mutually commuting matrices. These are obtained by extending (adding rows and columns to) certain noncommuting matrices A₁,...,A_d, related to the coordinate operators x₁,...,x_d, in R^d. We prove a correspondence between cubature formulae and “commuting extensions” of A₁,...,A_d, satisfying a compatibility condition which, in appropriate coordinates, constrains certain blocks in the extended matrices to be zero. Thus, the problem of finding cubature formulae can be transformed to the problem of computing (and then simultaneously diagonalizing) commuting extensions. We give a general discussion of existence and of the expected size of commuting extensions and briefly describe our attempts at computing them.     

Efficient computation of sparse approximate inverses

We investigate different methods for computing a sparse approximate inverse M for a given sparse matrix A by minimizing ∥AM − E∥ in the Frobenius norm. Such methods are very useful for deriving preconditioners in iterative solvers, especially in a parallel environment. We compare different strategies for choosing the sparsity structure of M and different ways for solving the small least squares problem that are related to the computation of each column of M. Especially we show how we can take full advantage of the sparsity of A. © 1998 John Wiley & Sons, Ltd.     

Some mathematical problems in seriation

Some mathematical aspects of seriation are studied in this paper. Certain conditions on an abundance or an incidence matrix have been given in the past which imply that there exists a permutation of its rows so that the resulting matrix is a Q matrix (in which case the original matrix is said to be a pre-Q). These types of results have applications to chronologically ordering archaeological provenances under certain circumstances. Unfortunately these conditions are deficient both theoretically and practically, in that for much archaeological data the conditions are not necessarily true yet the corresponding provenances do have chronological orderings. Here we are able to generalize these results in two ways. First we are able to establish necessary and sufficient conditions on the rows of a matrix for it to be pre-Q. These conditions are local in that they concern only certain triples and quadruples of the rows. Secondly, we are able to interpret seriation in terms of a ternary relation R on a set A and prove the results in this general context. In this form the theorem says that if only certain of the triples and quadruples are R-strings, then the whole set A is an R-string, and so has a linear order consistent with the ternary relation R. This would appear to generalize a theorem of P. C. Fishburn. Both aspects of the generalization mean that the results stated herein have a wider applicability than those given heretofore. Possibly more importantly than this is that they lead to numerical invariants, called the fixing number and the related linear rigidity, of such an R-string on A. The archaeological interpretation of these is given in the paper and data supplied which illustrates this point. Finally various other conditions on products and representations of relations are stated which imply that A is an R-string. One of these generalizes and completes a theorem of D. G. Kendall.     

Solution of large‐scale weighted least‐squares problems

A sequence of least‐squares problems of the form min_y∥G^1/2(A^T y?h)∥₂, where G is an n×n positive‐definite diagonal weight matrix, and A an m×n (m?n) sparse matrix with some dense columns; has many applications in linear programming, electrical networks, elliptic boundary value problems, and structural analysis. We suggest low‐rank correction preconditioners for such problems, and a mixed solver (a combination of a direct solver and an iterative solver). The numerical results show that our technique for selecting the low‐rank correction matrix is very effective. Copyright © 2002 John Wiley & Sons, Ltd.     

Least squares solutions to AX = B for bisymmetric matrices under a central principal submatrix constraint and the optimal approximation

A matrix A∈R^n×n is called a bisymmetric matrix if its elements a_i,j satisfy the properties a_i,j=a_j,i and a_i,j=a_n-j+1,n-i+1 for 1?i,j?n. This paper considers least squares solutions to the matrix equation AX=B for A under a central principal submatrix constraint and the optimal approximation. A central principal submatrix is a submatrix obtained by deleting the same number of rows and columns in edges of a given matrix. We first discuss the specified structure of bisymmetric matrices and their central principal submatrices. Then we give some necessary and sufficient conditions for the solvability of the least squares problem, and derive the general representation of the solutions. Moreover, we also obtain the expression of the solution to the corresponding optimal approximation problem.     

Stable signal recovery from incomplete and inaccurate measurements

Suppose we wish to recover a vector x₀ ∈ ?^?? (e.g., a digital signal or image) from incomplete and contaminated observations y = A x₀ + e; A is an ?? × ?? matrix with far fewer rows than columns (?? ? ??) and e is an error term. Is it possible to recover x₀ accurately based on the data y? To recover x₀, we consider the solution x^# to the ??₁‐regularization problem where ? is the size of the error term e. We show that if A obeys a uniform uncertainty principle (with unit‐normed columns) and if the vector x₀ is sufficiently sparse, then the solution is within the noise level As a first example, suppose that A is a Gaussian random matrix; then stable recovery occurs for almost all such A's provided that the number of nonzeros of x₀ is of about the same order as the number of observations. As a second instance, suppose one observes few Fourier samples of x₀; then stable recovery occurs for almost any set of ?? coefficients provided that the number of nonzeros is of the order of ??/(log ??)⁶. In the case where the error term vanishes, the recovery is of course exact, and this work actually provides novel insights into the exact recovery phenomenon discussed in earlier papers. The methodology also explains why one can also very nearly recover approximately sparse signals. © 2006 Wiley Periodicals, Inc.     

Finite Left-Distributive Algebras and Embedding Algebras

We consider algebras with one binary operation · and one generator (monogenic) and satisfying the left distributive lawa·(b·c)=(a·b)·(a·c). One can define a sequence of finite left-distributive algebrasA_n, and then take a limit to get an infinite monogenic left-distributive algebraA_∞. Results of Laver and Steel assuming a strong large cardinal axiom imply thatA_∞is free; it is open whether the freeness ofA_∞can be proved without the large cardinal assumption, or even in Peano arithmetic. The main result of this paper is the equivalence of this problem with the existence of a certain algebra of increasing functions on natural numbers, called anembedding algebra. Using this and results of the first author, we conclude that the freeness ofA_∞is unprovable in primitive recursive arithmetic.     

Component-wise perturbation analysis and error bounds for linear least squares solutions

Perturbation bounds for the linear least squares problem min _x Ax –b₂ corresponding tocomponent-wise perturbations in the data are derived. These bounds can be computed using a method of Hager and are often much better than the bounds derived from the standard perturbation analysis. In particular this is true for problems where the rows ofA are of widely different magnitudes. Generalizing a result by Oettli and Prager, we can use the bounds to compute a posteriori error bounds for computed least squares solutions.     

Sparse direct factorizations through unassembled hyper-matrices

We set out to efficiently compute the solution of a sequence of linear systems A_ix_i = b_i, where the matrix A_i is tightly related to matrix A_{i –1}. In the setting of an hp -adaptive Finite Element Method, the sequence of matrices A_i results from successive local refinements of the problem domain. At any step i > 1, a factorization already exists and it is the updated linear system relative to the refined mesh for which a factorization must be computed in the least amount of time. This observation holds the promise of a tremendous reduction in the cost of an individual refinement step. We argue that traditional matrix storage schemes, whether dense or sparse, are a bottleneck, limiting the potential efficiency of the solvers. We propose a new hierarchical data structure, the Unassembled Hyper-Matrix (UHM), which allows the matrix to be stored as a tree of unassembled element matrices, hierarchically ordered to mirror the refinement history of the physical domain. The factorization of such an UHM proceeds in terms of element matrices, only assembling nodes when they need to be eliminated. Efficiency comes in terms of both performance and space requirements. (© 2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim)     

具二次等式约束最小二乘问题的一种极端情形

本文研究具有二次等式约束的最小二乘问题(LSS): min‖Ax-b‖₂ s.t. ‖x‖₂=1, 其中A∈R^m×n, b∈R^m, 并假定‖A⁺b‖₂<1.重点关注一个极端情形: ‖A⁺b‖₂≈0. 敏度分析表明,这是一种病态问题. 基于Padé逼近给出了一种迭代解法. 数值算例表明,新方法在速度上较已有方法有优势.     

Gauss-seidel method for least-distance problems

In this paper, we reformulate the least-distance problems with bounded inequality constraints as an unconstrained convex minimization problem, which is equivalent to a system of piecewise linear equationsA(a+A ^T y) _c ^d =b. The proposed Gauss-Seidel method for solving the problems is easy to implement and behaves very well when the number of rows ofA is much less than the number of columns ofA. Moreover, we prove that the Gauss-Seidel method has a linear convergence rate.The authors would like to thank Jinshui Qin who did some preliminary numerical testing for their research.     

On co-positive,semi-monotoneQ-matrices

In this paper we consider not necessarily symmetric co-positive as well as semi-monotoneQ-matrices and give a set of sufficient conditions for such matrices to beR ₀-matrices. We give several examples to show the sharpness of our results. Construction of these examples is based on the following elementary proposition: IfA is a square matrix of ordern whose first two rows are identical and bothA ₁₁ andA ₂₂ areQ-matrices whereA _ii stands for the principal submatrix ofA obtained by deleting rowi and columni fromA, thenA is aQ-matrix.Dedicated to our colleague and friend B. Ramachandran on his sixtieth birthday.Corresponding author.     

Probabilistic Analysis of Condition Numbers for Linear Programming

In this paper, we provide bounds for the expected value of the log of the condition number C(A) of a linear feasibility problem given by a n × m matrix A (Ref. 1). We show that this expected value is O(min{n, m log n}) if n > m and is O(log n) otherwise. A similar bound applies for the log of the condition number C _R(A) introduced by Renegar (Ref. 2).     

A Faster and Simpler Recursive Algorithm for the LAPACK Routine DGELS

We present new algorithms for computing the linear least squares solution to overdetermined linear systems and the minimum norm solution to underdetermined linear systems. For both problems, we consider the standard formulation min AX – B _F and the transposed formulation min A ^T X – B _F , i.e, four different problems in all. The functionality of our implementation corresponds to that of the LAPACK routine DGELS. The new implementation is significantly faster and simpler. It outperforms the LAPACK DGELS for all matrix sizes tested. The improvement is usually 50–100% and it is as high as 400%. The four different problems of DGELS are essentially reduced to two, by use of explicit transposition of A. By explicit transposition we avoid computing Householder transformations on vectors with large stride. The QR factorization of block columns of A is performed using a recursive level-3 algorithm. By interleaving updates of B with the factorization of A, we reduce the number of floating point operations performed for the linear least squares problem. By avoiding redundant computations in the update of B we reduce the work needed to compute the minimum norm solution. Finally, we outline fully recursive algorithms for the four problems of DGELS as well as for QR factorization.This revised version was published online in October 2005 with corrections to the Cover Date.