Solution Of Least Squares Problem Research Articles

Linearly convergent adjoint free solution of least squares problems by random descent

We consider the problem of solving linear least squares problems in a framework where only evaluations of the linear map are possible. We derive randomized methods that do not need any other matrix operations than forward evaluations, especially no evaluation of the adjoint map is needed. Our method is motivated by the simple observation that one can get an unbiased estimate of the application of the adjoint. We show convergence of the method and then derive a more efficient method that uses an exact linesearch. This method, called random descent, resembles known methods in other context and has the randomized coordinate descent method as special case. We provide convergence analysis of the random descent method emphasizing the dependence on the underlying distribution of the random vectors. Furthermore we investigate the applicability of the method in the context of ill-posed inverse problems and show that the method can have beneficial properties when the unknown solution is rough. We illustrate the theoretical findings in numerical examples. One particular result is that the random descent method actually outperforms established transposed-free methods (TFQMR and CGS) in examples.

Global optimization for sparse solution of least squares problems

Finding solutions to least-squares problems with low cardinality has found many applications, including portfolio optimization, subset selection in statistics, and inverse problems in signal processing. Although most works consider local approaches that scale with high-dimensional problems, some others address its global optimization via mixed integer programming (MIP) reformulations. We propose dedicated branch-and-bound methods for the exact resolution of moderate-size, yet difficult, sparse optimization problems, through three possible formulations: cardinality-constrained and cardinality-penalized least-squares, and cardinality minimization under quadratic constraints. A specific tree exploration strategy is built. Continuous relaxation problems involved at each node are reformulated as -norm-based optimization problems, for which a dedicated algorithm is designed. The obtained certified solutions are shown to better estimate sparsity patterns than standard methods on simulated variable selection problems involving highly correlated variables. Problem instances selecting up to 24 components among 100 variables, and up to 15 components among 1000 variables, can be solved in less than 1000 s. Unguaranteed solutions obtained by limiting the computing time to 1s are also shown to provide competitive estimates. Our algorithms strongly outperform the CPLEX MIP solver as the dimension increases, especially for quadratically-constrained problems. The source codes are made freely available online.

Probabilistic error estimation for non-intrusive reduced models learned from data of systems governed by linear parabolic partial differential equations

This work derives a residual-based a posteriori error estimator for reduced models learned with non-intrusive model reduction from data of high-dimensional systems governed by linear parabolic partial differential equations with control inputs. It is shown that quantities that are necessary for the error estimator can be either obtained exactly as the solutions of least-squares problems in a non-intrusive way from data such as initial conditions, control inputs, and high-dimensional solution trajectories or bounded in a probabilistic sense. The computational procedure follows an offline/online decomposition. In the offline (training) phase, the high-dimensional system is judiciously solved in a black-box fashion to generate data and to set up the error estimator. In the online phase, the estimator is used to bound the error of the reduced-model predictions for new initial conditions and new control inputs without recourse to the high-dimensional system. Numerical results demonstrate the workflow of the proposed approach from data to reduced models to certified predictions.

Algebraic techniques for least squares problems in commutative quaternionic theory

Due to the rise of commutative quaternion in Hopfield neural networks, digital signal, and image processing, one encounters the approximate solution problems of the commutative quaternion linear equations and . This paper, by means of real representation and complex representation of commutative quaternion matrices, introduces concepts of norms of commutative quaternion matrices and derives two algebraic techniques for finding solutions of least squares problems in commutative quaternionic theory.

Condition numbers for the K-weighted pseudoinverse and their statistical estimation

The explicit expressions of normwise, mixed and componentwise condition numbers for the K-weighted pseudoinverse are first presented. With the intermediate result, i.e. the derivative of , we can recover the explicit expressions of condition numbers for solution of least squares problem with equality constraint. Then, we investigate the statistical estimation of condition numbers of using the probabilistic spectral norm estimator and the small-sample statistical condition estimation method, and devise three algorithms. Finally, numerical examples are provided to illustrate these results.

Total variation superiorized conjugate gradient method for image reconstruction

The conjugate gradient (CG) method is commonly used for the relatively-rapid solution of least squares problems. In image reconstruction, the problem can be ill-posed and also contaminated by noise; due to this, approaches such as regularization should be utilized. Total variation (TV) is a useful regularization penalty, frequently utilized in image reconstruction for generating images with sharp edges. When a non-quadratic norm is selected for regularization, as is the case for TV, then it is no longer possible to use CG. Non-linear CG is an alternative, but it does not share the efficiency that CG shows with least squares and methods such as fast iterative shrinkage-thresholding algorithms (FISTA) are preferred for problems with TV norm. A different approach to including prior information is superiorization. In this paper it is shown that the conjugate gradient method can be superiorized. Five different CG variants are proposed, including preconditioned CG. The CG methods superiorized by the total variation norm are presented and their performance in image reconstruction is demonstrated. It is illustrated that some of the proposed variants of the superiorized CG method can produce reconstructions of superior quality to those produced by FISTA and in less computational time, due to the speed of the original CG for least squares problems. In the Appendix we examine the behavior of one of the superiorized CG methods (we call it S-CG); one of its input parameters is a positive number ε. It is proved that, for any given ε that is greater than the half-squared-residual for the least squares solution, S-CG terminates in a finite number of steps with an output for which the half-squared-residual is less than or equal to ε. Importantly, it is also the case that the output will have a lower value of TV than what would be provided by unsuperiorized CG for the same value ε of the half-squared residual.

FPGA Implementation of Orthogonal Matching Pursuit for Compressive Sensing Reconstruction

In this paper, we present a novel architecture based on field-programmable gate arrays (FPGAs) for the reconstruction of compressively sensed signal using the orthogonal matching pursuit (OMP) algorithm. We have analyzed the computational complexities and data dependence between different stages of OMP algorithm to design its architecture that provides higher throughput with less area consumption. Since the solution of least square problem involves a large part of the overall computation time, we have suggested a parallel low-complexity architecture for the solution of the linear system. We have further modeled the proposed design using Simulink and carried out the implementation on FPGA using Xilinx system generator tool. We have presented here a methodology to optimize both area and execution time in Simulink environment. The execution time of the proposed design is reduced by maximizing parallelism by appropriate level of unfolding, while the FPGA resources are reduced by sharing the hardware for matrix–vector multiplication across the data-dependent sections of the algorithm. The hardware implementation on the Virtex6 FPGA provides significantly superior performance in terms of resource utilization measured in the number of occupied slices, and maximum usable frequency compared with the existing implementations. Compared with the existing similar design, the proposed structure involves 328 more DSP48s, but it involves $25\,802$ less slices and 1.85 times less computation time for signal reconstruction with $N = 1024$ , $ K = 256$ , and $m = 36$ , where $N$ is the number of samples, $K$ is the size of the measurement vector, and $m$ is the sparsity. It also provides a higher peak signal-to-noise ratio value of 38.9 dB with a reconstruction time of $0.34~\mu $ s, which is twice faster than the existing design. In addition, we have presented a performance metric to implement the OMP algorithm in resource constrained FPGA for the better quality of signal reconstruction.

Componentwise enclosure for solutions of least squares problems and underdetermined systems

Algorithms for enclosing solutions of least squares problems and underdetermined systems are proposed. The results obtained by these algorithms are “verified” in the sense that all the possible rounding errors have been taken into account. In order to develop these algorithms, theories for obtaining componentwise error bounds of numerical solutions are established. The error bounds by the proposed algorithms are vectors, as opposed that those by the fast algorithms proposed by Rump are scalars. The proposed algorithms require computational costs similar to those of Rumpʼs algorithms. It is moreover proved that each component of the error bounds by the proposed algorithms is equal or smaller than the error bounds by Rumpʼs algorithms. Numerical results show the properties of the algorithms.

Accurate Solution of Structured Least Squares Problems via Rank-Revealing Decompositions

Least squares problems $\min_x \|b - Ax\|_2$ where the matrix $A\in \mathbb{C}^{m\times n}$ ($m\geq n$) has some particular structure arise frequently in applications. Polynomial data fitting is a well-known instance of problems that yield highly structured matrices, but many other examples exist. Very often, structured matrices have huge condition numbers $\kappa_2 (A) = \|A\|_2 \, \|A^\dagger\|_2$ ($A^\dagger$ is the Moore--Penrose pseudoinverse of $A$) and therefore standard algorithms fail to compute accurate minimum 2-norm solutions of least squares problems. In this work, we introduce a framework that allows us to compute minimum 2-norm solutions of many classes of structured least squares problems accurately, i.e., with errors $\| \widehat{x}_0 - x_0 \|_2 / \|x_0 \|_2 = O({\tt u})$, where ${\tt u}$ is the unit roundoff, independently of the magnitude of $\kappa_2 (A)$ for most vectors $b$. The cost of these accurate computations is $O(n^2 m)$ flops, i.e., roughly the same cost as standard algorithm...

Least Squares Symmetrizable Solutions for a Class of Matrix Equations

In this paper, we discuss least squares symmetrizable solutions of matrix equations (AX = B, XC = D) and its optimal approximation solution. With the matrix row stacking, Kronecker product and special relations between two linear subspaces are topological isomorphism, and we derive the general solutions of least squares problem. With the invariance of the Frobenius norm under orthogonal transformations, we obtain the unique solution of optimal approximation problem. In addition, we present an algorithm and numerical experiment to obtain the optimal approximation solution.

Numerical investigations of linear least squares methods for derivative estimation

The results of a numerical investigation into the errors for least squares estimates of function gradients are presented. The underlying algorithm is obtained by constructing a least squares problem using a truncated Taylor expansion. An error bound associated with this method contains in its numerator terms related to the Taylor series remainder, while its denominator contains the smallest singular value of the least squares matrix. Perhaps for this reason the error bounds are often found to be pessimistic by several orders of magnitude. The circumstance under which these poor estimates arise is elucidated and an empirical correction of the theoretical error bounds is conjectured and investigated numerically. This is followed by an indication of how the conjecture is supported by a rigorous argument. References J. A. Belward, I. W. Turner, and M. Ilic. On derivative estimation and the solution of least squares problems. Journal of Computational and Applied Mathematics, 222:511--523, 2008. C. L. Lawson and R. J. Hansen. Solving Least Squares Problems. Prentice--Hall, 1974. P. Lancaster and K. Salkauskas. Curve and Surface Fitting, An Introduction. Academic Press, San Diego, 1986. M. N. Oqielat, J. A. Belward, I. W. Turner, and B. I. Loch. {A Hybrid Clough-Tocher Radial Basis Function Method for Modelling Leaf Surfaces}. In Oxley, L., and Kulasiri, D., (eds) MODSIM 2007 International Congress on Modelling and Simulation, pages 400--406, 2007. http://mssanz.org.au/MODSIM2007/papers/7_s50 M. N. Oqielat, I. W. Turner, and J. A. Belward. A hybrid clough-tocher method for surface fitting with application to leaf data. Applied Mathematical Modeling, 2008. doi:10.1016/j.apm.2008.07.023 I. W. Turner, J. A. Belward, and M. N. Oqielat. Error bounds for least squares gradient estimates. In preparation, 2008.

On derivative estimation and the solution of least squares problems

Surface interpolation finds application in many aspects of science and technology. Two specific areas of interest are surface reconstruction techniques for plant architecture and approximating cell face fluxes in the finite volume discretisation strategy for solving partial differential equations numerically. An important requirement of both applications is accurate local gradient estimation. In surface reconstruction this gradient information is used to increase the accuracy of the local interpolant, while in the finite volume framework accurate gradient information is essential to ensure second order spatial accuracy of the discretisation. In this work two different least squares strategies for approximating these local gradients are investigated and the errors associated with each analysed. It is shown that although the two strategies appear different, they produce the same least squares error. Some carefully chosen case studies are used to elucidate this finding.

Expansion of weighted pseudoinverse matrices with singular weights into matrix power products and iteration methods

We obtain expansions of weighted pseudoinverse matrices with singular weights into matrix power products with negative exponents and arbitrary positive parameters. We show that the rate of convergence of these expansions depends on a parameter. On the basis of the proposed expansions, we construct and investigate iteration methods with quadratic rate of convergence for the calculation of weighted pseudoinverse matrices and weighted normal pseudosolutions. Iteration methods for the calculation of weighted normal pseudosolutions are adapted to the solution of least-squares problems with constraints.

Fast Algorithms for Structured Least Squares and Total Least Squares Problems.

We consider the problem of solving least squares problems involving a matrix M of small displacement rank with respect to two matrices Z1 and Z2. We develop formulas for the generators of the matrix M HM in terms of the generators of M and show that the Cholesky factorization of the matrix M HM can be computed quickly if Z1 is close to unitary and Z2 is triangular and nilpotent. These conditions are satisfied for several classes of matrices, including Toeplitz, block Toeplitz, Hankel, and block Hankel, and for matrices whose blocks have such structure. Fast Cholesky factorization enables fast solution of least squares problems, total least squares problems, and regularized total least squares problems involving these classes of matrices.

Dynamic programming and minimal norm solutions of least squares problems

Dynamic programming and minimal norm solutions of least squares problems

Approximating minimum norm solutions of rank-deficient least squares problems

In this paper we consider the solution of linear least squares problems minx∥Ax - b∥22 where the matrix A ∈ Rm × n is rank deficient. Put p = min{m, n}, let σi, i = 1, 2,…, p, denote the singular values of A, and let ui and vi denote the corresponding left and right singular vectors. Then the minimum norm solution of the least squares problem has the form x* = ∫ri = 1(uTib/σi)vi, where r ≤ p is the rank of A. The Riley–Golub iteration, xk + 1 = arg minx{∥Ax - b∥22 + λ∥x − xk∥22} converges to the minimum norm solution if x0 is chosen equal to zero. The iteration is implemented so that it takes advantage of a bidiagonal decomposition of A. Thus modified, the iteration requires only O(p) flops (floating point operations). A further gain of using the bidiagonalization of A is that both the singular values σi and the scalar products uTib can be computed at marginal extra cost. Moreover, we determine the regularization parameter, λ, and the number of iterations, k, in a way that minimizes the difference x* − xk with respect to a certain norm. Explicit rules are derived for calculating these parameters. One advantage of our approach is that the numerical rank can be easily determined by using the singular values. Furthermore, by the iterative procedure, x* is approximated without computing the singular vectors of A. This gives a fast and reliable method for approximating minimum norm solutions of well-conditioned rank-deficient least squares problems. Numerical experiments illustrate the viability of our ideas, and demonstrate that the new method gives more accurate approximations than an approach based on a QR decomposition with column pivoting. © 1998 John Wiley & Sons, Ltd.

Approximating minimum norm solutions of rank‐deficient least squares problems

In this paper we consider the solution of linear least squares problems minx∥Ax - b∥22 where the matrix A ∈ Rm × n is rank deficient. Put p = min{m, n}, let σi, i = 1, 2,…, p, denote the singular values of A, and let ui and vi denote the corresponding left and right singular vectors. Then the minimum norm solution of the least squares problem has the form x* = ∫ri = 1(uTib/σi)vi, where r ≤ p is the rank of A. The Riley–Golub iteration, xk + 1 = arg minx{∥Ax - b∥22 + λ∥x − xk∥22} converges to the minimum norm solution if x0 is chosen equal to zero. The iteration is implemented so that it takes advantage of a bidiagonal decomposition of A. Thus modified, the iteration requires only O(p) flops (floating point operations). A further gain of using the bidiagonalization of A is that both the singular values σi and the scalar products uTib can be computed at marginal extra cost. Moreover, we determine the regularization parameter, λ, and the number of iterations, k, in a way that minimizes the difference x* − xk with respect to a certain norm. Explicit rules are derived for calculating these parameters. One advantage of our approach is that the numerical rank can be easily determined by using the singular values. Furthermore, by the iterative procedure, x* is approximated without computing the singular vectors of A. This gives a fast and reliable method for approximating minimum norm solutions of well-conditioned rank-deficient least squares problems. Numerical experiments illustrate the viability of our ideas, and demonstrate that the new method gives more accurate approximations than an approach based on a QR decomposition with column pivoting. © 1998 John Wiley & Sons, Ltd.

A Preconditioned Krylov Subspace Method for the Solution of Least Squares Problems in Inverse Scattering

We present an iterative method of preconditioned Krylov type for the solution of large least squares problems. We prove that the method is robust and investigate its rate of convergence. For an important application, originating from seismic inverse scattering, we derive a suitable preconditioner using asymptotic theory. Numerical experiments are used to compare the method with other iterative methods. It appears that the preconditioned Krylov method can be much more efficient than CG applied to the normal equations.

Iterative refinement enhances the stability ofQR factorization methods for solving linear equations

Iterative refinement is a well-known technique for improving the quality of an approximate solution to a linear system. In the traditional usage residuals are computed in extended precision, but more recent work has shown that fixed precision is sufficient to yield benefits for stability. We extend existing results to show that fixed precision iterative refinement renders anarbitrary linear equations solver backward stable in a strong, componentwise sense, under suitable assumptions. Two particular applications involving theQR factorization are discussed in detail: solution of square linear systems and solution of least squares problems. In the former case we show that one step of iterative refinement suffices to produce a small componentwise relative backward error. Our results are weaker for the least squares problem, but again we find that iterative refinement improves a componentwise measure of backward stability. In particular, iterative refinement mitigates the effect of poor row scaling of the coefficient matrix, and so provides an alternative to the use of row interchanges in the HouseholderQR factorization. A further application of the results is described to fast methods for solving Vandermonde-like systems.