Past Computational Mathematics and Applications Seminar

Many recent algorithmic approaches involve the construction of a differential equation model for computational purposes, typically by introducing an artificial time variable. The actual computational model involves a discretization of the now time-dependent differential system, usually employing forward Euler. The resulting dynamics of such an algorithm is then a discrete dynamics, and it is expected to be ''close enough'' to the dynamics of the continuous system (which is typically easier to analyze) provided that small -- hence many -- time steps, or iterations, are taken. Indeed, recent papers in inverse problems and image processing routinely report results requiring thousands of iterations to converge. This makes one wonder if and how the computational modeling process can be improved to better reflect the actual properties sought.

In this talk we elaborate on several problem instances that illustrate the above observations. Algorithms may often lend themselves to a dual interpretation, in terms of a simply discretized differential equation with artificial time and in terms of a simple optimization algorithm; such a dual interpretation can be advantageous. We show how a broader computational modeling approach may possibly lead to algorithms with improved efficiency.

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.)

Spectral methods are a class of methods for solving PDEs numerically.

If the solution is analytic, it is known that these methods converge

exponentially quickly as a function of the number of terms used.

The basic spectral method only works in regular geometry (rectangles/disks).

A huge amount of effort has gone into extending it to

domains with a complicated geometry. Domain decomposition/spectral

element methods partition the domain into subdomains on which the PDE

can be solved (after transforming each subdomain into a

regular one). We take the dual approach - embedding the domain into

a larger regular domain - known as the fictitious domain method or

domain embedding. This method is extremely simple to implement and

the time complexity is almost the same as that for solving the PDE

on the larger regular domain. We demonstrate exponential convergence

for Dirichlet, Neumann and nonlinear problems. Time permitting, we

shall discuss extension of this technique to PDEs with discontinuous

coefficients.

In this talk we will mainly analyze the vibrations of a simplified 1-d model for a multi-body structure consisting of a finite number of flexible strings distributed along a planar graph. In particular we shall analyze how solutions propagate along the graph as time evolves. The problem of the observation of waves is a natural framework to analyze this issue. Roughly, the question can be formulated as follows: Can we obtain complete information on the vibrations by making measurements in one single extreme of the network? This formulation is relevant both in the context of control and inverse problems.

Using the Fourier development of solutions and techniques of Nonharmonic Fourier Analysis, we give spectral conditions that guarantee the observability property to hold in any time larger than twice the total lengths of the network in a suitable Hilbert that can be characterized in terms of Fourier series by means of properly chosen weights. When the network graph is a tree these weights can be identified.

Once this is done these results can be transferred to other models as the Schroedinger, heat or beam-type equations.

This lecture is based on results obtained in collaboration with Rene Dager.

The "secular equation" is a special way of expressing eigenvalue

problems in a variety of applications. We describe the secular

equation for several problems, viz eigenvector problems with a linear

constraint on the eigenvector and the solution of eigenvalue problems

where the given matrix has been modified by a rank one matrix. Next we

show how the secular equation can be approximated by use of the

Lanczos algorithm. Finally, we discuss numerical methods for solving

the approximate secular equation.

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.

Several formulations of structural optimization problems based on linear and nonlinear semidefinite programming will be presented. SDP allows us to formulate and solve problems with difficult constraints that could hardly be solved before. We will show that sometimes it is advantageous to prefer a nonlinear formulation to a linear one. All the presented formulations result in large-scale sparse (nonlinear) SDPs. In the second part of the talk we will show how these problems can be solved by our augmented Lagrangian code PENNON. Numerical examples will illustrate the talk.

Joint work with Michael Stingl.

We consider the parallel solution of sparse linear systems of equations in a limited memory environment. A preliminary out-of core version of a sparse multifrontal code called MUMPS (MUltifrontal Massively Parallel Solver) has been developed as part of a collaboration between CERFACS, ENSEEIHT and INRIA (ENS-Lyon and Bordeaux).

We first briefly describe the current status of the out-of-core factorization phase. We then assume that the factors have been written on the hard disk during the factorization phase and we discuss the design of an efficient solution phase.Two different approaches are presented to read data from the disk, with a discussion on the advantages and the drawbacks of each one.

Our work differs and extends the work of Rothberg and Schreiber (1999) and of Rotkin and Toledo (2004) because firstly we consider a parallel out-of-core context, and secondly we also study the performance of the solve phase.

This is work on collaboration with E. Agullo, I.S Duff, A. Guermouche, J.-Y. L'Excellent, T. Slavova

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.

Radial basis functions have been used for decades for the interpolation of scattered,

high-dimensional data. Recently they have attracted interest as methods for simulating

partial differential equations as well. RBFs do not require a grid or triangulation, they

offer the possibility of spectral accuracy with local refinement, and their implementation

is very straightforward. A number of theoretical and practical breakthroughs in recent years

has improved our understanding and application of these methods, and they are currently being

tested on real-world applications in shallow water flow on the sphere and tear film evolution

in the human eye.

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.

Many large-scale optimization problems arise in the context of the discretization of infinite dimensional applications. In such cases, the description of the finite-dimensional problem is not unique, but depends on the discretization used, resulting in a natural multi-level description. How can such a problem structure be exploited, in discretized problems or more generally? The talk will focus on discussing this issue in the context of unconstrained optimization and in relation with the classical multigrid approach to elliptic systems of partial differential equations. Both theoretical convergence properties of special purpose algorithms and their numerical performances will be discussed. Perspectives will also be given.

Collaboration with S. Gratton, A. Sartenaer and M. Weber.

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.

The aim of this talk is to describe several methods for numerically approximating

the integral of a multivariate highly oscillatory function. We begin with a review

of the asymptotic and Filon-type methods developed by Iserles and Nørsett. Using a

method developed by Levin as a point of departure we will construct a new method that

uses the same information as the Filon-type method, and obtains the same asymptotic

order, while not requiring moments. This allows us to integrate over nonsimplicial

domains, and with complicated oscillators.

High-performance computing is an important tool for computational science.

Oxford University has recently decided to invest £3M in a new supercomputing

facility which is under development now. In this seminar I will give an overview

of some activities in Oxford and provide a vision for supercomputing here.

I will discuss some of the numerical analysis software and tools,

such as Distributed Matlab and indicate some of the challenges at

the intersection of numerical analysis and high-performance computing.

The aim of this talk is to give some understanding of the theory of matrix $p$'th roots (solutions to the nonlinear matrix equation $X^{p} = A$), to explain how and how not to compute roots, and to describe some applications. In particular, an application in finance will be described concerning roots of transition matrices from Markov models.