Past

Image processing is an area with many important applications, as well as challenging problems for mathematicians. In particular, Fourier/wavelets analysis and stochastic/statistical methods have had major impact in this area. Recently, there has been increased interest in a new and complementary approach, using partial differential equations (PDEs) and differential-geometric models. It offers a more systematic treatment of geometric features of mages, such as shapes, contours and curvatures, etc., as well as allowing the wealth of techniques developed for PDEs and Computational Fluid Dynamics (CFD) to be brought to bear on image processing tasks.

I'll use two examples from my recent work to illustrate this synergy:

1. A unified image restoration model using Total Variation (TV) which can be used to model denoising, deblurring, as well as image inpainting (e.g. restoring old scratched photos). The TV idea can be traced to shock capturing methods in CFD and was first used in image processing by Rudin, Osher and Fatemi.

2. An "active contour" model which uses a variational level set method for object detection in scalar and vector-valued images. It can detect objects not necessarily defined by sharp edges, as well as objects undetectable in each channel of a vector-valued image or in the combined intensity. The contour can go through topological changes, and the model is robust to noise. The level set method was originally developed by Osher and Sethian for tracking interfaces in CFD.

(The above are joint works with Jackie Shen at the Univ. of Minnesota and Luminita Vese in the Math Dept at UCLA.)

2.00 pm Professor Iain Duff (RAL) Opening remarks

2.15 pm Professor M J D Powell (University of Cambridge)

Some developments of work with Alan on cubic splines

3.00 pm Professor Kevin Burrage (University of Queensland)

Stochastic models and simulations for chemically reacting systems

3.30 pm Tea/Coffee

4.00 pm Professor John Reid (RAL)

Sparse matrix research at Harwell and the Rutherford Appleton Laboratory

4.30 pm Dr Ian Jones (AEA PLC)

Computational fluid dynamics and the role of stiff solvers

5.00 pm Dr Lawrence Daniels (Hyprotech UK Ltd)

Current work with Alan on ODE solvers for HSL

Long-step primal-dual path-following algorithms constitute the

framework of practical interior point methods for

solving linear programming problems. We consider

such an algorithm and a second order variant of it.

We address the problem of the convergence of

the sequences of iterates generated by the two algorithms

to the analytic centre of the optimal primal-dual set.

Several real Lie and Jordan algebras, along with their associated

automorphism groups, can be elegantly expressed in the quaternion tensor

algebra. The resulting insight into structured matrices leads to a class

of simple Jacobi algorithms for the corresponding $n \times n$ structured

eigenproblems. These algorithms have many desirable properties, including

parallelizability, ease of implementation, and strong stability.

A method for computing periodic orbits for the Navier-Stokes

equations will be presented. The method uses a finite-element Galerkin

discretisation for the spatial part of the problem and a spectral

Galerkin method for the temporal part of the problem. The method will

be illustrated by calculations of the periodic flow behind a circular

cylinder in a channel. The problem has a simple reflectional symmetry

and it will be explained how this can be exploited to reduce the cost

of the computations.

The convergence of Krylov subspace methods like conjugate gradients

depends on the eigenvalues of the underlying matrix. In many cases

the exact location of the eigenvalues is unknown, but one has some

information about the distribution of eigenvalues in an asymptotic

sense. This could be the case for linear systems arising from a

discretization of a PDE. The asymptotic behavior then takes place

when the meshsize tends to zero.

\\

\\

We discuss two possible approaches to study the convergence of

conjugate gradients based on such information.

The first approach is based on a straightforward idea to estimate

the condition number. This method is illustrated by means of a

comparison of preconditioning techniques.

The second approach takes into account the full asymptotic

spectrum. It gives a bound on the asymptotic convergence factor

which explains the superlinear convergence observed in many situations.

This method is mathematically more involved since it deals with

potential theory. I will explain the basic ideas.

We develop an algorithm for estimating the local Sobolev regularity index

of a given function by monitoring the decay rate of its Legendre expansion

coefficients. On the basis of these local regularities, we design and

implement an hp--adaptive finite element method based on employing

discontinuous piecewise polynomials, for the approximation of nonlinear

systems of hyperbolic conservation laws. The performance of the proposed

adaptive strategy is demonstrated numerically.

The study of the finite precision behaviour of numerical algorithms dates back at least as far as Turing and Wilkinson in the 1940s. At the start of the 21st century, this area of research is still very active.

We focus on some topics of current interest, describing recent developments and trends and pointing out future research directions. The talk will be accessible to those who are not specialists in numerical analysis.

Specific topics intended to be addressed include

• Floating point arithmetic: correctly rounded elementary functions, and the fused multiply-add operation.
• The use of extra precision for key parts of a computation: iterative refinement in fixed and mixed precision.
• Gaussian elimination with rook pivoting and new error bounds for Gaussian elimination.
• Automatic error analysis.
• Application and analysis of hyperbolic transformations.

We describe "avoidance of crossing" and its explanation by von

Neumann and Wigner. We show Lax's criterion for degeneracy and then

discover matrices whose determinants give the discriminant of the

given matrix. This yields a simple proof of the bound given by

Ilyushechkin on the number of terms in the expansion of the discriminant

as a sum of squares. We discuss the 3 x 3 case in detail.

The talk will discuss unsymmetric sparse LU factorization based on

the Markowitz pivot selection criterium. The key question for the

author is the following: Is it possible to implement a sparse

factorization where the overhead is limited to a constant times

the actual numerical work? In other words, can the work be bounded

by o(sum(k, M(k)), where M(k) is the Markowitz count in pivot k.

The answer is probably NO, but how close can we get? We will give

several bad examples for traditional methods and suggest alternative

methods / data structure both for pivot selection and for the sparse

update operations.

We discuss two filters that are frequently used to smooth data.

One is the (nonlinear) median filter, that chooses the median

of the sample values in the sliding window. This deals effectively

with "outliers" that are beyond the correct sample range, and will

never be chosen as the median. A straightforward implementation of

the filter is expensive for large windows, particularly in two dimensions

(for images).

\\

\\

The second filter is linear, and known as "Savitzky-Golay". It is

frequently used in spectroscopy, to locate positions and peaks and

widths of spectral lines. This filter is based on a least-squares fit

of the samples in the sliding window to a polynomial of relatively

low degree. The filter coefficients are unlike the equiripple filter

that is optimal in the maximum norm, and the "maxflat" filters that

are central in wavelet constructions. Should they be better known....?

\\

\\

We will discuss the analysis and the implementation of both filters.

The asymptotic rate of convergence of the trapezium rule is

defined, for smooth functions, by the Euler-Maclaurin expansion.

The extension to other methods, such as Gauss rules, is straightforward;

this talk begins with some special cases, such as Periodic functions, and

functions with various singularities.

\\

\\

Convergence rates for ODEs (Initial and Boundary value problems)

and for PDEs are available, but not so well known. Extension to singular

problems seems to require methods specific to each situation. Some of

the results are unexpected - to me, anyway.

In the past few years we have developed some expertise in solving optimization

problems that involve large scale simulations in various areas of Computational

Geophysics and Engineering. We will discuss some of those applications here,

namely: inversion of seismic data, characterization of piezoelectrical crystals

material properties, optimal design of piezoelectrical transducers and

opto-electronic devices, and the optimal design of steel structures.

\\

\\

A common theme among these different applications is that the goal functional

is very expensive to evaluate, often, no derivatives are readily available, and

some times the dimensionality can be large.

\\

\\

Thus parallelism is a need, and when no derivatives are present, search type

methods have to be used for the optimization part. Additional difficulties can

be ill-conditioning and non-convexity, that leads to issues of global

optimization. Another area that has not been extensively explored in numerical

optimization and that is important in real applications is that of

multiobjective optimization.

\\

\\

As a result of these varied experiences we are currently designing a toolbox

to facilitate the rapid deployment of these techniques to other areas of

application with a minimum of retooling.

A-Posteriori Error estimates for high order Godunov finite

volume methods are presented which exploit the two solution

representations inherent in the method, viz. as piecewise

constants $u_0$ and cell-wise $q$-th order reconstructed

functions $R^0_q u_0$. The analysis provided here applies

directly to schemes such as MUSCL, TVD, UNO, ENO, WENO or any

other scheme that is a faithful extension of Godunov's method

to high order accuracy in a sense that will be made precise.

Using standard duality arguments, we construct exact error

representation formulas for derived functionals that are

tailored to the class of high order Godunov finite volume

methods with data reconstruction, $R^0_q u_0$. We then consider

computable error estimates that exploit the structure of higher

order Godunov finite volume methods. The analysis technique used

in this work exploits a certain relationship between higher

order Godunov methods and the discontinuous Galerkin method.

Issues such as the treatment of nonlinearity and the optional

post-processing of numerical dual data are also discussed.

Numerical results for linear and nonlinear scalar conservation

laws are presented to verify the analysis. Complete details can

be found in a paper appearing in the proceedings of FVCA3,

Porquerolles, France, June 24-28, 2002.

SMP (Symmetric Multi-Processors) hardware technologies are very popular

with vendors and end-users alike for a number of reasons. However, true

shared memory parallelism has experienced somewhat slower to take up

amongst the scientific-programming community. NAG has been at the

forefront of SMP technology for a number of years, and the NAG SMP

Library has shown the potential of SMP systems.

\\

\\

At the very high end, SMP hardware technologies are used as building

blocks of modern supercomputers, which truly consist of clusters of SMP

systems, for which no dedicated model of parallelism yet exists.

\\

\\

The aim of this talk is to introduce SMP systems and their potential.

Results from our work at NAG will also be introduced to show how SMP

parallelism, based on a shared memory paradigm, can be used to very

good effect and can produce high performance, scalable software. The

talk also aims to discuss some aspects of the apparent slow take up of

shared memory parallelism and the potential competition from PC (i.e.

Intel)-based cluster technology. The talk then aims to explore the

potential of SMP technology within "hybrid parallelism", i.e. mixed

distributed and shared memory modes, illustrating the point with some

preliminary work carried out by the author and others. Finally, a

number of potential future challenges to numerical analysts will be

discussed.

\\

\\

The talk is aimed at all who are interested in SMP technologies for

numerical computing, irrespective of any previous experience in the

field. The talk aims to stimulate discussion, by presenting some ideas,

backing these with data, not to stifle it in an ocean of detail!

It is well known that discrete solutions to the convection-diffusion

equation contain nonphysical oscillations when boundary layers are present

but not resolved by the discretisation. For the Galerkin finite element

method with linear elements on a uniform 1D grid, a precise statement as

to exactly when such oscillations occur can be made, namely, that for a

problem with mesh size h, constant advective velocity and different values

at the left and right boundaries, oscillations will occur if the mesh

P\'{e}clet number $P_e$ is greater than one. In 2D, however, the situation

is not so well understood. In this talk, we present an analysis of a 2D

model problem on a square domain with grid-aligned flow which enables us

to clarify precisely when oscillations occur, and what can be done to

prevent them. We prove the somewhat surprising result that there are

oscillations in the 2D problem even when $P_e$ is less than 1. Also, we show that there

are distinct effects arising from differences in the top and bottom

boundary conditions (equivalent to those seen in 1D), and the non-zero

boundaries parallel to the flow direction.

Algebra based modeling systems are becoming essential elements in the

application of large and complex mathematical programs. These systems

enable the abstraction, expression and translation of practical

problems into reliable and effective operational systems. They provide

the bridged between algorithms and real world problems by automating

the problem analysis and translation into specific data structures and

provide computational services required by different solvers. The

modeling system GAMS will be used to illustrate the design goals and

main features of such systems. Applications in use and under

development will be used to provide the context for discussing the

changes in user focus and future requirements. This presents new sets

of opportunities and challenges to the supplier and implementer of

mathematical programming solvers and modeling systems.

A systematic approach to error control and mesh adaptation for

optimal control of systems governed by PDEs is presented.

Starting from a coarse mesh, the finite element spaces are successively

enriched in order to construct suitable discrete models.

This process is guided by an a posteriori error estimator which employs

sensitivity factors from the adjoint equation.

We consider different examples with the stationary Navier-Stokes

equations as state equation.

In this presentation we examine the convergence characteristics of a

Krylov subspace solver preconditioned by a new indefinite

constraint-type preconditioner, when applied to discrete systems

arising from low-order mixed finite element approximation of the

classical biharmonic problem. The preconditioning operator leads to

preconditioned systems having an eigenvalue distribution consisting of

a tightly clustered set together with a small number of outliers. We

compare the convergence characteristics of a new approach with the

convergence characteristics of a standard block-diagonal Schur

complement preconditioner that has proved to be extremely effective in

the context of mixed approximation methods.

\\

\\

In the second part of the presentation we are concerned with the

efficient parallel implementation of proposed algorithm on modern

shared memory architectures. We consider use of the efficient parallel

"black-box'' solvers for the Dirichlet Laplacian problems based on

sparse Cholesky factorisation and multigrid, and for this purpose we

use publicly available codes from the HSL library and MGNet collection.

We compare the performance of our algorithm with sparse direct solvers

from the HSL library and discuss some implementation related issues.

(Joint work with Felipe Cucker and Dennis Cheung, City University of Hong Kong.)

\\

\\

Condition numbers are important complexity-theoretic tools to capture

a "distillation" of the input aspects of a computational problem that

determine the running time of algorithms for its solution and the

sensitivity of the computed output. The motivation for our work is the

desire to understand the average case behaviour of linear programming

algorithms for a large class of randomly generated input data in the

computational model of a machine that computes with real numbers. In

this model it is not known whether linear programming is polynomial

time solvable, or so-called "strongly polynomial". Closely related to

linear programming is the problem of either proving non-existence of

or finding an explicit example of a point in a polyhedral cone defined

in terms of certain input data. A natural condition number for this

computational problem was developed by Cheung and Cucker, and we analyse

its distributions under a rather general family of input distributions.

We distinguish random sampling of primal and dual constraints

respectively, two cases that necessitate completely different techniques

of analysis. We derive the exact exponents of the decay rates of the

distribution tails and prove various limit theorems of complexity

theoretic importance. An interesting result is that the existence of

the k-th moment of Cheung-Cucker's condition number depends only very

mildly on the distribution of the input data. Our results also form

the basis for a second paper in which we analyse the distributions of

Renegar's condition number for the randomly generated linear programming

problem.

Toeplitz matrices enjoy the dual virtues of ubiquity and beauty. We begin this talk by surveying some of the interesting spectral properties of such matrices, emphasizing the distinctions between infinite-dimensional Toeplitz matrices and the large-dimensional limit of the corresponding finite matrices. These basic results utilize the algebraic Toeplitz structure, but in many applications, one is forced to spoil this structure with some perturbations (e.g., by imposing boundary conditions upon a finite difference discretization of an initial-boundary value problem). How do such

perturbations affect the eigenvalues? This talk will address this question for "localized" perturbations, by which we mean perturbations that are restricted to a single entry, or a block of entries whose size remains fixed as the matrix dimension grows. One finds, for a broad class of matrices, that sufficiently small perturbations fail to alter the spectrum, though the spectrum is exponentially sensitive to other perturbations. For larger real single-entry

perturbations, one observes the perturbed eigenvalues trace out curves in the complex plane. We'll show a number of illustrations of this phenomenon for tridiagonal Toeplitz matrices.

\\

\\

This talk describes collaborative work with Albrecht Boettcher, Marko Lindner, and Viatcheslav Sokolov of TU Chemnitz.

We investigate the use of interior-point and semismooth methods for solving

quadratic programming problems with a small number of linear constraints,

where the quadratic term consists of a low-rank update to a positive

semi-definite matrix. Several formulations of the support vector machine

fit into this category. An interesting feature of these particular problems

is the volume of data, which can lead to quadratic programs with between 10

and 100 million variables and, if written explicitly, a dense $Q$ matrix.

Our codes are based on OOQP, an object-oriented interior-point code, with the

linear algebra specialized for the support vector machine application.

For the targeted massive problems, all of the data is stored out of core and

we overlap computation and I/O to reduce overhead. Results are reported for

several linear support vector machine formulations demonstrating that the

methods are reliable and scalable and comparing the two approaches.

Non-selfadjoint singular differential equation eigenproblems arise in a number of contexts, including scattering theory, the study of quantum-mechanical resonances, and hydrodynamic and magnetohydrodynamic stability theory.

\\

\\

It is well known that the spectra of non-selfadjoint operators can be pathologically sensitive to perturbation of the

operator. Wilkinson provides matrix examples in his classical text, while Trefethen has studied the phenomenon extensively through pseudospectra, which he argues are often of more physical relevance than the spectrum itself. E.B. Davies has studied the phenomenon particularly in the context of Sturm-Liouville operators and has shown that the eigenfunctions and associated functions of non-selfadjoint singular Sturm-Liouville operators may not even form a complete set in $L^2$.

\\

\\

In this work we ask the question: under what conditions can one expect the regularization process used for selfadjoint singular Sturm-Liouville operators to be successful for non-selfadjoint operators? The answer turns out to depend in part on the so-called Sims Classification of the problem. For Sims Case I the process is not guaranteed to work, and indeed Davies has very recently described the way in which spurious eigenvalues may be generated and converge to certain curves in the complex plane.

\\

\\

Using the Titchmarsh-Weyl theory we develop a very simple numerical procedure which can be used a-posteriori to distinguish genuine eigenvalues from spurious ones. Numerical results indicate that it is able to detect not only the spurious eigenvalues due to the regularization process, but also spurious eigenvalues due to the numerics on an already-regular problem. We present applications to quantum mechanical resonance calculations and to the Orr-Sommerfeld equation.

\\

\\

This work, in collaboration with B.M. Brown in Cardiff, has recently been generalized to Hamiltonian systems.

Stiff systems of ODEs arise commonly when solving PDEs by spectral methods,

so conventional explicit time-stepping methods require very small time steps.

The stiffness arises predominantly through the linear terms, and these

terms can be handled implicitly or exactly, permitting larger time steps.

This work develops and investigates a class of methods known as

'exponential time differencing'. These methods are shown to have a

number of advantages over the more well-known linearly implicit

methods and integrating factor methods.

The Kestrel interface for submitting optimization problems to the NEOS Server augments the established e-mail, socket, and web interfaces by enabling easy usage of remote solvers from a local modeling environment.

\\

\\

Problem generation, including the run-time detection of syntax errors, occurs on the local machine using any available modeling language facilities. Finding a solution to the problem takes place on a remote machine, with the result returned in the native modeling language format for further processing. A byproduct of the Kestrel interface is the ability to solve multiple problems generated by a modeling language in parallel.

\\

\\

This mechanism is used, for example, in the GAMS/AMPL solver available through the NEOS Server, which internally translates a submitted GAMS problem into AMPL. The resulting AMPL problem is then solved through the NEOS Server via the Kestrel interface. An advantage of this design is that the GAMS to AMPL translator does not need to be collocated with the AMPL solver used, removing restrictions on solver choice and reducing administrative costs.

\\

\\

This talk is joint work with Elizabeth Dolan.

Tridiagonal matrices and three term recurrences and second order equations appear amazingly often, throughout all of mathematics. We won't try to review this subject. Instead we look in two less familiar directions.

\\

\\

Here is a tridiagonal matrix problem that waited surprisingly long for a solution. Forward elimination factors T into LDU, with the pivots in D as usual. Backward elimination, from row n to row 1, factors T into U_D_L_. Parlett asked for a proof that diag(D + D_) = diag(T) + diag(T^-1).^-1. In an excellent paper (Lin Alg Appl 1997) Dhillon and Parlett extended this four-diagonal identity to block tridiagonal matrices, and also applied it to their "Holy Grail" algorithm for the eigenproblem. I would like to make a different connection, to the Kalman filter.

\\

\\

The second topic is a generalization of tridiagonal to "tree-diagonal". Unlike the interval, the tree can branch. In the matrix T, each vertex is connected only to its neighbors (but a branch point has more than two neighbors). The continuous analogue is a second order differential equation on a tree. The "non-jump" conditions at a meeting of N edges are continuity of the potential (N-1 equations) and Kirchhoff's Current Law (1 equation). Several important properties of tridiagonal matrices, including O(N) algorithms, survive on trees.

No seminar this week

Let the thin plate spline radial basis function method be applied to

interpolate values of a smooth function $f(x)$, $x \!\in\! {\cal R}^d$.

It is known that, if the data are the values $f(jh)$, $j \in {\cal Z}^d$,

where $h$ is the spacing between data points and ${\cal Z}^d$ is the

set of points in $d$ dimensions with integer coordinates, then the

accuracy of the interpolant is of magnitude $h^{d+2}$. This beautiful

result, due to Buhmann, will be explained briefly. We will also survey

some recent findings of Bejancu on Lagrange functions in two dimensions

when interpolating at the integer points of the half-plane ${\cal Z}^2

\cap \{ x : x_2 \!\geq\! 0 \}$. Most of our attention, however, will

be given to the current research of the author on interpolation in one

dimension at the points $h {\cal Z} \cap [0,1]$, the purpose of the work

being to establish theoretically the apparent deterioration in accuracy

at the ends of the range from ${\cal O} ( h^3 )$ to ${\cal O} ( h^{3/2}

)$ that has been observed in practice. The analysis includes a study of

the Lagrange functions of the semi-infinite grid ${\cal Z} \cap \{ x :

x \!\geq\! 0 \}$ in one dimension.

In this talk we review experiences of using the Harwell Subroutine

Library and other numerical software codes in implementing large scale

solvers for commercial industrial process simulation packages. Such

packages are required to solve problems in an efficient and robust

manner. A core requirement is the solution of sparse systems of linear

equations; various HSL routines have been used and are compared.

Additionally, the requirement for fast small dense matrix solvers is

examined.

The sphere is an important setting for applied mathematics, yet the underlying approximation theory and numerical analysis needed for serious applications (such as, for example, global weather models) is much less developed than, for example, for the cube.

\\

\\

This lecture will apply recent developments in approximation theory on the sphere to two different problems in scientific computing.

\\

\\

First, in geodesy there is often the need to evaluate integrals using data selected from the vast amount collected by orbiting satellites. Sometimes the need is for quadrature rules that integrate exactly all spherical polynomials up to a specified degree $n$ (or equivalently, that integrate exactly all spherical harmonies $Y_{\ell ,k}(\theta ,\phi)$ with $\ell \le n).$ We shall demonstrate (using results of M. Reimer, I. Sloan and R. Womersley in collaboration with

W. Freeden) that excellent quadrature rules of this kind can be obtained from recent results on polynomial interpolation on the sphere, if the interpolation points (and thus the quadrature points) are chosen to be points of a so-called extremal fundamental system.

\\

\\

The second application is to the scattering of sound by smooth three-dimensional objects, and to the inverse problem of finding the shape of a scattering object by observing the pattern of the scattered sound waves. For these problems a methodology has been developed, in joint work with I.G. Graham, M. Ganesh and R. Womersley, by applying recent results on constructive polynomial approximation on the sphere. (The scattering object is treated as a deformed sphere.)

Continuing advances in computing technology provide the power not only to solve

increasingly large and complex process modeling and optimization problems, but also

to address issues concerning the reliability with which such problems can be solved.

For example, in solving process optimization problems, a persistent issue

concerning reliability is whether or not a global, as opposed to local,

optimum has been achieved. In modeling problems, especially with the

use of complex nonlinear models, the issue of whether a solution is unique

is of concern, and if no solution is found numerically, of whether there

actually exists a solution to the posed problem. This presentation

focuses on an approach, based on interval mathematics,

that is capable of dealing with these issues, and which

can provide mathematical and computational guarantees of reliability.

That is, the technique is guaranteed to find all solutions to nonlinear

equation solving problems and to find the global optimum in nonlinear

optimization problems. The methodology is demonstrated using several

examples, drawn primarily from the modeling of phase behavior, the

estimation of parameters in models, and the modeling, using lattice

density-functional theory, of phase transitions in nanoporous materials.