In this paper we develop new Newton and conjugate gradient algorithms on the Grassmann and Stiefel manifolds. These manifolds represent the constraints that arise in such areas as the symmetric eigenvalue problem, nonlinear eigenvalue problems, electronic structures computations, and signal processing. In addition to the new algorithms, we show how the geometrical framework gives penetrating new insights allowing us to create, understand, and compare algorithms. The theory proposed here provides a taxonomy for numerical linear algebra algorithms that provide a top level mathematical view of previously unrelated algorithms. It is our hope that developers of new algorithms and perturbation theories will benefit from the theory, methods, and examples in this paper.

. We consider least-squares problems where the coefficient matrices A; b are unknown-butbounded. We minimize the worst-case residual error using (convex) second-order cone programming, yielding an algorithm with complexity similar to one singular value decomposition of A. The method can be interpreted as a Tikhonov regularization procedure, with the advantage that it provides an exact bound on the robustness of solution, and a rigorous way to compute the regularization parameter. When the perturbation has a known (e.g., Toeplitz) structure, the same problem can be solved in polynomial-time using semidefinite programming (SDP). We also consider the case when A; b are rational functions of an unknown-but-bounded perturbation vector. We show how to minimize (via SDP) upper bounds on the optimal worst-case residual. We provide numerical examples, including one from robust identification and one from robust interpolation. Key Words. Least-squares, uncertainty, robustness, second-order cone...

The L-curve is a log-log plot of the norm of a regularized solution versus the norm of the corresponding residual norm. It is a convenient graphical tool for displaying the trade-off between the size of a regularized solution and its fit to the given data, as the regularization parameter varies. The L-curve thus gives insight into the regularizing properties of the underlying regularization method, and it is an aid in choosing an appropriate regularization parameter for the given data. In this chapter we summarize the main properties of the L-curve, and demonstrate by examples its usefulness and its limitations both as an analysis tool and as a method for choosing the regularization parameter. 1 Introduction Practically all regularization methods for computing stable solutions to inverse problems involve a trade-off between the "size" of the regularized solution and the quality of the fit that it provides to the given data. What distinguishes the various regularization methods is how...

. The L-curve criterion is often applied to determine a suitable value of the regularization parameter when solving ill-conditioned linear systems of equations with a right-hand side contaminated by errors of unknown norm. However, the computation of the L-curve is quite costly for large problems; the determination of a point on the L-curve requires that both the norm of the regularized approximate solution and the norm of the corresponding residual vector be available. Therefore, usually only a few points on the L-curve are computed, and these values rather than the L-curve, are used to determine a value of the regularization parameter. We propose a new approach to determine a value of the regularization parameter based on computing an L-ribbon that contains the L-curve in its interior. An L-ribbon can be computed fairly inexpensively by partial Lanczos bidiagonalization of the matrix of the given linear system of equations. A suitable value of the regularization parameter is then det...

An optimization problem that does not have an unique local minimum is often very difficult to solve. For a nonlinear least squares problem this is the case when the Jacobian is rank deficient in a neighborhood of a local minimum. Moreover, a Gauss-Newton method such as Levenberg-Marquardt will have very slow convergence for such a problem. We analyze these problems where the Jacobian is rank deficient and suggest other problem formulations more suitable for Gauss-Newton methods. The two methods we propose are a truncated Gauss-Newton method and a Gauss-Newton method based on the Tikhonov regularized nonlinear least squares problem. We test the methods on artificial problems where the rank of the Jacobian and the nonlinearity of the problem may be chosen making it possible to show the different features of the problem and the methods. The conclusion from the analysis and the tests is that the two methods have similar local convergence properties. The method based on Tikhononv regulari...

We consider the problem of computing the solution of large-scale discrete ill-posed problems when there is noise in the data. These problems arise in important areas such as seismic inversion, medical imaging and signal processing. We pose the problem as a quadratically constrained least squares problem and develop a method for the solution of such problem. Our method does not require factorization of the coefficient matrix, it has very low storage requirements and handles the high degree of singularities arising in discrete ill-posed problems. We present numerical results on test problems and an application of the method to a practical problem with real data.

. First-kind Volterra problems arise in numerous applications, from inverse problems in mathematical biology to inverse heat conduction problems. Unfortunately, such problems are also ill-posed due to lack of continuous dependence of solutions on data. Consequently, numerical methods to solve first-kind Volterra equations are only effective when regularizing features are built into the algorithms or used to control stepsize. Classical methods often combine numerical discretization with Tikhonov regularization, but in doing so the underlying Volterra (or causal) nature of the original problem is often destroyed. Instead, a "predictor-corrector" type of numerical method is proposed which combines at each step "local regularization" ideas with the use of small intervals of future data. The result is a regularized numerical method which retains much of the causal nature of the Volterra problem and may be solved in fast sequential steps, often improving upon the performance of classical alg...

Abstract We present a new method for solving large-scale quadratic problems with quadratic and nonnegativity constraints. Such problems arise for example in the regularization of ill-posed problems in image restoration where, in addition, some of the matrices involved are very ill-conditioned. The new method uses recently developed techniques for the large-scale trust-region subproblem.

Image restoration applications often result in ill-posed least squares problems involving large, structured matrices. One approach used extensively is to restore the image in the frequency domain, thus providing fast algorithms using ffts. This is equivalent to using a circulant approximation to a given matrix. Iterative methods may also be used effectively by exploiting the structure of the matrix. While iterative schemes are more expensive than fft-based methods, it has been demonstrated that they are capable of providing better restorations. As an alternative, we propose an approximate singular value decomposition, which can be used in a variety of applications. Used as a direct method, the computed restorations are comparable to iterative methods but are computationally less expensive. In addition, the approximate svd may be used with the generalized cross validation method to choose regularization parameters. It is also demonstrated that the approximate svd can be an ef...

An orthogonal Procrustes problem on the Stiefel manifold is studied, where a matrix Q with orthonormal columns is to be found that minimizes kAQ \Gamma BkF for an l \Theta m matrix A and an l \Theta n matrix B with l m and m ? n. Based on the normal and secular equations and the properties of the Stiefel manifold, necessary conditions for a global minimum, as well as necessary and sufficient conditions for a local minimum, are derived. Key words. constraint, Lagrange multiplier, least squares problems, normal equations, orthogonal Procrustes problem, secular equations, singular value decomposition, Stiefel manifold Department of Mathematics, Linkoping University, S-581 83 Linkoping, Sweden. Email: elden@math.liu.se. y Computer Science Department, University of Minnesota, Minneapolis, MN 55455, U.S.A. E-mail: hpark@cs.umn.edu. The work of this author was supported in part by the National Science Foundation grant CCR-9507307. 1 Introduction Let two matrices A 0 2 R l\Thetam w...