Rutherford Appleton Laboratory, nr Didcot

We consider the solution of a linear system of equations using the GMRES iterative method. In some applications, performing inexact matrix-vector products in this method may be interesting, provided that a reasonable convergence of GMRES is achieved. A GMRES algorithm where the matrix vector product is performed inexactly is termed ”inexact GMRES algorithm”. An application of this idea occurs in computational electromagnetics, where the fast multipole method provides approximations of the matrix-vector product within a user-defined precision, and where these inaccurate matrix-vector products are all the more cheaper (in terms of CPU time) as the user-defined precision is low. The key point is then to design a control strategy of the accuracy of the matrix-vector product so that the GMRES converges better in the sense that 1) the inexact method achieves a satisfactory limiting accuracy, 2) within a reasonable number of steps.

\\

In [1], a relaxation strategy is proposed for general systems and validated on a large set of numerical experiments. This work is based on heuristic considerations and proposes a strategy that enables a convergence of the GMRES iterates $x_{k}$ within a relative normwise backward error $\frac{\|b−Ax_{k}\|}{\|A\| \|x_{k}\| + \|b\|}$ less than a prescribed quantity $\eta$ > 0, on a significant number of numerical experiments. Similar strategies have been applied to the solution of device simulation problems using domain decomposition [2] and to the preconditioning of a radiation diffusion problem in [5].

\\

A step toward a theoretical explanation of the observed behaviour of the inexact GMRES is proposed in [3, 4]. In this talk, we show that in spite of this considerable theoretical study, the experimental work of [1] is not fully understood yet. We give an overview of the questions that still remains open both in exact arithmetic and in floating-point arithmetic, and we provide some insights into the solution of some of them. Possible applications of this work for the preconditioned GMRES method, when the matrix-vector product is accurate but the preconditioning operation is approximate, are also investigated, based on [3].

\\

\\

References

\\

[1] A. Bouras and V. Frayss´e. Inexact matrix-vector products in Krylov methods for solving linear systems: a relaxation strategy. SIAM Journal on Matrix Analysis and Applications, 2004. To appear.

\\

[2] A. Bouras, V. Frayss´e, and L. Giraud. A relaxation strategy for inner-outer linear solvers in domain decomposition methods. Technical Report TR/PA/00/17, CERFACS, Toulouse, France, 2000.

\\

[3] V. Simoncini and D. B. Szyld. Theory of inexact Krylov subspace methods and applications to scientific computing. SIAM Journal Scientific Computing, 25:454–477, 2003.

\\

[4] J. van den Eshof and G. L. G. Sleijpen. Inexact Krylov subspace methods for linear systems. SIAM Journal on Matrix Analysis and Applications, February 2004. To appear.

\\

[5] J. S. Warsa, M. Benzi, T. A. Warein, and J. E. Morel. Preconditioning a mixed discontinuous finite element method for radiation diffusion. Numerical Linear Algebra with Applications, 2004. To appear.

Consider approximating a set of discretely defined values $f_{1}, \ldots , f_{m}$ say at $x=x_{1}, x_{2}, \ldots, x_{m}$, with a chosen approximating form. Given prior knowledge that noise is present and that some might be outliers, a standard least squares approach based on $l_{2}$ norm of the error $\epsilon$ may well provide poor estimates. We instead consider a least squares approach based on a modified measure of the form $\tilde{\epsilon} = \epsilon (1+c^{2}\epsilon^{2})^{-\frac{1}{2}}$, where $c$ is a constant to be fixed.

\\

The choice of the constant $c$ in this estimator has a significant effect on the performance of the estimator both in terms of its algorithmic convergence to a solution and its ability to cope effectively with outliers. Given a prior estimate of the likely standard deviation of the noise in the data, we wish to determine a value of $c$ such that the estimator behaves like a robust estimator when outliers are present but like a least squares estimator otherwise.

\\

We describe approaches to determining suitable values of $c$ and illustrate their effectiveness on approximation with polynomial and radial basis functions. We also describe algorithms for computing the estimates based on an iteratively weighted linear least squares scheme.

The use of weighted matchings is becoming increasingly standard in the

solution of sparse linear systems. While non-symmetric permutations based on these

matchings have been the state-of-the-art for

several years (especially for direct solvers), approaches for symmetric

matrices have only recently gained attention.

\\

In this talk we discuss results of our work on using weighted matchings in

the preconditioning of symmetric indefinite linear systems, following ideas

introduced by Duff and Gilbert. In order to maintain symmetry,

the weighted matching is symmetrized and the cycle structure of the

resulting matching is used to build reorderings that form small diagonal

blocks from the matched entries.

\\

For the preconditioning we investigated two approaches. One is an

incomplete $LDL^{T}$ preconditioning, that chooses 1x1 or 2x2 diagonal pivots

based on a simple tridiagonal pivoting criterion. The second approach

targets distributed computing, and is based on factorized sparse approximate

inverses, whose existence, in turn, is based on the existence of an $LDL^{T}$

factorization. Results for a number of comprehensive test sets are given,

including comparisons with sparse direct solvers and other preconditioning

approaches.

An interior-point method for solving mathematical programs with

equilibrium constraints (MPECs) is proposed. At each iteration of the

algorithm, a single primal-dual step is computed from each subproblem of

a sequence. Each subproblem is defined as a relaxation of the MPEC with

a nonempty strictly feasible region. In contrast to previous

approaches, the proposed relaxation scheme preserves the nonempty strict

feasibility of each subproblem even in the limit. Local and superlinear

convergence of the algorithm is proved even with a less restrictive

strict complementarity condition than the standard one. Moreover,

mechanisms for inducing global convergence in practice are proposed.

Numerical results on the MacMPEC test problem set demonstrate the

fast-local convergence properties of the algorithm.

The development and quantitative validation of complex nonlinear differential equation models is a difficult task that requires the support by numerical methods for sensitivity analysis, parameter estimation, and the optimal design of experiments. The talk first presents particularly efficient "simultaneous" boundary value problems methods for parameter estimation in nonlinear differential algebraic equations, which are based on constrained Gauss-Newton-type methods and a time domain decomposition by multiple shooting. They include a numerical analysis of the well-posedness of the problem and an assessment of the error of the resulting parameter estimates. Based on these approaches, efficient optimal control methods for the determination of one, or several complementary, optimal experiments are developed, which maximize the information gain subject to constraints such as experimental costs and feasibility, the range of model validity, or further technical constraints.

Special emphasis is placed on issues of robustness, i.e. how to reduce the sensitivity of the problem solutions with respect to uncertainties - such as outliers in the measurements for parameter estimation, and in particular the dependence of optimum experimental designs on the largely unknown values of the model parameters. New numerical methods will be presented, and applications will be discussed that arise in satellite orbit determination, chemical reaction kinetics, enzyme kinetics and robotics. They indicate a wide scope of applicability of the methods, and an enormous potential for reducing the experimental effort and improving the statistical quality of the models.

(Based on joint work with H. G. Bock, S. Koerkel, and J. P. Schloeder.)

The numerical study of lattice quantum chromodynamics (QCD) is an attempt to extract predictions about the world around us from the standard model of physics. Worldwide, there are several large collaborations on lattice QCD methods, with terascale computing power dedicated to these problems. Central to the computation in lattice QCD is the inversion of a series of fermion matrices, representing the interaction of quarks on a four-dimensional space-time lattice. In practical computation, this inversion may be approximated based on the solution of a set of linear systems.

In this talk, I will present a basic description of the linear algebra problems in lattice QCD and why we believe that multigrid methods are well-suited to effectively solving them. While multigrid methods are known to be efficient solution techniques for many operators, those arising in lattice QCD offer new challenges, not easily handled by classical multigrid and algebraic multigrid approaches. The role of adaptive multigrid techniques in addressing the fermion matrices will be highlighted, along with preliminary results for several model problems.

In this talk we shall be looking at recent and forthcoming developments in the widely used LAPACK and ScaLAPACK numerical linear algebra libraries.

Improvements include the following: Faster algorithms, better numerical methods, memory hierarchy optimizations, parallelism, and automatic performance tuning to accommodate new architectures; more accurate algorithms, and the use of extra precision; expanded functionality, including updating and downdating and new eigenproblems; putting more of LAPACK into ScaLAPACK; and improved ease of use with friendlier interfaces in multiple languages. To accomplish these goals we are also relying on better software engineering techniques and contributions from collaborators at many institutions.

After an overview, this talk will highlight new more accurate algorithms; faster algorithms, including those for pivoted Cholesky and updating of factorizations; and hybrid data formats.

This is joint work with Jim Demmel, Jack Dongarra and the LAPACK/ScaLAPACK team.

A strict adherence to threshold pivoting in the direct solution of symmetric indefinite problems can result in substantially more work and storage than forecast by an sparse analysis of the symmetric problem. One way of avoiding this is to use static pivoting where the data structures and pivoting sequence generated by the analysis are respected and pivots that would otherwise be very small are replaced by a user defined quantity. This can give a stable factorization but of a perturbed matrix.

The conventional way of solving the sparse linear system is then to use iterative refinement (IR) but there are cases where this fails to converge. We will discuss the use of more robust iterative methods, namely GMRES and its variant FGMRES and their backward stability when the preconditioning is performed by HSL_M57 with a static pivot option.

Several examples under Matlab will be presented.

Joint work with Yogi Erlangga and Kees Vuik.

Shifted Laplace preconditioners have attracted considerable attention as a technique to speed up convergence of iterative solution methods for the Helmholtz equation. In this paper we present a comprehensive spectral analysis of the Helmholtz operator preconditioned with a shifted Laplacian. Our analysis is valid under general conditions. The propagating medium can be heterogeneous, and the analysis also holds for different types of damping, including a radiation condition for the boundary of the computational domain. By combining the results of the spectral analysis of the preconditioned Helmholtz operator with an upper bound on the GMRES-residual norm we are able to provide an optimal value for the shift, and to explain the mesh-depency of the convergence of GMRES preconditioned with a shifted Laplacian. We illustrate our results with a seismic test problem.

The talk starts with a general introduction of the convex

quadratic semidefinite programming problem (QSDP), followed by a number of

interesting examples arising from finance, statistics and computer sciences.

We then discuss the concept of primal nondegeneracy for QSDP and show that

some QSDPs are nondegenerate and others are not even under the linear

independence assumption. The talk then moves on to the algorithmic side by

introducing the dual approach and how it naturally leads to Newton's method,

which is quadratically convergent for degenerate problems. On the

implementation side of the Newton method, we stress that direct methods for

the linear equations in Newton's method are impossible simply because the

equations are quite large in size and dense. Our numerical experiments use

the conjugate gradient method, which works quite well for the nearest

correlation matrix problem. We also discuss difficulties for us to find

appropriate preconditioners for the linear system encountered. The talk

concludes in discussing some other available methods and some future topics.