Search

Purpose: Developing a computationally efficient automated method for the optimal choice of regularization parameter in diffuse optical tomography. Methods: The least-squares QR (LSQR)-type method that uses Lanczos bidiagonalization is known to be computationally efficient in performing the reconstruction procedure in diffuse optical tomography. The same is effectively deployed via an optimization procedure that uses the simplex method to find the optimal regularization parameter. The proposed LSQR-type method is compared with the traditional methods such as L-curve, generalized cross-validation (GCV), and recently proposed minimal residual method (MRM)-based choice of regularization parameter using numerical and experimental phantom data. Results: The results indicate that the proposed LSQR-type and MRM-based methods performance in terms of reconstructed image quality is similar and superior compared to L-curve and GCV-based methods. The proposed method computational complexity is at least five times lower compared to MRM-based method, making it an optimal technique. Conclusions: The LSQR-type method was able to overcome the inherent limitation of computationally expensive nature of MRM-based automated way finding the optimal regularization parameter in diffuse optical tomographic imaging, making this method more suitable to be deployed in real-time. (C) 2013 American Association of Physicists in Medicine. http://dx.doi.org/10.1118/1.4792459]

Publisher Summary In linear regularization problems of the form min value the quotient singular value decomposition (QSVD) of the matrix pair (A, B) plays an essential role as analysis tool. In particular, the largest quotient singular values and the associated singular vectors are important. This chapter describes an algorithm based on Lanczos bidiagonalization that approximates these quantities. Regularization is needed in order to compute a stabilized solution which is less sensitive to the errors. The chapter discusses the most popular—Tikhnonov/Phillips regularization method.

Issues raised in the implementation of dense linear algebra algorithms on a distributed memory architecture are discussed along with considerations that arise when the algorithms are implemented as general purpose library routines. We focus on a few methods which have proven to be important in the CMSSL implementation of elimination based algorithms such as LU and QR decomposition and bidiagonalization.

We present a novel parameter choice strategy for the conjugate gradient regularization algorithm which does not assume a priori information about the magnitude of the measurement error. Our approach imitates the truncated singular value decomposition within the Krylov subspaces associated with the normal equations. Conjugate gradient is implemented using the Lanczos bidiagonalization process with reorthogonalization. We compare our method with one proposed by Hanke and Raus and illustrate its performance with numerical experiments, including an inverse problem of acoustic source detection.

In this paper, we consider ill-posed problems which discretize to linear least squares problems with matrices K of high dimensions. The algorithm proposed uses K only as an operator and does not need to explicitly store or modify it. A method related to one of Lanczos is used to project the problem onto a subspace for which K is bidiagonal. It is then an easy matter to solve the projected problem by standard regularization techniques. These ideas are illustrated with some integral equations of the first kind with convolution kernels, and sample numerical results are given.

The Lanczos method with the Golub and Kahan bidiagonalization algorithm is used to solve non‐symmetric diagonal dominant simultaneous equations. The method is very suitable for sparse matrix and vector computer. The problem with loss of orthogonal property is dealt with by restating the iteration. Numerical examples together with FORTRAN 77 routine are given to illustrate the algorithm. Both scalar and vector CPU times are given for comparison.

Techniques to diagonalize and to bidiagonalize a matrix are discussed.

An iterative method is given for solving Ax ~ffi b and minU Ax b 112, where the matrix A is large and sparse. The method is based on the bidiagonalization procedure of Golub and Kahan. It is analytically equivalent to the standard method of conjugate gradients, but possesses more favorable numerical properties. Reliable stopping criteria are derived, along with estimates of standard errors for x and the condition number of A. These are used in the FORTRAN implementation of the method, subroutine LSQR. Numerical tests are described comparing I~QR with several other conjugate-gradient algorithms, indicating that I~QR is the most reliable algorithm when A is ill-conditioned.

Two algorithms derived by Paige (1974), to solve the linear least squares problems based on Golub and Kahan's bidiagonalization algorithm to reduce a general matrix to a bidiagonal form. A preconditioning of the bidiagonalization algorithm is shown to be very effective to improve the rate of convergence of the bidiagonalization algorithm with a general application to solve a nonlinear system with non-symmetric Jacobian by linearization. Two preconditioning algorithms have been proposed, a comparison taken place to prove the efficiency of these algorithms.

Manne, R., 1987. Analysis of two partial-least-squares algorithms for multivariate calibration. Chemometrics and Intelligent Laboratory Systems , 2: 187-197. Two algorithms for multivariate calibration are analysed in terms of standard linear regression theory. The matrix inversion problem of linear regression is shown to be solved by transformations to a bidiagonal form in PLS1 and to a triangular form in PLS2. PLS1 gives results identical with the bidiagonalization algorithm by Golub and Kahan, similar to the method of conjugate gradients. The general efficiency of the algorithms is discussed.

An algorithm given by Golub and Kahan [2] for reducing a general matrix to bidiagonal form is shown to be very important for large sparse matrices. The singular values of the matrix are those of the bidiagonal form, and these can be easily computed. The bidiagonalization algorithm is shown to be the basis of important methods for solving the linear least squares problem for large sparse matrices. Eigenvalues of certain 2-cyclic matrices can also be efficiently computed using this bidiagonalization.

An iterative method based on Lanczos bidiagonalization is developed for computing regularized solutions of large and sparse linear systems, which arise from discretizations of ill-posed problems in partial differential or integral equations. Determination of the regularization parameter and termination criteria are discussed. Comments are given on the computational implementation of the algorithm.

A numerically stable and fairly fast scheme is described to compute the unitary matrices U and V which transform a given matrix A into a diagonal form $\Sigma = U^ * AV$, thus exhibiting A’s singular values on $\Sigma $’s diagonal. The scheme first transforms A to a bidiagonal matrix J, then diagonalizes J. The scheme described here is complicated but does not suffer from the computational difficulties which occasionally afflict some previously known methods. Some applications are mentioned, in particular the use of the pseudo-inverse $A^I = V\Sigma ^I U^* $ to solve least squares problems in a way which dampens spurious oscillation and cancellation.

This paper presents a fast algorithm for bidiagonalizing a Hankel matrix. An m × n Hankel matrix is reduced to a real bidiagonal matrix in O ( ( m + n ) n log ( m + n ) ) floating-point operations (flops) using the Lanczos method with modified partial orthogonalization and reset schemes to improve its stability. Performance improvement is achieved by exploiting the Hankel structure, as fast Hankel matrix–vector multiplication is used. The accuracy and efficiency of the algorithm are demonstrated by our numerical experiments.

In this paper, an improved algorithm PSVD for computing the singular subspace of a matrix corresponding to its smallest singular values is presented. As only a basis of the desired singular subspace is needed, the classical Singular Value Decomposition (SVD) algorithm can be modified in three ways. First, the Householder transformations of the bidiagonalization need only to be applied on the base vectors of the desired singular subspace. Second, the bidiagonal must only be partially diagonalized and third, the convergence rate of the iterative diagonalization can be improved by an appropriate choice between QR and QL iteration steps. An analysis of the operation counts, as well as computational results, show the relative efficiency of PSVD with respect to the classical SVD algorithm. Depending on the gap, the desired numerical accuracy and the dimension of the desired subspace, PSVD can be three times faster than the classical SVD algorithm while the same accuracy can be maintained. The new algorithm can be successfully used in total least squares applications, in the computation of the null space of a matrix and in solving (non) homogeneous linear equations. Based on PSVD a very efficient and reliable algorithm is also derived for solving nonhomogeneous equations.

This paper presents a comparison of certain direct algorithms for computing the Moore-Penrose inverse, for matrices of full column rank, from the point of view of numerical stability. It is proved that the algorithm using Householder QR decomposition, implemented in floating point arithmetic, is forward stable but only conditionally mixed forward-backward stable. A similar result holds also for the Classical Gram-Schmidt algorithm with reorthogonalization (CGS2). This algorithm was developed and analyzed by Abdelmalek (BIT, 11(4):354–367, 1971) and its detailed error analysis was given in Giraud et al. (Numer. Math. 101(1):87–100, 2005). We incorporate the recent result of Byers and Xu (SIAM J. Matrix Anal. Appl. 30:822–843, 2008) to show that the algorithm based on the Golub-Kahan bidiagonalization is mixed forward-backward stable under natural assumptions.

{The conjugate gradient method applied to the normal equations ATAx=ATb (CGLS) is often used for solving large sparse linear least squares problems. The mathematically equivalent algorithm LSQR based on the Lanczos bidiagonalization process is an often recommended alternative. In this paper, the achievable accuracy of different conjgate gradient and Lanczos methods in finite precision is studied. It is shown that an implementation of algorithm CGLS in which the residual sk=AT(b-Axk) of the normal equations is recurred will not in general achieve accurate solutions. The same conclusion holds for the method based on Lanczos bidiagonalization with starting vector ATb. For the preferred implementation of CGLS we bound the error ||r-rk|| of the computed residual rk. Numerical tests are given that confirm a conjecture of backward stability. The achievable accuracy of LSQR is shown to be similar. The analysis essentially also covers the preconditioned case.

A matrix-free algorithm, IRLANB, for the efficient computation of the smallest singular triplets of large and possibly sparse matrices is described. Key characteristics of the approach are its use of Lanczos bidiagonalization, implicit restarting, and harmonic Ritz values. The algorithm also uses a deflation strategy that can be applied directly on Lanczos bidiagonalization. A refinement postprocessing phase is applied to the converged singular vectors. The computational costs of the above techniques are kept small as they make direct use of the bidiagonal form obtained in the course of the Lanczos factorization. Several numerical experiments with the method are presented that illustrate its effectiveness and indicate that it performs well compared to existing codes.