Showing papers on "Bidiagonalization published in 2018"

Efficient generalized Golub–Kahan based methods for dynamic inverse problems

Abstract: We consider efficient methods for computing solutions to and estimating uncertainties in dynamic inverse problems, where the parameters of interest may change during the measurement procedure. Compared to static inverse problems, incorporating prior information in both space and time in a Bayesian framework can become computationally intensive, in part, due to the large number of unknown parameters. In these problems, explicit computation of the square root and/or inverse of the prior covariance matrix is not possible. In this work, we develop efficient, iterative, matrix-free methods based on the generalized Golub-Kahan bidiagonalization that allow automatic regularization parameter and variance estimation. We demonstrate that these methods can be more flexible than standard methods and develop efficient implementations that can exploit structure in the prior, as well as possible structure in the forward model. Numerical examples from photoacoustic tomography, deblurring, and passive seismic tomography demonstrate the range of applicability and effectiveness of the described approaches. Specifically, in passive seismic tomography, we demonstrate our approach on both synthetic and real data. To demonstrate the scalability of our algorithm, we solve a dynamic inverse problem with approximately $43,000$ measurements and $7.8$ million unknowns in under $40$ seconds on a standard desktop.

Solution methods for linear discrete ill-posed problems for color image restoration

Abstract: This work discusses four algorithms for the solution of linear discrete ill-posed problems with several right-hand side vectors. These algorithms can be applied, for instance, to multi-channel image restoration when the image degradation model is described by a linear system of equations with multiple right-hand sides that are contaminated by errors. Two of the algorithms are block generalizations of the standard Golub–Kahan bidiagonalization method with the block size equal to the number of channels. One algorithm uses standard Golub–Kahan bidiagonalization without restarts for all right-hand sides. These schemes are compared to standard Golub–Kahan bidiagonalization applied to each right-hand side independently. Tikhonov regularization is used to avoid severe error propagation. Numerical examples illustrate the performance of these algorithms. Applications include the restoration of color images.

An iterative generalized Golub-Kahan algorithm for problems in structural mechanics

Abstract: This paper studies the Craig variant of the Golub-Kahan bidiagonalization algorithm as an iterative solver for linear systems with saddle point structure. Such symmetric indefinite systems in 2x2 block form arise in many applications, but standard iterative solvers are often found to perform poorly on them and robust preconditioners may not be available. Specifically, such systems arise in structural mechanics, when a semidefinite finite element stiffness matrix is augmented with linear multi-point constraints via Lagrange multipliers. Engineers often use such multi-point constraints to introduce boundary or coupling conditions into complex finite element models. The article will present a systematic convergence study of the Golub-Kahan algorithm for a sequence of test problems of increasing complexity, including concrete structures enforced with pretension cables and the coupled finite element model of a reactor containment building. When the systems are suitably transformed using augmented Lagrangians on the semidefinite block and when the constraint equations are properly scaled, the Golub-Kahan algorithm is found to exhibit excellent convergence that depends only weakly on the size of the model. The new algorithm is found to be robust in practical cases that are otherwise considered to be difficult for iterative solvers.

The Low Rank Approximations and Ritz Values in LSQR For Linear Discrete Ill-Posed Problems

Abstract: LSQR and its mathematically equivalent CGLS have been popularly used over the decades for large-scale linear discrete ill-posed problems, where the iteration number $k$ plays the role of the regularization parameter. It has been long known that if the Ritz values in LSQR converge to the large singular values of $A$ in natural order until its semi-convergence then LSQR must have the same the regularization ability as the truncated singular value decomposition (TSVD) method and can compute a 2-norm filtering best possible regularized solution. However, hitherto there has been no definitive rigorous result on the approximation behavior of the Ritz values in the context of ill-posed problems. In this paper, for severely, moderately and mildly ill-posed problems, we give accurate solutions of the two closely related fundamental and highly challenging problems on the regularization of LSQR: (i) How accurate are the low rank approximations generated by Lanczos bidiagonalization? (ii) Whether or not the Ritz values involved in LSQR approximate the large singular values of $A$ in natural order? We also show how to judge the accuracy of low rank approximations reliably during computation without extra cost. Numerical experiments confirm our results.

SPMR: A Family of Saddle-Point Minimum Residual Solvers

Abstract: We introduce a new family of saddle-point minimum residual methods for iteratively solving saddle-point systems using a minimum or quasi-minimum residual approach. No symmetry assumptions are made. The basic mechanism underlying the method is a novel simultaneous bidiagonalization procedure that yields a simplified saddle-point matrix on a projected Krylov-like subspace and allows for a monotonic short-recurrence iterative scheme. We develop a few variants, demonstrate the advantages of our approach, derive optimality conditions, and discuss connections to existing methods. Numerical experiments illustrate the merits of this new family of methods.

Weighted Golub-Kahan-Lanczos bidiagonalization algorithms

Abstract: We present weighted Golub-Kahan-Lanczos algorithms. We demonstrate their applications to the eigenvalue problem of a product of two symmetric positive definite matrices and an eigenvalue problem for the linear response problem. A convergence analysis is provided and numerical test results are reported. As another application we make a connection between the proposed algorithms and the preconditioned conjugate gradient (PCG) method.

LNLQ: an iterative method for least-norm problems with an error minimization property

Abstract: We describe LNLQ for solving the least-norm problem min $\|x\|$ subject to $Ax=b$, using the Golub--Kahan bidiagonalization of $[b\ A]$. Craig's method is known to be equivalent to applying the con...

Quantifying uncertainties in Large-scale Bayesian linear inverse problems using Krylov subspace methods

Abstract: For linear inverse problems with a large number of unknown parameters, uncertainty quantification remains a challenging task. In this work, we use Krylov subspace methods to approximate the posterior covariance matrix and describe efficient methods for exploring the posterior distribution. Assuming that Krylov methods (e.g., based on the generalized Golub-Kahan bidiagonalization) have been used to compute an estimate of the solution, we get an approximation of the posterior covariance matrix for `free.' We provide theoretical results that quantify the accuracy of the approximation and of the resulting posterior distribution. Then, we describe efficient methods that use the approximation to compute measures of uncertainty, including the Kullback-Liebler divergence. We present two methods that use preconditioned Lanczos methods to efficiently generate samples from the posterior distribution. Numerical examples from tomography demonstrate the effectiveness of the described approaches.

LSMB: Minimizing the Backward Error for Least-Squares Problems

Abstract: An iterative method, LSMB, is given for solving $\min_x \|Ax-b\|_2$. LSMB is based on the Golub--Kahan bidiagonalization process and is constructed so that an objective function closely related to ...

Uncertainty quantification in large Bayesian linear inverse problems using Krylov subspace methods

Abstract: For linear inverse problems with a large number of unknown parameters, uncertainty quantification remains a challenging task. In this work, we use Krylov subspace methods to approximate the posterior covariance matrix and describe efficient methods for exploring the posterior distribution. Assuming that Krylov methods (e.g., based on the generalized Golub-Kahan bidiagonalization) have been used to compute an estimate of the solution, we get an approximation of the posterior covariance matrix for `free.' We provide theoretical results that quantify the accuracy of the approximation and of the resulting posterior distribution. Then, we describe efficient methods that use the approximation to compute measures of uncertainty, including the Kullback-Liebler divergence. We present two methods that use preconditioned Lanczos methods to efficiently generate samples from the posterior distribution. Numerical examples from tomography demonstrate the effectiveness of the described approaches.

A new implementation of LSMR algorithm for the quaternionic least squares problem

Abstract: This paper is endeavored to present a new version of the LSMR algorithm for solving the linear least squares problem in quaternion field, by means of direct quaternion arithmetics rather than the usually used real or complex representation methods. The present new algorithm is based on the classical Golub-Kahan bidiagonalization process, but is instead of using two QR factorizations. It has several advantages as follows: (i) does not make the scale of the problem dilate exponentially, compared to the conventional complex representation or real representation methods, (ii) has monotonic and smooth convergence behavior, compared to the Q-LSQR algorithm, and (iii) the new algorithm is more straightforward, and there is no expensive matrix inversion or decomposition. It may reduce the number of iterations in some cases. The performances of the algorithm are illustrated by some numerical experiments.

Fast computation of the principal components of genotype matrices in Julia

Abstract: Finding the largest few principal components of a matrix of genetic data is a common task in genome-wide association studies (GWASs), both for dimensionality reduction and for identifying unwanted factors of variation. We describe a simple random matrix model for matrices that arise in GWASs, showing that the singular values have a bulk behavior that obeys a Marchenko-Pastur distributed with a handful of large outliers. We also implement Golub-Kahan-Lanczos (GKL) bidiagonalization in the Julia programming language, providing thick restarting and a choice between full and partial reorthogonalization strategies to control numerical roundoff. Our implementation of GKL bidiagonalization is up to 36 times faster than software tools used commonly in genomics data analysis for computing principal components, such as EIGENSOFT and FlashPCA, which use dense LAPACK routines and randomized subspace iteration respectively.

Low-complexity geometric mean value decomposition precoding implementation method based on bidiagonalization

Abstract: The present invention provides a low-complexity geometric mean value decomposition precoding implementation method based on bidiagonalization. The method comprises the following steps: (1) calculatinga conjugate transpose of a channel matrix and own product; (2) based on a given Hermitian matrix bidiagonalization method, changing the channel matrix into a bidiagonal matrix through Givens rotation; (3) based on a given geometric mean value decomposition method, changing the bidiagonal matrix into an upper triangular matrix where the diagonal elements all equal to the geometric mean value of the channel matrix eigenvalues through Givens rotation; (4) constructing a precoding matrix with the geometric mean value decomposed, that is, the product of all Givens right rotation matrices. The technical scheme can determine the number of iterations, and uses the Hermitian matrix for solutions further to reduce the implementation complexity and reduce the use of the CORDIC (Coordinate Rotation Digital Computer) module. The method utilizes the properties of the Hermitian matrix to effectively reduce the implementation complexity of geometric mean value decomposition precoding based on bidiagonalization.

Large-scale Inversion of Magnetic Data Using Golub-Kahan Bidiagonalization with Truncated Generalized Cross Validation for Regularization Parameter Estimation

Abstract: In this paper a fast method for large-scale sparse inversion of magnetic data is considered. The L1-norm stabilizer is used to generate models with sharp and distinct interfaces. To deal with the non-linearity introduced by the L1-norm, a model-space iteratively reweighted least squares algorithm is used. The original model matrix is factorized using the Golub-Kahan bidiagonalization that projects the problem onto a Krylov subspace with a significantly reduced dimension. The model matrix of the projected system inherits the ill-conditioning of the original matrix, but the spectrum of the projected system accurately captures only a portion of the full spectrum. Equipped with the singular value decomposition of the projected system matrix, the solution of the projected problem is expressed using a filtered singular value expansion. This expansion depends on a regularization parameter which is determined using the method of Generalized Cross Validation (GCV), but here it is used for the truncated spectrum. This new technique, Truncated GCV (TGCV), is more effective compared with the standard GCV method. Numerical results using a synthetic example and real data demonstrate the efficiency of the presented algorithm.

A staggered eigensolver based on sparse matrix bidiagonalization

Abstract: We present a method for calculating eigenvectors of the staggered Dirac operator based on the Golub-Kahan-Lanczos bidiagonalization algorithm. Instead of using orthogonalization during the bidiagonalization procedure to increase stability, we choose to stabilize the method by combining it with an outer iteration that refines the approximate eigenvectors obtained from the inner bidiagonalization procedure. We discuss the performance of the current implementation using QEX and compare with other methods.

Fast exponential filtering regularization photoacoustic imaging reconstruction method based on Lanczos bidiagonalization

Abstract: The invention discloses a fast exponential filtering regularization photoacoustic imaging reconstruction method based on Lanczos bidiagonalization, and belongs to the field of medical image processing. When short-pulse laser irradiates to-be-tested biological tissue, the tissue absorbs the heat generated by the energy of the laser, heat conduction is neglected and an approximate heat equation of photoacoustic imaging is constructed; based on an open-source tool box of the MATLAB, a photoacoustic equation is solved. A photoacoustic signal wave at the edge of an imaging region is collected through a set of energy converters. The process for collecting these signals is expressed as a time-varying causal system. A k-Wave tool box is used for constructing a system matrix of the system, and impulse response (IR) is recorded by pixels one by one for a complete imaging domain. The method can ensure the imaging quality, meanwhile the imaging speed is improved, the photoacoustic signals can be accurately reconstructed, and meanwhile the computational efficiency is high.

Singular Value Decomposition Approximation via Kronecker Summations for Imaging Applications

Abstract: In this paper we propose an approach to approximate a truncated singular value decomposition of a large structured matrix. By first decomposing the matrix into a sum of Kronecker products, our approach can be used to approximate a large number of singular values and vectors more efficiently than other well known schemes, such as randomized matrix algorithms or iterative algorithms based on Golub-Kahan bidiagonalization. We provide theoretical results and numerical experiments to demonstrate the accuracy of our approximation and show how the approximation can be used to solve large scale ill-posed inverse problems, either as an approximate filtering method, or as a preconditioner to accelerate iterative algorithms.