Multiplier convergence in trust-region methods with application to convergence of decomposition methods for MPECs

We study piecewise decomposition methods for mathematical programs with equilibrium constraints (MPECs) for which all constraint functions are linear. At each iteration of a decomposition method, one step of a nonlinear programming scheme is applied to one piece of the MPEC to obtain the next iterate. Our goal is to understand global convergence to B-stationary points of these methods when the embedded nonlinear programming solver is a trust-region scheme, and the selection of pieces is determined using multipliers generated by solving the trust-region subproblem. To this end we study global convergence of a linear trust-region scheme for linearly-constrained NLPs that we call a trust-search method. The trust-search has two features that are critical to global convergence of decomposition methods for MPECs: a robustness property with respect to switching pieces, and a multiplier convergence result that appears to be quite new for trust-region methods. These combine to clarify and strengthen global convergence of decomposition methods without resorting either to additional conditions such as eventual inactivity of the trust-region constraint, or more complex methods that require a separate subproblem for multiplier estimation.      

Block s‐step Krylov iterative methods

Block (including s‐step) iterative methods for (non)symmetric linear systems have been studied and implemented in the past. In this article we present a (combined) block s‐step Krylov iterative method for nonsymmetric linear systems. We then consider the problem of applying any block iterative method to solve a linear system with one right‐hand side using many linearly independent initial residual vectors. We present a new algorithm which combines the many solutions obtained (by any block iterative method) into a single solution to the linear system. This approach of using block methods in order to increase the parallelism of Krylov methods is very useful in parallel systems. We implemented the new method on a parallel computer and we ran tests to validate the accuracy and the performance of the proposed methods. It is expected that the block s‐step methods performance will scale well on other parallel systems because of their efficient use of memory hierarchies and their reduction of the number of global communication operations over the standard methods. Copyright © 2009 John Wiley & Sons, Ltd.     

Birkhoff normal form for splitting methods applied to semilinear Hamiltonian PDEs. Part II. Abstract splitting

We consider Hamiltonian PDEs that can be split into a linear unbounded operator and a regular non linear part. We consider abstract splitting methods associated with this decomposition where no discretization in space is made. We prove a normal form result for the corresponding discrete flow under generic non resonance conditions on the frequencies of the linear operator and on the step size, and under a condition of zero momentum on the nonlinearity. This result implies the conservation of the regularity of the numerical solution associated with the splitting method over arbitrary long time, provided the initial data is small enough. This result holds for rounded numerical schemes avoiding at each step possible high frequency energy drift. We apply these results to nonlinear Schrödinger equations as well as the nonlinear wave equation.     

Birkhoff normal form for splitting methods applied to semilinear Hamiltonian PDEs. Part I. Finite-dimensional discretization

We consider discretized Hamiltonian PDEs associated with a Hamiltonian function that can be split into a linear unbounded operator and a regular nonlinear part. We consider splitting methods associated with this decomposition. Using a finite-dimensional Birkhoff normal form result, we show the almost preservation of the actions of the numerical solution associated with the splitting method over arbitrary long time and for asymptotically large level of space approximation, provided the Sobolev norm of the initial data is small enough. This result holds under generic non-resonance conditions on the frequencies of the linear operator and on the step size. We apply these results to nonlinear Schrödinger equations as well as the nonlinear wave equation.     

Operator splitting methods for pricing American options under stochastic volatility

We consider the numerical pricing of American options under Heston’s stochastic volatility model. The price is given by a linear complementarity problem with a two-dimensional parabolic partial differential operator. We propose operator splitting methods for performing time stepping after a finite difference space discretization. The idea is to decouple the treatment of the early exercise constraint and the solution of the system of linear equations into separate fractional time steps. With this approach an efficient numerical method can be chosen for solving the system of linear equations in the first fractional step before making a simple update to satisfy the early exercise constraint. Our analysis suggests that the Crank–Nicolson method and the operator splitting method based on it have the same asymptotic order of accuracy. The numerical experiments show that the operator splitting methods have comparable discretization errors. They also demonstrate the efficiency of the operator splitting methods when a multigrid method is used for solving the systems of linear equations.     

Symmetric general linear methods

The article considers symmetric general linear methods, a class of numerical time integration methods which, like symmetric Runge–Kutta methods, are applicable to general time-reversible differential equations, not just those derived from separable second-order problems. A definition of time-reversal symmetry is formulated for general linear methods, and criteria are found for the methods to be free of linear parasitism. It is shown that symmetric parasitism-free methods cannot be explicit, but such a method of order 4 is constructed with only one implicit stage. Several characterizations of symmetry are given, and connections are made with G-symplecticity. Symmetric methods are shown to be of even order, a suitable symmetric starting method is constructed and shown to be essentially unique. The underlying one-step method is shown to be time-symmetric.     

Local convergence of quasi-Newton methods for B-differentiable equations

We study local convergence of quasi-Newton methods for solving systems of nonlinear equations defined by B-differentiable functions. We extend the classical linear and superlinear convergence results for general quasi-Newton methods as well as for Broyden's method. We also show how Broyden's method may be applied to nonlinear complementarity problems and illustrate its computational performance on two small examples.     

The inexact, inexact perturbed, and quasi-Newton methods are equivalent models

A classical model of Newton iterations which takes into account some error terms is given by the quasi-Newton method, which assumes perturbed Jacobians at each step. Its high convergence orders were characterized by Dennis and Moré [Math. Comp. 28 (1974), 549-560]. The inexact Newton method constitutes another such model, since it assumes that at each step the linear systems are only approximately solved; the high convergence orders of these iterations were characterized by Dembo, Eisenstat and Steihaug [SIAM J. Numer. Anal. 19 (1982), 400-408]. We have recently considered the inexact perturbed Newton method [J. Optim. Theory Appl. 108 (2001), 543-570] which assumes that at each step the linear systems are perturbed and then they are only approximately solved; we have characterized the high convergence orders of these iterates in terms of the perturbations and residuals.

In the present paper we show that these three models are in fact equivalent, in the sense that each one may be used to characterize the high convergence orders of the other two. We also study the relationship in the case of linear convergence and we deduce a new convergence result.

     

On superlinear convergence of quasi-Newton methods for nonsmooth equations

We show that strong differentiability at solutions is not necessary for superlinear convergence of quasi-Newton methods for solving nonsmooth equations. We improve the superlinear convergence result of Ip and Kyparisis for general quasi-Newton methods as well as the Broyden method. For a special example, the Newton method is divergent but the Broyden method is superlinearly convergent.     

Evaluating the Evans function: Order reduction in numerical methods

We consider the numerical evaluation of the Evans function, a Wronskian-like determinant that arises in the study of the stability of travelling waves. Constructing the Evans function involves matching the solutions of a linear ordinary differential equation depending on the spectral parameter. The problem becomes stiff as the spectral parameter grows. Consequently, the Gauss-Legendre method has previously been used for such problems; however more recently, methods based on the Magnus expansion have been proposed. Here we extensively examine the stiff regime for a general scalar Schrödinger operator. We show that although the fourth-order Magnus method suffers from order reduction, a fortunate cancellation when computing the Evans matching function means that fourth-order convergence in the end result is preserved. The Gauss-Legendre method does not suffer from order reduction, but it does not experience the cancellation either, and thus it has the same order of convergence in the end result. Finally we discuss the relative merits of both methods as spectral tools.

     

Partitioned quasi-Newton methods for sparse nonlinear equations

In this paper, we present two partitioned quasi-Newton methods for solving partially separable nonlinear equations. When the Jacobian is not available, we propose a partitioned Broyden’s rank one method and show that the full step partitioned Broyden’s rank one method is locally and superlinearly convergent. By using a well-defined derivative-free line search, we globalize the method and establish its global and superlinear convergence. In the case where the Jacobian is available, we propose a partitioned adjoint Broyden method and show its global and superlinear convergence. We also present some preliminary numerical results. The results show that the two partitioned quasi-Newton methods are effective and competitive for solving large-scale partially separable nonlinear equations.     

Convergence properties of trust region methods for linear and convex constraints

We develop a convergence theory for convex and linearly constrained trust region methods which only requires that the step between iterates produce a sufficient reduction in the trust region subproblem. Global convergence is established for general convex constraints while the local analysis is for linearly constrained problems. The main local result establishes that if the sequence converges to a nondegenerate stationary point then the active constraints at the solution are identified in a finite number of iterations. As a consequence of the identification properties, we develop rate of convergence results by assuming that the step is a truncated Newton method. Our development is mainly geometrical; this approach allows the development of a convergence theory without any linear independence assumptions.Work supported in part by the Applied Mathematical Sciences subprogram of the Office of Energy Research of the U.S. Department of Energy under Contract W-31-109-Eng-38.Work supported in part by the National Science Foundation grant DMS-8803206 and by the Air Force Office of Scientific Research grant AFSOR-860080.     

The convergence of inexact Chebyshev and Richardson iterative methods for solving linear systems

Summary The Chebyshev and second-order Richardson methods are classical iterative schemes for solving linear systems. We consider the convergence analysis of these methods when each step of the iteration is carried out inexactly. This has many applications, since a preconditioned iteration requires, at each step, the solution of a linear system which may be solved inexactly using an inner iteration. We derive an error bound which applies to the general nonsymmetric inexact Chebyshev iteration. We show how this simplifies slightly in the case of a symmetric or skew-symmetric iteration, and we consider both the cases of underestimating and overestimating the spectrum. We show that in the symmetric case, it is actually advantageous to underestimate the spectrum when the spectral radius and the degree of inexactness are both large. This is not true in the case of the skew-symmetric iteration. We show how similar results apply to the Richardson iteration. Finally, we describe numerical experiments which illustrate the results and suggest that the Chebyshev and Richardson methods, with reasonable parameter choices, may be more effective than the conjugate gradient method in the presence of inexactness.This work was supported in part by National Science Foundation Grants DCR-8412314 and DCR-8502014The work of this author was completed while he was on sabbatical leave at the Centre for Mathematical Analysis and Mathematical Sciences Research Institute at the Australian National University, Canberra, Australia     

Low rank methods for a class of generalized Lyapunov equations and related issues

In this paper, we study possible low rank solution methods for generalized Lyapunov equations arising in bilinear and stochastic control. We show that under certain assumptions one can expect a strong singular value decay in the solution matrix allowing for low rank approximations. Since the theoretical tools strongly make use of a connection to the standard linear Lyapunov equation, we can even extend the result to the $d$ -dimensional case described by a tensorized linear system of equations. We further provide some reasonable extensions of some of the most frequently used linear low rank solution techniques such as the alternating directions implicit (ADI) iteration and the Krylov-Plus-Inverted-Krylov (K-PIK) method. By means of some standard numerical examples used in the area of bilinear model order reduction, we will show the efficiency of the new methods.     

A class of pseudo Runge-Kutta methods

This paper develops a general theory for a class of Runge-Kutta methods which are based, in addition to the stages of the current step, also on the stages of the previous step. Such methods have been introduced previously for the case of one and two stages. We show that for any numbers of stages methods of orderp withs+1 p 2s can be constructed. The paper terminates with a study of step size change and stability.     

Bifurcations of hyperbolic fixed points for explicit Runge-Kutta methods

Discretization of autonomous ordinary differential equationsby numerical methods might, for certain step sizes, generatesolution sequences not corresponding to the underlying flow—so-called‘spurious solutions’ or ‘ghost solutions’.In this paper we explain this phenomenon for the case of explicitRunge-Kutta methods by application of bifurcation theory fordiscrete dynamical systems. An important tool in our analysisis the domain of absolute stability, resulting from the applicationof the method to a linear test problem. We show that hyperbolicfixed points of the (nonlinear) differential equation are inheritedby the difference scheme induced by the numerical method whilethe stability type of these inherited genuine fixed points iscompletely determined by the method's domain of absolute stability.We prove that, for small step sizes, the inherited fixed pointsexhibit the correct stability type, and we compute the correspondinglimit step size. Moreover, we show in which way the bifurcationsoccurring at the limit step size are connected to the valuesof the stability function on the boundary of the domain of absolutestability, where we pay special attention to bifurcations leadingto spurious solutions. In order to explain a certain kind ofspurious fixed points which are not connected to the set ofgenuine fixed points, we interprete the domain of absolute stabilityas a Mandeibrot set and generalize this approach to nonlinearproblems.     

The existence of symplectic general linear methods

We derive a criterion that any general linear method must satisfy if it is symplectic. It is shown, by considering the method over several steps, that the satisfaction of this condition leads to a reducibility in the method. Linking the symplectic criterion here to that for Runge–Kutta methods, we demonstrate that a general linear method is symplectic only if it can be reduced to a method with a single input value.      

Accelerating large partial EVD/SVD calculations by filtered block Davidson methods

Partial eigenvalue decomposition (PEVD) and partial singular value decomposition (PSVD) of large sparse matrices are of fundamental importance in a wide range of applications, including latent semantic indexing, spectral clustering, and kernel methods for machine learning. The more challenging problems are when a large number of eigenpairs or singular triplets need to be computed. We develop practical and efficient algorithms for these challenging problems. Our algorithms are based on a filter-accelerated block Davidson method. Two types of filters are utilized, one is Chebyshev polynomial filtering, the other is rational-function filtering by solving linear equations. The former utilizes the fastest growth of the Chebyshev polynomial among same degree polynomials; the latter employs the traditional idea of shift-invert, for which we address the important issue of automatic choice of shifts and propose a practical method for solving the shifted linear equations inside the block Davidson method. Our two filters can efficiently generate high-quality basis vectors to augment the projection subspace at each Davidson iteration step, which allows a restart scheme using an active projection subspace of small dimension. This makes our algorithms memory-economical, thus practical for large PEVD/PSVD calculations. We compare our algorithms with representative methods, including ARPACK, PROPACK, the randomized SVD method, and the limited memory SVD method. Extensive numerical tests on representative datasets demonstrate that, in general, our methods have similar or faster convergence speed in terms of CPU time, while requiring much lower memory comparing with other methods. The much lower memory requirement makes our methods more practical for large-scale PEVD/PSVD computations.     

Sufficient criteria are necessary for monotone control volume methods

Control volume methods are prevailing for solving the potential equation arising in porous media flow. The continuous form of this equation is known to satisfy a maximum principle, and it is desirable that the numerical approximation shares this quality. Recently, sufficient criteria were derived guaranteeing a discrete maximum principle for a class of control volume methods. We show that most of these criteria are also necessary. An implication of our work is that no linear nine-point control volume method can be constructed for quadrilateral grids in 2D that is exact for linear solutions while remaining monotone for general problems.     

On the convergence of alternating minimization methods in variational PGD

The approximation of solutions to partial differential equations by tensorial separated representations is one of the most efficient numerical treatment of high dimensional problems. The key step of such methods is the computation of an optimal low-rank tensor to enrich the obtained iterative tensorial approximation. In variational problems, this step can be carried out by alternating minimization (AM) technics, but the convergence of such methods presents a real challenge. In the present work, the convergence of rank-one AM algorithms for a class of variational linear elliptic equations is studied. More precisely, we show that rank-one AM-sequences are in general bounded in the ambient Hilbert tensor space and are compact if a uniform non-orthogonality condition between iterates and the reaction term is fulfilled. In particular, if a rank-one AM-sequence is weakly convergent then it converges strongly and the common limit is a solution of the rank-one optimization problem.