A fast semi-direct least squares algorithm for hierarchically block separable matrices

We present a fast algorithm for linear least squares problems governed by hierarchically block separable (HBS) matrices. Such matrices are generally dense but data-sparse and can describe many important operators including those derived from asymptotically smooth radial kernels that are not too oscillatory. The algorithm is based on a recursive skeletonization procedure that exposes this sparsity and solves the dense least squares problem as a larger, equality-constrained, sparse one. It relies on a sparse QR factorization coupled with iterative weighted least squares methods. In essence, our scheme consists of a direct component, comprised of matrix compression and factorization, followed by an iterative component to enforce certain equality constraints. At most two iterations are typically required for problems that are not too ill-conditioned. For an $M \times N$ HBS matrix with $M \geq N$ having bounded off-diagonal block rank, the algorithm has optimal $\mathcal{O} (M + N)$ complexity. If the rank increases with the spatial dimension as is common for operators that are singular at the origin, then this becomes $\mathcal{O} (M + N)$ in 1D, $\mathcal{O} (M + N^{3/2})$ in 2D, and $\mathcal{O} (M + N^{2})$ in 3D. We illustrate the performance of the method on both over- and underdetermined systems in a variety of settings, with an emphasis on radial basis function approximation and efficient updating and downdating.Comment: 24 pages, 8 figures, 6 tables; to appear in SIAM J. Matrix Anal. App

Using parallel computation to apply the singular value decomposition (SVD) in solving for large Earth gravity fields based on satellite data

textUsing satellite data only to estimate for an Earth gravity field introduces the problem of an ill-conditioned system of equations. This mathematical difficulty amplifies as the number of unknown gravity field parameters increases, requiring a stabilization of the inversion for solution. But the number of parameters to be estimated can also be too large to allow inversion using a sequential algorithm (one computer processor). Therefore the challenge is two-fold. A stabilized inversion must be performed with a parallel (multi-processor) algorithm. Thus, new code was developed in the parallel computing infrastructure of Parallel Linear Algebra Package (PLAPACK) to achieve the task of applying the Singular Value Decomposition (SVD) to invert for (and stabilize) very large gravity fields of well over 25,000 unknown parameters. This new code is given the name (Parallel LArge Svd Solver) PLASS. The choice of the SVD was made because it offers multiple opportunities of stabilization techniques. Poorly observed parameter corrections are removed from the culpable eigenspace of the normal matrix of CHAMP or the singular vector space of the upper R triangular matrix of GRACE. Solutions were stabilized based on the removal of either eigenvalues or singular values using four different standard optimization criteria: Inspection, Relative Error, Norm Norm minimization, trace of the Mean Square Error (MSE) matrix, and with a fifth method, independently introduced for this investigation, that optimizes removal of eigenvalues or singular values based on Kaula’s power rule of thumb. This method is given the name “Kaula Eigenvalue (KEV) or Kaula Singular Value (KSV) relation”. For the gravity fields of this investigation, orbital fits, geodetic evaluations and error propagations of the best of the resulting SVD gravity fields were performed, and shown to be comparable to the CHAMP solution obtained by the GeoForschungsZentrum (GFZ) and to the full rank GRACE solution obtained by the Center for Space Research (CSR).Aerospace Engineering and Engineering Mechanic

A Householder-based algorithm for Hessenberg-triangular reduction

The QZ algorithm for computing eigenvalues and eigenvectors of a matrix pencil $A - \lambda B$ requires that the matrices first be reduced to Hessenberg-triangular (HT) form. The current method of choice for HT reduction relies entirely on Givens rotations regrouped and accumulated into small dense matrices which are subsequently applied using matrix multiplication routines. A non-vanishing fraction of the total flop count must nevertheless still be performed as sequences of overlapping Givens rotations alternately applied from the left and from the right. The many data dependencies associated with this computational pattern leads to inefficient use of the processor and poor scalability. In this paper, we therefore introduce a fundamentally different approach that relies entirely on (large) Householder reflectors partially accumulated into block reflectors, by using (compact) WY representations. Even though the new algorithm requires more floating point operations than the state of the art algorithm, extensive experiments on both real and synthetic data indicate that it is still competitive, even in a sequential setting. The new algorithm is conjectured to have better parallel scalability, an idea which is partially supported by early small-scale experiments using multi-threaded BLAS. The design and evaluation of a parallel formulation is future work

GMRES implementations and residual smoothing techniques for solving ill-posed linear systems

AbstractThere are verities of useful Krylov subspace methods to solve nonsymmetric linear system of equations. GMRES is one of the best Krylov solvers with several different variants to solve large sparse linear systems. Any GMRES implementation has some advantages. As the solution of ill-posed problems are important. In this paper, some GMRES variants are discussed and applied to solve these kinds of problems. Residual smoothing techniques are efficient ways to accelerate the convergence speed of some iterative methods like CG variants. At the end of this paper, some residual smoothing techniques are applied for different GMRES methods to test the influence of these techniques on GMRES implementations