Iterative refinement for constrained and weighted linear least squares

We present an algorithm for mixed precision iterative refinement on the constrained and weighted linear least squares problem, the CWLSQ problem. The approximate solution is obtained by solving the CWLSQ problem with the weightedQR factorization [6]. With backward errors for the weightedQR decomposition together with perturbation bounds for the CWLSQ problem we analyze the convergence behaviour of the iterative refinement procedure.In the unweighted case the initial convergence rate of the error of the iteratively refined solution is determined essentially by the condition number. For the CWLSQ problem the initial convergence behaviour is more complicated. The analysis shows that the initial convergence is dependent both on the condition of the problem related to the solution,x, and the vector =Wr, whereW is the weight matrix andr is the residual.We test our algorithm on two examples where the solution is known and the condition number of the problem can be varied. The computational test confirms the theoretical results and verifies that mixed precision iterative refinement, using the system matrix and the weightedQR decomposition, is an effective way of improving an approximate solution to the CWLSQ problem.     

Solution of augmented linear systems using orthogonal factorizations

We consider solvingx+Ay=b andA ^T x=c for givenb, c andm ×n A of rankn. This is called the augmented system formulation (ASF) of two standard optimization problems, which include as special cases the minimum 2-norm of a linear underdetermined system (b=0) and the linear least squares problem (c=0), as well as more general problems. We examine the numerical stability of methods (for the ASF) based on the QR factorization ofA, whether by Householder transformations, Givens rotations, or the modified Gram-Schmidt (MGS) algorithm, and consider methods which useQ andR, or onlyR. We discuss the meaning of stability of algorithms for the ASF in terms of stability of algorithms for the underlying optimization problems.We prove the backward stability of several methods for the ASF which useQ andR, inclusing a new one based on MGS, and also show under what circumstances they may be regarded as strongly stable. We show why previous methods usingQ from MGS were not backward stable, but illustrate that some of these methods may be acceptable-error stable. We point out that the numerical accuracy of methods that do not useQ does not depend to any significant extent on which of of the above three QR factorizations is used. We then show that the standard methods which do not useQ are not backward stable or even acceptable-error stable for the general ASF problem, and discuss how iterative refinement can be used to counteract these deficiencies.Dedicated to Carl-Eric Fröberg on the occasion of his 75th birthdayThis research was partially supported by NSERC of Canada Grant No. A9236.     

The structured sensitivity of Vandermonde-like systems

Summary We consider a general class of structured matrices that includes (possibly confluent) Vandermonde and Vandermonde-like matrices. Here the entries in the matrix depend nonlinearly upon a vector of parameters. We define, condition numbers that measure the componentwise sensitivity of the associated primal and dual solutions to small componentwise perturbations in the parameters and in the right-hand side. Convenient expressions are derived for the infinity norm based condition numbers, and order-of-magnitude estimates are given for condition numbers defined in terms of a general vector norm. We then discuss the computation of the corresponding backward errors. After linearising the constraints, we derive an exact expression for the infinity norm dual backward error and show that the corresponding primal backward error is given by the minimum infinity-norm solution of an underdetermined linear system. Exact componentwise condition numbers are also derived for matrix inversion and the least squares problem, and the linearised least squares backward error is characterised.     

Complexity reduction of least squares problems involving special vandermonde matrices

This paper presents a new QRD factorization of a rectangular Vandermonde matrix for a special point distribution, including the symmetric case, based on ak-dimensional block decomposition of the matrix and some properties of the Kronecker product. The computational reduction factor with respect to any QR method isk ², in the general case, and 4 in the symmetric one. By the resulting matrix factorization, new formulas are devised for the least squares system solution, whose implementation produces an algorithm of reduced computational cost and computer storage. Finally the perturbation bounds of this new factorization are devised.     

Computation of orthonormal factors for fundamental solution matrices

Summary. In this work, we introduce and analyze two new techniques for obtaining the Q factor in the QR factorization of some (or all) columns of a fundamental solution matrix Y of a linear differential system. These techniques are based on elementary Householder and Givens transformations. We implement and compare these new techniques with existing approaches on some examples. Received October 27, 1997 / Revised version received September 21, 1998 / Published online August 19, 1999     

A stationary iterative pseudoinverse algorithm

Iterative methods applied to the normal equationsA ^T Ax=A ^T b are sometimes used for solving large sparse linear least squares problems. However, when the matrix is rank-deficient many methods, although convergent, fail to produce the unique solution of minimal Euclidean norm. Examples of such methods are the Jacobi and SOR methods as well as the preconditioned conjugate gradient algorithm. We analyze here an iterative scheme that overcomes this difficulty for the case of stationary iterative methods. The scheme combines two stationary iterative methods. The first method produces any least squares solution whereas the second produces the minimum norm solution to a consistent system. This work was supported by the Swedish Research Council for Engineering Sciences, TFR.     

Constrained least squares methods for linear timevarying DAE systems

Summary This paper presents a family of methods for accurate solution of higher index linear variable DAE systems, . These methods use the DAE system and some of its first derivatives as constraints to a least squares problem that corresponds to a Taylor series ofy, or an approximative equality derived from a Pade' approximation of the exponential function. Accuracy results for systems transformable to standard canonical form are given. Advantages, disadvantages, stability properties and implementation of these methods are discussed and two numerical examples are given, where we compare our results with results from more traditional methods.     

Computing Row and Column Counts for Sparse QR and LU Factorization

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.This revised version was published online in October 2005 with corrections to the Cover Date.     

Componentwise error analysis for the block LU factorization of totally nonnegative matrices

For the first time, perturbation bounds including componentwise perturbation bounds for the block LU factorization have been provided by Dopico and Molera (2005) [5]. In this paper, componentwise error analysis is presented for computing the block LU factorization of nonsingular totally nonnegative matrices. We present a componentwise bound on the equivalent perturbation for the computed block LU factorization. Consequently, combining with the componentwise perturbation results we derive componentwise forward error bounds for the computed block factors.     

Reliable parallel solution of bidiagonal systems

Summary. This paper presents a new efficient algorithm for solving bidiagonal systems of linear equations on massively parallel machines. We use a divide and conquer approach to compute a representative subset of the solution components after which we solve the complete system in parallel with no communication overhead. We address the numerical properties of the algorithm in two ways: we show how to verify the à posteriori backward stability at virtually no additional cost, and prove that the algorithm is à priori forward stable. We then show how we can use the algorithm in order to bound the possible perturbations in the solution components. Received March 13, 1998 / Revised version received December 21, 1999 / Published online June 20, 2001     

Numerical behaviour of the modified gram-schmidt GMRES implementation

In [6] the Generalized Minimal Residual Method (GMRES) which constructs the Arnoldi basis and then solves the transformed least squares problem was studied. It was proved that GMRES with the Householder orthogonalization-based implementation of the Arnoldi process (HHA), see [9], is backward stable. In practical computations, however, the Householder orthogonalization is too expensive, and it is usually replaced by the modified Gram-Schmidt process (MGSA). Unlike the HHA case, in the MGSA implementation the orthogonality of the Arnoldi basis vectors is not preserved near the level of machine precision. Despite this, the MGSA-GMRES performs surprisingly well, and its convergence behaviour and the ultimately attainable accuracy do not differ significantly from those of the HHA-GMRES. As it was observed, but not explained, in [6], it is thelinear independence of the Arnoldi basis, not the orthogonality near machine precision, that is important. Until the linear independence of the basis vectors is nearly lost, the norms of the residuals in the MGSA implementation of GMRES match those of the HHA implementation despite the more significant loss of orthogonality. In this paper we study the MGSA implementation. It is proved that the Arnoldi basis vectors begin to lose their linear independenceonly after the GMRES residual norm has been reduced to almost its final level of accuracy, which is proportional to κ(A)ε, where κ(A) is the condition number ofA and ε is the machine precision. Consequently, unless the system matrix is very ill-conditioned, the use of the modified Gram-Schmidt GMRES is theoretically well-justified. This work was supported by the GA AS CR under grants 230401 and A 2030706, by the NSF grant Int 921824, and by the DOE grant DEFG0288ER25053. Part of the work was performed while the third author visited Department of Mathematics and Computer Science, Emory University, Atlanta, USA.     

Preconditionings and splittings for rectangular systems

Summary Recently Eiermann, Marek, and Niethammer have shown how to applysemiiterative methods to a fixed point systemx=Tx+c which isinconsistent or in whichthe powers of the fixed point operator T have no limit, to obtain iterative methods which converge to some approximate solution to the fixed point system. In view of their results we consider here stipulations on apreconditioning QAx=Qb of the systemAx=b and, separately, on asplitting A=M–N which lead to fixed point systems such that, with the aid of a semiiterative method, the iterative scheme converges to a weighted Moore-Penrose solution to the systemAx=b. We show in several ways that to obtain a meaningful limit point from a semiiterative method requires less restrictions on the splittings or the reconditionings than those which have been required in the classical Picard iterative method (see, e.g., the works of Berman and Plemmons, Berman and Neumann, and Tanabe).We pay special attention to the case when the weighted Moore-Penrose solution which is sought is the minimal norm least squares solution toAx=b.Research supported by the Deutsche ForschungsgemeinschaftPartially supported by AFOSR research grant 88-0047     

On the modified Gram-Schmidt algorithm for weighted and constrained linear least squares problems

A framework and an algorithm for using modified Gram-Schmidt for constrained and weighted linear least squares problems is presented. It is shown that a direct implementation of a weighted modified Gram-Schmidt algorithm is unstable for heavily weighted problems. It is shown that, in most cases it is possible to get a stable algorithm by a simple modification free from any extra computational costs. In particular, it is not necessary to perform reorthogonalization.Solving the weighted and constrained linear least squares problem with the presented weighted modified Gram-Schmidt algorithm is seen to be numerically equivalent to an algorithm based on a weighted Householder-likeQR factorization applied to a slightly larger problem. This equivalence is used to explain the instability of the weighted modified Gram-Schmidt algorithm. If orthogonality, with respect to a weighted inner product, of the columns inQ is important then reorthogonalization can be used. One way of performing such reorthogonalization is described.Computational tests are given to show the main features of the algorithm.     

Parallel factorization of banded linear matrices using a systolic array processor

A parallel algorithm for calculating theQR factorization of a banded system of linear equations using a systolic array processor is presented and an application to spline fitting is given. The major advantage of the method is that the size of the processor array is fixed by the size of the bandwidth of the system to be solved. This allows the factorization of large systems using small systolic arrays. The cost of the method, in terms of storage and time, is optimal.     

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.     

Bases for kernel-based spaces

Since it is well-known (De Marchi and Schaback (2001) [4]) that standard bases of kernel translates are badly conditioned while the interpolation itself is not unstable in function space, this paper surveys the choices of other bases. All data-dependent bases turn out to be defined via a factorization of the kernel matrix defined by these data, and a discussion of various matrix factorizations (e.g. Cholesky, QR, SVD) provides a variety of different bases with different properties. Special attention is given to duality, stability, orthogonality, adaptivity, and computational efficiency. The “Newton” basis arising from a pivoted Cholesky factorization turns out to be stable and computationally cheap while being orthonormal in the “native” Hilbert space of the kernel. Efficient adaptive algorithms for calculating the Newton basis along the lines of orthogonal matching pursuit conclude the paper.     

Improved componentwise verified error bounds for least squares problems and underdetermined linear systems

Recently Miyajima presented algorithms to compute componentwise verified error bounds for the solution of full-rank least squares problems and underdetermined linear systems. In this paper we derive simpler and improved componentwise error bounds which are based on equalities for the error of a given approximate solution. Equalities are not improvable, and the expressions are formulated in a way that direct evaluation yields componentwise and rigorous estimates of good quality. The computed bounds are correct in a mathematical sense covering all sources of errors, in particular rounding errors. Numerical results show a gain in accuracy compared to previous results.     

Extra Precise Preconditioning for Non-Hermitian Eigenvalue Problems

For badly conditioned Non-Hermite eigenvalue problems known iterative refinement procedures often fail. It is proposed to use extra precise iterative refinement in a precondition step. As preconditioners the matrix Q from the QR decomposition or at least the matrix of - even not very accurate - eigenvalues can be used. (© 2006 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim)     

A generalized successive overrelaxation method for least squares problems

In this paper a new iterative method is given for solving large sparse least squares problems and computing the minimum norm solution to underdetermined consistent linear systems. The new scheme is called the generalized successive overrelaxation (GSOR) method and is shown to be convergent ifA is full column rank. The GSOR method involves a parameter ρ and an auxiliary matrixP. One can choose matrix P so that the GSOR method only involves matrix and vector operations; therefore the GSOR method is suitable for parallel computations. Besides, the GSOR method can be combined with preconditioning techniques, and therefore can be expected to be more effective. This author's work was supported by Natural Science Foundation of Liaoning Province, China.     

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||₂.