Comlab

From the point of view of a numerical analyst, I will describe some algorithms for:

• clustering data points based on pairwise similarity,
• reordering a sparse matrix to reduce envelope, two-sum or bandwidth,
• reordering nodes in a range-dependent random graph to reflect the range-dependency,

and point out some connections between seemingly disparate solution techniques. These datamining problems arise across a range of disciplines. I will mention a particularly new and important application from bioinformatics concerning the analysis of gene or protein interaction data.

Immersed interface methods have been developed for a variety of

differential equations on domains containing interfaces or irregular

boundaries. The goal is to use a uniform Cartesian grid (or other fixed

grid on simple domain) and to allow other boundaries or interfaces to

cut through this grid. Special finite difference formulas are developed

at grid points near an interface that incorporate the appropriate jump

conditions across the interface so that uniform second-order accuracy

(or higher) can be obtained. For fluid flow problems with an immersed

deformable elastic membrane, the jump conditions result from a balance

between the singular force imposed by the membrane, inertial forces if

the membrane has mass, and the jump in pressure across the membrane.

A second-order accurate method of this type for Stokes flow was developed

with Zhilin Li and more recently extended to the full incompressible

Navier-Stokes equations in work with Long Lee.

Feedback design for a second order control system leads to an

eigenstructure assignment problem for a quadratic matrix polynomial. It is

desirable that the feedback controller not only assigns specified

eigenvalues to the second order closed loop system, but also that the

system is robust, or insensitive to perturbations. We derive here new

sensitivity measures, or condition numbers, for the eigenvalues of the

quadratic matrix polynomial and define a measure of robustness of the

corresponding system. We then show that the robustness of the quadratic

inverse eigenvalue problem can be achieved by solving a generalized linear

eigenvalue assignment problem subject to structured perturbations.

Numerically reliable methods for solving the structured generalized linear

problem are developed that take advantage of the special properties of the

system in order to minimize the computational work required.

Electricity markets facilitate pricing and delivery of wholesale power.

Generators submit bids to an Independent System Operator (ISO) to indicate

how much power they can produce depending on price. The ISO takes these bids

with demand forecasts and minimizes the total cost of power production

subject to feasibility of distribution in the electrical network.

\\

\\

Each generator can optimise its bid using a bilevel program or

mathematical program with equilibrium (or complementarity) constraints, by

taking the ISOs problem, which contains all generators bid information, at

the lower level. This leads immediately to a game between generators, where

a Nash equilibrium - at which each generator's bid maximises its profit

provided that none of the other generators changes its bid - is sought.

\\

\\

In particular, we examine the idealised model of Berry et al (Utility

Policy 8, 1999), which gives a bilevel game that can be modelled as an

"equilibrium problem with complementarity constraints" or EPCC.

Unfortunately, like bilevel games, EPCCs on networks may not have Nash

equilibria in the (common) case when one or more of links of the network is

saturated (at maximum capacity). Nevertheless we explore some theory and

algorithms for this problem, and discuss the economic implications of

numerical examples where equilibria are found for small electricity

networks.

The seminar will address issues related to numerical simulation

of non-linear behaviour of solid materials to impact loading.

The kinematic and constitutive aspects of the transition from

continuum to discontinuum will be presented as utilised

within an explicit finite element development framework.

Material softening, mesh sensitivity and regularisation of

solutions will be discussed.

To numerical analysts and other applied mathematicians Jacobians and Hessians

are matrices, i.e. rectangular arrays of numbers or algebraic expressions.

Possibly taking account of their sparsity such arrays are frequently passed

into library routines for performing various computational tasks.

\\

\\

A central goal of an activity called automatic differentiation has been the

accumulation of all nonzero entries from elementary partial derivatives

according to some variant of the chainrule. The elementary partials arise

in the user-supplied procedure for evaluating the underlying vector- or

scalar-valued function at a given argument.

\\

\\

We observe here that in this process a certain kind of structure that we

call "Jacobian scarcity" might be lost. This loss will make the subsequent

calculation of Jacobian vector-products unnecessarily expensive.

Instead we advocate the representation of the Jacobian as a linear computational

graph of minimal complexity. Many theoretical and practical questions remain unresolved.

Condition numbers are a central concept in numerical analysis.

They provide a natural parameter for studying the behavior of

algorithms, as well as sensitivity and geometric properties of a problem.

The condition number of a problem instance is usually a measure

of the distance to the set of ill-posed instances. For instance, the

classical Eckart and Young identity characterizes the condition

number of a square matrix as the reciprocal of its relative distance

to singularity.

\\

\\

We present concepts of conditioning for optimization problems and

for more general variational problems. We show that the Eckart and

Young identity has natural extension to much wider contexts. We also

discuss conditioning under the presence of block-structure, such as

that determined by a sparsity pattern. The latter has interesting

connections with the mu-number in robust control and with the sign-real

spectral radius.

The simulation of sedimentary basins aims at reconstructing its historical

evolution in order to provide quantitative predictions about phenomena

leading to hydrocarbon accumulations. The kernel of this simulation is the

numerical solution of a complex system of time dependent, three

dimensional partial differential equations of mixed parabolic-hyperbolic

type in highly heterogeneous media. A discretisation and linearisation of

this system leads to large ill-conditioned non-symmetric linear systems

with three unknowns per mesh element.

\\

\\

In the seminar I will look at different preconditioning approaches for

these systems and at their parallelisation. The most effective

preconditioner which we developed so far consists in three stages: (i) a

local decoupling of the equations which (in addition) aims at

concentrating the elliptic part of the system in the "pressure block'';

(ii) an efficient preconditioning of the pressure block using AMG; (iii)

the "recoupling'' of the equations. Numerical results on real case

studies, exhibit (i) a significant reduction of sequential CPU times, up

to a factor 5 with respect to the current ILU(0) preconditioner, (ii)

robustness with respect to physical and numerical parameters, and (iii) a

speedup of up to 4 on 8 processors.

Rapidly oscillating problems, whether differential equations or

integrals, ubiquitous in applications, are allegedly difficult to

compute. In this talk we will endeavour to persuade the audience that

this is false: high oscillation, properly understood, is good for

computation! We describe methods for differential equations, based on

Magnus and Neumann expansions of modified systems, whose efficacy

improves in the presence of high oscillation. Likewise, we analyse

generalised Filon quadrature methods, showing that also their error

sharply decreases as the oscillation becomes more rapid.

In many applications of interest, such as the conformational

dynamics of molecules, large deterministic systems can exhibit

stochastic behaviour in a relative small number of coarse-grained

variables. This kind of dimension reduction, from a large deterministic

system to a smaller stochastic one, can be very useful in understanding

the problem. Whilst the subject of statistical mechanics provides

a wealth of explicit examples where stochastic models for coarse

variables can be found analytically, it is frequently the case

that applications of interest are not amenable to analytic

dimension reduction. It is hence of interest to pursue computational

algorithms for such dimension reduction. This talk will be devoted

to describing recent work on parameter estimation aimed at

problems arising in this context.

\\

\\

Joint work with Raz Kupferman (Jerusalem) and Petter Wiberg (Warwick)