Past

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.

This talk reviews some recent joint work with Michael Eiermann and Olaf

Schneider which introduced a framework for analyzing some popular

techniques for accelerating restarted Krylov subspace methods for

solving linear systems of equations. Such techniques attempt to compensate

for the loss of information due to restarting methods like GMRES, the

memory demands of which are usually too high for it to be applied to

large problems in unmodified form. We summarize the basic strategies which

have been proposed and present both theoretical and numerical comparisons.

Support Vector Machines are a new and very promising approach to

machine learning. They can be applied to a wide range of tasks such as

classification, regression, novelty detection, density estimation,

etc. The approach is motivated by statistical learning theory and the

algorithms have performed well in practice on important applications

such as handwritten character recognition (where they currently give

state-of-the-art performance), bioinformatics and machine vision. The

learning task typically involves optimisation theory (linear, quadratic

and general nonlinear programming, depending on the algorithm used).

In fact, the approach has stimulated new questions in optimisation

theory, principally concerned with the issue of how to handle problems

with a large numbers of variables. In the first part of the talk I will

overview this subject, in the second part I will describe some of the

speaker's contributions to this subject (principally, novelty

detection, query learning and new algorithms) and in the third part I

will outline future directions and new questions stimulated by this

research.

The flow in a cylinder with a rotating endwall has continued to

attract much attention since Vogel (1968) first observed the vortex

breakdown of the central core vortex that forms. Recent experiments

have observed a multiplicity of unsteady states that coexist over a

range of the governing parameters. In spite of numerous numerical and

experimental studies, there continues to be considerable controversy

with fundamental aspects of this flow, particularly with regards to

symmetry breaking. Also, it is not well understood where these

oscillatory states originate from, how they are interrelated, nor how

they are related to the steady, axisymmetric basic state.

\\

\\

In the aspect ratio (height/radius) range 1.6 2.8. An efficient and

accurate numerical scheme is presented for the three-dimensional

Navier-Stokes equations in primitive variables in a cylinder. Using

these code, primary and secondary bifurcations breaking the SO(2)

symmetry are analyzed.

\\

\\

We have located a double Hopf bifurcation, where an axisymmetric limit

cycle and a rotating wave bifurcate simultaneously. This codimension-2

bifurcation is very rich, allowing for several different scenarios. By

a comprehensive two-parameter exploration about this point we have

identified precisely to which scenario this case corresponds. The mode

interaction generates an unstable two-torus modulate rotating wave

solution and gives a wedge-shaped region in parameter space where the

two periodic solutions are both stable.

\\

\\

For aspect ratios around three, experimental observations suggest that

the first mode of instability is a precession of the central vortex

core, whereas recent linear stability analysis suggest a Hopf

bifurcation to a rotating wave at lower rotation rates. This apparent

discrepancy is resolved with the aid of the 3D Navier-Stokes

solver. The primary bifurcation to an m=4 traveling wave, detected by

the linear stability analysis, is located away from the axis, and a

secondary bifurcation to a modulated rotating wave with dominant modes

m=1 and 4, is seen mainly on the axis as a precessing vortex breakdown

bubble. Experiments and the linear stability analysis detected

different aspects of the same flow, that take place in different

spatial locations.

We combine linear algebra techniques with finite element techniques to obtain a reliable stopping criterion for the Conjugate Gradient algorithm. The finite element method approximates the weak form of an elliptic partial differential equation defined within a Hilbert space by a linear system of equations A x = b, where A is a real N by N symmetric and positive definite matrix. The conjugate gradient method is a very effective iterative algorithm for solving such systems. Nevertheless, our experiments provide very good evidence that the usual stopping criterion based on the Euclidean norm of the residual b - Ax can be totally unsatisfactory and frequently misleading. Owing to the close relationship between the conjugate gradient behaviour and the variational properties of finite element methods, we shall first summarize the principal properties of the latter. Then, we will use the recent results of [1,2,3,4]. In particular, using the conjugate gradient, we will compute the information which is necessary to evaluate the energy norm of the difference between the solution of the continuous problem, and the approximate solution obtained when we stop the iterations by our criterion.

Finally, we will present the numerical experiments we performed on a selected ill-conditioned problem.

References

• [1] M. Arioli, E. Noulard, and A. Russo, Vector Stopping Criteria for Iterative Methods: Applications to PDE's, IAN Tech. Rep. N.967, 1995.
• [2] G.H. Golub and G. Meurant, Matrices, moments and quadrature II; how to compute the norm of the error in iterative methods, BIT., 37 (1997), pp.687-705.
• [3] G.H. Golub and Z. Strakos, Estimates in quadratic formulas, Numerical Algorithms, 8, (1994), pp.~241--268.
• [4] G. Meurant, The computation of bounds for the norm of the error in the conjugate gradient algorithm, Numerical Algorithms, 16, (1997), pp.~77--87.

Boundary layers are often studied with no pressure gradient

or with an imposed pressure gradient. Either of these assumptions

can lead to difficulty in obtaining solutions. A major advance in fluid

dynamics last century (1969) was the development of a triple deck

formulation for boundary layers where the pressure is not

specified but emerges through an interaction between

boundary layer and the inviscid outer flow. This has given rise to

new computational problems and computations have in turn

fed ideas back into theoretical developments. In this survey talk

based on my new book, I will look at three problems:

flow past a plate, flow separation and flow in channels

and discuss the interaction between theory and computation

in advancing boundary layer theory.

This talk reviews some theoretical and practical aspects

of incompressible flow modelling using finite element approximations

of the (Navier-) Stokes equations.

The infamous Q1-P0 velocity/pressure mixed finite element approximation

method is discussed. Two practical ramifications of the inherent

instability are focused on, namely; the convergence of the approximation

with and without regularisation, and the behaviour of fast iterative

solvers (of multigrid type) applied to the pressure Poisson system

that arises when solving time-dependent Navier-Stokes equations

using classical projection methods.

\\

\\

This is joint work with David Griffiths from the University of Dundee.

We examine the convergence characteristics of iterative methods based

on a new preconditioning operator for solving the linear systems

arising from discretization and linearization of the Navier-Stokes

equations. With a combination of analytic and empirical results, we

study the effects of fundamental parameters on convergence. We

demonstrate that the preconditioned problem has an eigenvalue

distribution consisting of a tightly clustered set together with a

small number of outliers. The structure of these distributions is

independent of the discretization mesh size, but the cardinality of

the set of outliers increases slowly as the viscosity becomes smaller.

These characteristics are directly correlated with the convergence

properties of iterative solvers.

The talk will focus on solution methods for augmented linear systems of

the form

\\

\\

$[ A B ][x] = [b] [ B' 0 ][y] [0]$.

\\

\\

Augmented linear systems of this type arise in several areas of

numerical applied mathematics including mixed finite element / finite

difference discretisations of flow equations (Darcy flow and Stokes

flow), electrical network simulation and optimisation. The general

properties of such systems are that they are large, sparse and

symmetric, and efficient solution techniques should make use of the

block structure inherent in the system as well as of these properties.

\\

\\

Iterative linear solution methods will be described that

attempt to take advantage of the structure of the system, and

observations on augmented

systems, in particular the distribution of their eigenvalues, will be

presented which lead to further iterative methods and also to

preconditioners for existing solution methods. For the particular case

of Darcy flow, comments on properties of domain decomposition methods

of additive Schwarz type and similarities to incomplete factorisation

preconditioners will be made.

I am going to show that all self-scaled barriers for the

cone of symmetric positive semidefinite matrices are of the form

$X\mapsto -c_1\ln\det X +c_0$ for some constants $c_1$ > $0,c_0 \in$ \RN.

Equivalently one could state say that all such functions may be

obtained via a homothetic transformation of the universal barrier

functional for this cone. The result shows that there is a certain

degree of redundancy in the axiomatic theory of self-scaled barriers,

and hence that certain aspects of this theory can be simplified. All

relevant concepts will be defined. In particular I am going to give

a short introduction to the notion of self-concordance and the

intuitive ideas that motivate its definition.

Interval arithmetic is a way to produce guaranteed enclosures of the

results of numerical calculations. Suppose $f(x)$ is a real

expression in real variables $x= (x_1, \ldots, x_n)$, built up from

the 4 basic arithmetic operations and other 'standard functions'. Let

$X_1, \ldots, X_n$ be (compact) real intervals. The process of {\em

interval evaluation} of $f(X_1, ..., X_n)$ replaces each real

operation by the corresponding interval operation wherever it occurs

in $f$, e.g. $A \times B$ is the smallest interval containing $\{a

\times b \mid a \in A, b \in B\}$.

\\

\\

As is well known, it yields a guaranteed enclosure for the true range

$\{f(x_1, \ldots, x_n) \mid x_1 \in X_1, \ldots, x_n \in X_n\}$,

provided no exceptions such as division by (an interval containing)

zero occur during evaluation.

\\

\\

Interval arithmetic takes set inputs and produces set outputs. Noting

this, we show there is a consistent way to extend arithmetic to $R^* =

R \cup \{-\infty, +\infty\}$, such that interval evaluation continues

to give enclosures, and there are {\em no exceptions}. The basic

ideas are: the usual set-theory meaning of evaluating a relation at a

set; and taking topological closure of the graph of a function in a

suitable $(R^{*})^n$. It gives rigorous meaning to intuitively

sensible statements like $1/0 = \{-\infty, +\infty\}$, $0/0 = R^*$

(but $(x/x)_{|x=0} = 1$), $\sin(+\infty) = [-1,1]$, etc.

\\

\\

A practical consequence is that an exception-free floating-point

interval arithmetic system is possible. Such a system is implemented

at hardware level in the new Sun Fortran compiler, currently on

beta-release.

The Cholesky factorization approach to solving the symmetric definite generalized eigenvalue problem

$Ax = \lambda Bx$, where $A$ is symmetric and $B$ is symmetric positive definite, computes a Cholesky factorization $B = LL^T$ and solves the equivalent standard symmetric eigenvalue problem $C y = \l y$ where $C = L^{-1} A L^{-T}$. Provided that a stable eigensolver is used, standard error analysis says that the computed eigenvalues are exact for $A+\dA$ and $B+\dB$ with $\max( \normt{\dA}/\normt{A}, \normt{\dB}/\normt{B} )$

bounded by a multiple of $\kappa_2(B)u$, where $u$ is the unit roundoff. We take the Jacobi method as the eigensolver and explain how backward error bounds potentially much smaller than $\kappa_2(B)u$ can be obtained.

To show the practical utility of our bounds we describe a vibration problem from structural engineering in which $B$ is ill conditioned yet the error bounds are small. We show how, in cases of instability, iterative refinement based on Newton's method can be used to produce eigenpairs with small backward errors.

Our analysis and experiments also give insight into the popular Cholesky--QR method used in LAPACK, in which the QR method is used as the eigensolver.

In contrast to spectra, pseudospectra of large Toeplitz matrices

behave as nicely as one could ever expect. We demonstrate some

basic phenomena of the asymptotic distribution of the spectra

and pseudospectra of Toeplitz matrices and show how by employing

a few simple $C^*$-algebra arguments one can prove rigorous

convergence results for the pseudospectra. The talk is a survey

of the development since a 1992 paper by Reichel and Trefethen

and is not addressed to specialists, but rather to a general

mathematically interested audience.

This lecture is about the extension of our CAD-Free platform to the simulation and sensitivity analysis for design and control of multi-model configurations. We present the different ingredients of the platform using simple models for the physic of the problem. This presentation enables for an easy evaluation of coupling, design and control strategies.

\\

\\

Sensitivity analysis has been performed using automatic differentiation in direct or reverse modes and finite difference or complex variable methods. This former approach is interesting for cases where only one control parameter is involved as we can evaluate the state and sensitivity in real time with only one evaluation.

\\

\\

We show that for some class of applications, incomplete sensitivities can be evaluated dropping the state dependency which leads to a drastic cost reduction.

\\

\\

The design concerns the structural characteristic of the model and our control approach is based in perturbating the inflow incidence by a time dependent flap position described by a dynamic system encapsulating a gradient based minimization algorithm expressed as dynamic system. This approach enables for the definition of various control laws for different minimization algorithm.

\\

\\

We present dynamic minimization algorithms based on the coupling of several heavy ball method. This approach enables for global minimization at a cost proportional to the number of balls times the cost of one steepest descent minimization.

\\

\\

Two and three dimension flutter problem simulation and control are presented and the sensitivity of the different parameters in the model is discussed.

This seminar has been cancelled.

There is a growing interest in the study of periodic phenomena in

large-scale nonlinear dynamical systems. Often the high-dimensional

system has only low-dimensional dynamics, e.g., many reaction-diffusion

systems or Navier-Stokes flows at low Reynolds number. We present an

approach that exploits this property in order to compute branches of

periodic solutions of the large system of ordinary differential

equations (ODEs) obtained after a space discretisation of the PDE. We

call our approach the Newton-Picard method. Our method is based on the

recursive projection method (RPM) of Shroff and Keller but extends this

method in many different ways. Our technique tries to combine the

performance of straightforward time integration with the advantages of

solving a nonlinear boundary value problem using Newton's method and a

direct solver. Time integration works well for very stable limit

cycles. Solving a boundary value problem is expensive, but works also

for unstable limit cycles.

\\

\\

We will present some background material on RPM. Next we will explain

the basic features of the Newton-Picard method for single shooting. The

linearised system is solved by a combination of direct and iterative

techniques. First, we isolate the low-dimensional subspace of unstable

and weakly stable modes (using orthogonal subspace iteration) and

project the linearised system on this subspace and on its

(high-dimensional) orthogonal complement. In the high-dimensional

subspace we use iterative techniques such as Picard iteration or GMRES.

In the low-dimensional (but "hard") subspace, direct methods such as

Gaussian elimination or a least-squares are used. While computing the

projectors, we also obtain good estimates for the dominant,

stability-determining Floquet multipliers. We will present a framework

that allows us to monitor and steer the convergence behaviour of the

method.

\\

\\

RPM and the Newton-Picard technique have been developed for PDEs that

reduce to large systems of ODEs after space discretisation. In fact,

both methods can be applied to any large system of ODEs. We will

indicate how these methods can be applied to the discretisation of the

Navier-Stokes equations for incompressible flow (which reduce to an

index-2 system of differential-algebraic equations after space

discretisation when written in terms of velocity and pressure.)

\\

\\

The Newton-Picard method has already been extended to the computation

of bifurcation points on paths of periodic solutions and to multiple

shooting. Extension to certain collocation and finite difference

techniques is also possible.

Krylov subspace methods offer good possibilities for the solution of

large sparse linear systems of equations.For general systems, some of

the popular methods often show an irregular type of convergence

behavior and one may wonder whether that could be improved or not. Many

suggestions have been made for improvement and the question arises

whether these corrections are cosmetic or not. There is also the

question whether the irregularity shows inherent numerical instability.

In such cases one should take extra care in the application of

smoothing techniques. We will discuss strategies that work well and

strategies that might have been expected to work well.

This work forms part of a larger research project to develop efficient

low-dissipative high-order numerical techniques for high-speed

turbulent flow simulation, including shock wave interactions with

turbulence. The requirements on a numerical method are stringent.For

the turbulence the method must be capable of resolving accurately a

wide range of length scales, whilst for shock waves the method must be

stable and not generate excessive local oscillations. Conventional

methods are either too dissipative, or incapable of shock capturing.

Higher-order ENO, WENO or hybrid schemes are too expensive for

practical computations. Previous work of Yee, Sandham & Djomehri

(1999) developed high-order shock-capturing schemes which minimize the

use of numerical dissipation away from shock

waves. The objective of the present study is to further minimize the

use of numerical dissipation for shock-free compressible turbulence

simulations.

In 1983 Pantoja described a stagewise construction of the exact Newton

direction for a discrete time optimal control problem. His algorithm

requires the solution of linear equations with coefficients given by

recurrences involving second derivatives, for which accurate values are

therefore required.

\\

\\

Automatic differentiation is a set of techniques for obtaining derivatives

of functions which are calculated by a program, including loops and

subroutine calls, by transforming the text of the program.

\\

\\

In this talk we show how automatic differentiation can be used to

evaluate exactly the quantities required by Pantoja's algorithm,

thus avoiding the labour of forming and differentiating adjoint

equations by hand.

\\

\\

The cost of calculating the newton direction amounts to the cost of

solving one set of linear equations, of the order of the number of

control variables, for each time step. The working storage cost can be made

smaller than that required to hold the solution.