A subspace implementation of quasi-Newton trust region methods for unconstrained optimization

This paper studies subspace properties of trust region methods for unconstrained optimization, assuming the approximate Hessian is updated by quasi- Newton formulae and the initial Hessian approximation is appropriately chosen. It is shown that the trial step obtained by solving the trust region subproblem is in the subspace spanned by all the gradient vectors computed. Thus, the trial step can be defined by minimizing the quasi-Newton quadratic model in the subspace. Based on this observation, some subspace trust region algorithms are proposed and numerical results are also reported.     

Krylov eigenvalue strategy using the FEAST algorithm with inexact system solves

The FEAST eigenvalue algorithm is a subspace iteration algorithm that uses contour integration to obtain the eigenvectors of a matrix for the eigenvalues that are located in any user‐defined region in the complex plane. By computing small numbers of eigenvalues in specific regions of the complex plane, FEAST is able to naturally parallelize the solution of eigenvalue problems by solving for multiple eigenpairs simultaneously. The traditional FEAST algorithm is implemented by directly solving collections of shifted linear systems of equations; in this paper, we describe a variation of the FEAST algorithm that uses iterative Krylov subspace algorithms for solving the shifted linear systems inexactly. We show that this iterative FEAST algorithm (which we call IFEAST) is mathematically equivalent to a block Krylov subspace method for solving eigenvalue problems. By using Krylov subspaces indirectly through solving shifted linear systems, rather than directly using them in projecting the eigenvalue problem, it becomes possible to use IFEAST to solve eigenvalue problems using very large dimension Krylov subspaces without ever having to store a basis for those subspaces. IFEAST thus combines the flexibility and power of Krylov methods, requiring only matrix–vector multiplication for solving eigenvalue problems, with the natural parallelism of the traditional FEAST algorithm. We discuss the relationship between IFEAST and more traditional Krylov methods and provide numerical examples illustrating its behavior.     

The effect of non-optimal bases on the convergence of Krylov subspace methods

Summary There are many examples where non-orthogonality of a basis for Krylov subspace methods arises naturally. These methods usually require less storage or computational effort per iteration than methods using an orthonormal basis (optimal methods), but the convergence may be delayed. Truncated Krylov subspace methods and other examples of non-optimal methods have been shown to converge in many situations, often with small delay, but not in others. We explore the question of what is the effect of having a non-optimal basis. We prove certain identities for the relative residual gap, i.e., the relative difference between the residuals of the optimal and non-optimal methods. These identities and related bounds provide insight into when the delay is small and convergence is achieved. Further understanding is gained by using a general theory of superlinear convergence recently developed. Our analysis confirms the observed fact that in exact arithmetic the orthogonality of the basis is not important, only the need to maintain linear independence is. Numerical examples illustrate our theoretical results.This revised version was published online in June 2005 due to a typesetting mistake in the footnote on page 7.     

Linear M-estimation with bounded variables

A subproblem in the trust region algorithm for non-linear M-estimation by Ekblom and Madsen is to find the restricted step. It is found by calculating the M-estimator of the linearized model, subject to anL ₂-norm bound on the variables. In this paper it is shown that this subproblem can be solved by applying Hebden-iterations to the minimizer of the Lagrangian function. The new method is compared with an Augmented Lagrange implementation.     

Breakdowns and stagnation in iterative methods

One of the disadvantages of Krylov subspace iterative methods is the possibility of breakdown. This occurs when it is impossible to get the next approximation of the solution to the linear system of equationsAu=f. There are two different situations: lucky breakdown, when we have found the solution and hard breakdown, when the next Krylov subspace cannot be generated and/or the next approximate solution (iterate) cannot be computed. We show that some breakdowns depend on the chosen method of generating the basis vectors. Another undesirable feature of the iterative methods is stagnation. This occurs when the error does not change for several iterative steps. We investigate when iterative methods can stagnate and describe conditions which characterize stagnation. We show that in some cases stagnation can imply breakdown.     

CMRH: A new method for solving nonsymmetric linear systems based on the Hessenberg reduction algorithm

The Generalized Minimal Residual (GMRES) method and the Quasi-Minimal Residual (QMR) method are two Krylov methods for solving linear systems. The main difference between these methods is the generation of the basis vectors for the Krylov subspace. The GMRES method uses the Arnoldi process while QMR uses the Lanczos algorithm for constructing a basis of the Krylov subspace. In this paper we give a new method similar to QMR but based on the Hessenberg process instead of the Lanczos process. We call the new method the CMRH method. The CMRH method is less expensive and requires slightly less storage than GMRES. Numerical experiments suggest that it has behaviour similar to GMRES. This revised version was published online in June 2006 with corrections to the Cover Date.     

An Arnoldi-Extrapolation algorithm for computing PageRank

The Arnoldi-type algorithm proposed by Golub and Greif [G. Golub, C. Greif, An Arnoldi-type algorithm for computing PageRank, BIT 46 (2006) 759-771] is a restarted Krylov subspace method for computing PageRank. However, this algorithm may not be efficient when the damping factor is high and the dimension of the search subspace is small. In this paper, we first develop an extrapolation method based on Ritz values. We then consider how to periodically knit this extrapolation method together with the Arnoldi-type algorithm. The resulting algorithm is the Arnoldi-Extrapolation algorithm. The convergence of the new algorithm is analyzed. Numerical experiments demonstrate the numerical behavior of this algorithm.     

On the construction of restricted-denominator exponential W-methods

In this paper we consider the practical construction of exponential W-methods for the solution of large stiff nonlinear initial value problems, based on the restricted-denominator rational approach for the computation of the functions of matrices required. This approach is employed together with the Krylov subspace method based on the Arnoldi algorithm. Two integrators are constructed and tested on some classical stiff equations arising from the semidiscretization of parabolic problems.     

Augmented conjugate gradient. Application in an iterative process for the solution of scattering problems

We discuss the application of an augmented conjugate gradient to the solution of a sequence of linear systems of the same matrix appearing in an iterative process for the solution of scattering problems. The conjugate gradient method applied to the first system generates a Krylov subspace, then for the following systems, a modified conjugate gradient is applied using orthogonal projections on this subspace to compute an initial guess and modified descent directions leading to a better convergence. The scattering problem is treated via an Exact Controllability formulation and a preconditioned conjugate gradient algorithm is introduced. The set of linear systems to be solved are associated to this preconditioning. The efficiency of the method is tested on different 3D acoustic problems. This revised version was published online in August 2006 with corrections to the Cover Date.     

Hybrid enriched bidiagonalization for discrete ill‐posed problems

The regularizing properties of the Golub–Kahan bidiagonalization algorithm are powerful when the associated Krylov subspace captures the dominating components of the solution. In some applications the regularized solution can be further improved by enrichment, that is, by augmenting the Krylov subspace with a low‐dimensional subspace that represents specific prior information. Inspired by earlier work on GMRES, we demonstrate how to carry these ideas over to the bidiagonalization algorithm, and we describe how to incorporate Tikhonov regularization. This leads to a hybrid iterative method where the choice of regularization parameter in each iteration also provides a stopping rule.     

An ODE-based trust region method for unconstrained optimization problems

In this paper, a new trust region algorithm is proposed for solving unconstrained optimization problems. This method can be regarded as a combination of trust region technique, fixed step-length and ODE-based methods. A feature of this proposed method is that at each iteration, only a system of linear equations is solved to obtain a trial step. Another is that when a trial step is not accepted, the method generates an iterative point whose step-length is defined by a formula. Under some standard assumptions, it is proven that the algorithm is globally convergent and locally superlinear convergent. Preliminary numerical results are reported.     

A global harmonic Arnoldi method for large non-Hermitian eigenproblems with an application to multiple eigenvalue problems

The global Arnoldi method can be used to compute exterior eigenpairs of a large non-Hermitian matrix A, but it does not work well for interior eigenvalue problems. Based on the global Arnoldi process that generates an F-orthonormal basis of a matrix Krylov subspace, we propose a global harmonic Arnoldi method for computing certain harmonic F-Ritz pairs that are used to approximate some interior eigenpairs. We propose computing the F-Rayleigh quotients of the large non-Hermitian matrix with respect to harmonic F-Ritz vectors and taking them as new approximate eigenvalues. They are better and more reliable than the harmonic F-Ritz values. The global harmonic Arnoldi method inherits convergence properties of the harmonic Arnoldi method applied to a larger matrix whose distinct eigenvalues are the same as those of the original given matrix. Some properties of the harmonic F-Ritz vectors are presented. As an application, assuming that A is diagonalizable, we show that the global harmonic Arnoldi method is able to solve multiple eigenvalue problems both in theory and in practice. To be practical, we develop an implicitly restarted global harmonic Arnoldi algorithm with certain harmonic F-shifts suggested. In particular, this algorithm can be adaptively used to solve multiple eigenvalue problems. Numerical experiments show that the algorithm is efficient for the eigenproblem and is reliable for quite ill-conditioned multiple eigenproblems.     

信赖域策略的内点投影既约Hessian方法解非线性约束优化

A interior point scaling projected reduced Hessian method with combination of nonmonotonic backtracking technique and trust region strategy for nonlinear equality constrained optimization with nonegative constraint on variables is proposed. In order to deal with large problems,a pair of trust region subproblems in horizontal and vertical subspaces is used to replace the general full trust region subproblem. The horizontal trust region subproblem in the algorithm is only a general trust region subproblem while the vertical trust region subproblem is defined by a parameter size of the vertical direction subject only to an ellipsoidal constraint. Both trust region strategy and line search technique at each iteration switch to obtaining a backtracking step generated by the two trust region subproblems. By adopting the l1 penalty function as the merit function, the global convergence and fast local convergence rate of the proposed algorithm are established under some reasonable conditions. A nonmonotonic criterion and the second order correction step are used to overcome Maratos effect and speed up the convergence progress in some ill-conditioned cases.     

The block grade of a block Krylov space

The aim of the paper is to compile and compare basic theoretical facts on Krylov subspaces and block Krylov subspaces. Many Krylov (sub)space methods for solving a linear system Ax=b have the property that in exact computer arithmetic the true solution is found after ν iterations, where ν is the dimension of the largest Krylov subspace generated by A from r₀, the residual of the initial approximation x₀. This dimension is called the grade of r₀ with respect to A. Though the structure of block Krylov subspaces is more complicated than that of ordinary Krylov subspaces, we introduce here a block grade for which an analogous statement holds when block Krylov space methods are applied to linear systems with multiple, say s, right-hand sides. In this case, the s initial residuals are bundled into a matrix R₀ with s columns. The possibility of linear dependence among columns of the block Krylov matrix , which in practical algorithms calls for the deletion (or, deflation) of some columns, requires extra care. Relations between grade and block grade are also established, as well as relations to the corresponding notions of a minimal polynomial and its companion matrix.     

A review of formal orthogonality in Lanczos-based methods

Krylov subspace methods and their variants are presently the favorite iterative methods for solving a system of linear equations. Although it is a purely linear algebra problem, it can be tackled by the theory of formal orthogonal polynomials. This theory helps to understand the origin of the algorithms for the implementation of Krylov subspace methods and, moreover, the use of formal orthogonal polynomials brings a major simplification in the treatment of some numerical problems related to these algorithms. This paper reviews this approach in the case of Lanczos method and its variants, the novelty being the introduction of a preconditioner.     

Numerical solution of electromagnetic scattering from a large partly covered cavity

The paper focuses on the numerical study of electromagnetic scattering from two-dimensional (2D) large partly covered cavities, which is described by the Helmholtz equation with a nonlocal boundary condition on the aperture. The classical five-point finite difference method is applied for the discretization of the Helmholtz equation and a linear approximation is used for the nonlocal boundary condition. We prove the existence and uniqueness of the numerical solution when the medium in the cavity is y-direction layered or the number of the mesh points on the aperture is large enough. The fast algorithm proposed in Bao and Sun (2005) [2] for open cavity models is extended to solving the partly covered cavity problem with (vertically) layered media. A preconditioned Krylov subspace method is proposed to solve the partly covered cavity problem with a general medium, in which a layered medium model is used as a preconditioner of the general model. Numerical results for several types of partly covered cavities with different wave numbers are reported and compared with those by ILU-type preconditioning algorithms. Our numerical experiments show that the proposed preconditioning algorithm is more efficient for partly covered cavity problems, particularly with large wave numbers.     

Accuracy and effectiveness of the Lanczos algorithm for the symmetric eigenproblem

Eigenvalues and eigenvectors of a large sparse symmetric matrix A can be found accurately and often very quickly using the Lanczos algorithm without reorthogonalization. The algorithm gives essentially correct information on the eigensystem of A, although it does not necessarily give the correct multiplicity of multiple, or even single, eigenvalues. It is straightforward to determine a useful bound on the accuracy of every eigenvalue given by the algorithm. The initial behavior of the algorithm is surprisingly good: it produces vectors spanning the Krylov subspace of a matrix very close to A until this subspace contains an exact eigenvector of a matrix very close to A, and up to this point the effective behavior of the algorithm for the eigenproblem is very like that of the Lanczos algorithm using full reorthogonalization. This helps to explain the remarkable behavior of the basic Lanczos algorithm.     

Eigenvalue bounds for the Schur complement with a pressure convection–diffusion preconditioner in incompressible flow computations

If the stationary Navier–Stokes system or an implicit time discretization of the evolutionary Navier–Stokes system is linearized by a Picard iteration and discretized in space by a mixed finite element method, there arises a saddle point system which may be solved by a Krylov subspace method or an Uzawa type approach. For each of these resolution methods, it is necessary to precondition the Schur complement associated to the saddle point problem in question. In the work at hand, we give upper and lower bounds of the eigenvalues of this Schur complement under the assumption that it is preconditioned by a pressure convection–diffusion matrix.     

Multigrid preconditioned Krylov subspace methods for three-dimensional numerical solutions of the incompressible Navier–Stokes equations

We consider numerical methods for the incompressible Reynolds averaged Navier–Stokes equations discretized by finite difference techniques on non-staggered grids in body-fitted coordinates. A segregated approach is used to solve the pressure–velocity coupling problem. Several iterative pressure linear solvers including Krylov subspace and multigrid methods and their combination have been developed to compare the efficiency of each method and to design a robust solver. Three-dimensional numerical experiments carried out on scalar and vector machines and performed on different fluid flow problems show that a combination of multigrid and Krylov subspace methods is a robust and efficient pressure solver. This revised version was published online in June 2006 with corrections to the Cover Date.     

A finite-element capacitance matrix method for exterior Helmholtz problems

Summary. We introduce an algorithm for the efficient numerical solution of exterior boundary value problems for the Helmholtz equation. The problem is reformulated as an equivalent one on a bounded domain using an exact non-local boundary condition on a circular artificial boundary. An FFT-based fast Helmholtz solver is then derived for a finite-element discretization on an annular domain. The exterior problem for domains of general shape are treated using an imbedding or capacitance matrix method. The imbedding is achieved in such a way that the resulting capacitance matrix has a favorable spectral distribution leading to mesh independent convergence rates when Krylov subspace methods are used to solve the capacitance matrix equation. Received May 2, 1995