The Australian National University
Mathematical Sciences Institute (MSI)
Events - Special Year 2005
document location: http://wwwmaths.anu.edu.au/events/sy2005/program2.html

2005 Special Year on Computational Mathematics

Wednesday Thursday Friday
Hutchinson
Womersley
Osborne
Hudson
Ryan
Golub
Smola
Griewank
Bock
Bakin
Saunders
Björck
Turlach
Moré
Smyth
Haviv
Filar
Brent
Burrage
Howlett
McLean
Prvan
Gordon
Stokes
Roberts
Manton
Sloan
Kuo
Butcher
Hegland
Sidje
Podhaisky

Newsam
Nuyens


Wednesday, September 21


 8:30 – 8:50  Registration


 8:50 – 9:00  Opening remarks


 9:00 – 9:30  Mike Hutchinson
Australian National University
Locally adaptive gridding of elevation data
 9:30 – 10:00  Malcolm Hudson
Macquarie University
Block Fisher scoring optimization of penalized likelihoods in emission tomography
 10:00 – 10:40  Alex Smola
NICTA
Nonparametric tests for distributions
 10:40 – 11:00  Morning Tea
 11:00 – 11:40  Sergey Bakin
Suncorp
Regression tree ensemble models
 11:40 – 12:20  Berwin Turlach
University of Western Australia
On homotopy algorithms in statistics
 12:20 – 13:30  Lunch break
 13:30 – 14:00  Moshe Haviv
The Hebrew University of Jerusalem
On singularly perturbed Markov chains
 14:00 – 14:30  Kevin Burrage
University of Queensland
Multiscale modelling of cellular processes
 14:30 – 15:00 
Tania Prvan
Macquarie University
Thin film models in a stochastic setting
 15:00 – 15:30 
Steve Roberts
Australian National University
An application of sparse grids to the estimation of probability density functions
 15:30 – 16:00  Afternoon Tea
 16:00 – 16:30 
Frances Kuo
University of New South Wales
Importance sampling plays a crucial role for quasi-Monte Carlo methods
 16:30 – 17:00 
Roger Sidje
University of Queensland
A QR-based tridiagonalization algorithm for nonsymmetric matrices
 17:00 – 17:30 
Garry Newsam
Defence, Science & Technology Org. (DSTO)
Domain decomposition relations as a source of fast algorithms for evaluating integral operators

Thursday, September 22



 9:00 – 9:30  Rob Womersley
University of New South Wales
Distributing points on manifolds
 9:30 – 10:00  David Ryan
University of Auckland
The use of severely limited subsequence in set partitioning optimisation
 10:00 – 10:40  Andreas Griewank
Humboldt Universität Berlin (member of MATHEON)
On the computational complexity of Jacobians and Hessians
 10:40 – 11:00  Morning Tea
 11:00 – 11:40 
Jorge Moré
Argonne National Laboratory
Evaluating optimization
 11:40 – 12:20 
Michael Saunders
Stanford University
Least-squares methods in optimization
 12:20 – 13:30  Lunch break
 13:30 – 14:00 
Jerzy Filar
University of South Australia
Doubly stochastic matrices, the Hamiltonian cycle problem & a control theoretic perspective of static optimization problems.
 14:00 – 14:30 
Phil Howlett
University of South Australia
Calculation of optimal driving strategies for freight trains
 14:30 – 14:50 
Dan Gordon
Constraint Technologies International
Airline optimisation problems – an industry perspective
 14:50 – 15:20 
Jonathan Manton
Australian National University
On the theory of optimisation on manifolds
 15:20 – 15:50 Afternoon Tea
 15:50 – 16:20 
John Butcher
University of Auckland
Towards efficient general linear methods
 16:20 – 16:40 
Helmut Podhaisky
Martin Luther Universität
On variable order for general linear methods in Nordsieck form
 16:40 – 17:10 
Dirk Nuyens
Katholieke Universiteit Leuven, Belgium
On the contruction of rank-1 lattice rules

Following thursday's program, there will be a symposium dinner at University House at 18:30 for 19:00.


Friday, September 23


                                                                                                     

 9:00 – 9:30  Mike Osborne
Australian National University
Multiple shooting revisited
 9:30 – 10:00  Gene Golub
Stanford University
Matrices, moments and quadrature with emphasis on least squares problems
 10:00 – 10:40  Georg Bock
University of Heidelberg
From parameter estimation to the design of optimal experiments in differential algebraic equations
 10:40 – 11:00  Morning Tea
 11:00 – 11:40 
Åke Björck
Linköping University
Recent developments in least squares computations — Bidiagonal decomposition and least squares
 11:40 – 12:20 
Gordon Smyth
Walter & Eliza Hall Institute of Medical Research
Random effects, heteroscedastic regression, algorithms and microarray data analysis
 12:20 – 13:30  Lunch break
 13:30 – 14:00 
Richard Brent
Australian National University
Challenges in solving large sparse linear systems over finite fields
 14:00 – 14:30 
William McLean
University of New South Wales
The conditioning of boundary element equations
 14:30 – 14:50 
Nick Stokes
CSIRO
Material flows — matching micro-scale computations with macroscopic effects
 15:00 – 15:30 Afternoon Tea
 15:30 – 16:00 
Ian Sloan
University of New South Wales
The curse of dimensionality
 16:00 – 16:30 
Markus Hegland
Australian National University
Resolution enhancement of spectra using differentiation



Abstracts





Mike Hutchinson, Australian National University
Locally adaptive gridding of elevation data

Finite difference discretisation, based on a regular two-dimensional grid, has provided the means for efficiently calculating close approximations to bivariate thin plate smoothing splines. When a simple nested grid SOR iteration is used, the number of operations can be made optimal, in the sense that they are proportional to the number of grid points, regardless of grid spacing (Hutchinson 1989, 2000, 2003). Problems involving millions of data points and grid points can be routinely solved on modest workstations. A particular additional advantage of the finite difference approach is that a variety of locally adaptive constraints can be straightforwardly applied to satisfy known process based conditions. This is particularly appropriate for generating two dimensional grids of elevation based on a variety of topographic data, including point heights, elevation contours, streamlines, lakes and cliff lines. The latter four of these data types have implicit and explicit information that can be used to help define the shape of the interpolated elevation grid. Moreover a locally adaptive drainage enforcement algorithm can be added to produce digital elevation models with connected drainage structure. Theoretical and practical issues in constructing this method are described, including the vectorisation of most of the calculations. The method also simultaneously minimises a weighted residual sum of squares that respects the natural discretisation error that arises with a finite difference discretisation. This has led to a method for optimising the spatial resolution of digital elevation models to the information content of source topographic data. A locally adaptive modification to the usual minimum curvature roughness penalty has also been suggested to better respect drainage structure (Hutchinson 1996, 2000).

References:
Hutchinson, M. F. (1989) A new method for gridding elevation and streamline data with automatic removal of spurious pits. Journal of Hydrology 106: 211-232.
Hutchinson, M. F. (1996) A locally adaptive approach to the interpolation of digital elevation models.
Hutchinson, M. F. (2000) Optimising the degree of data smoothing for locally adaptive finite element bivariate smoothing splines. ANZIAM Journal 42(E), C774-C796.
Hutchinson, M. F. (2003) ANUDEM Version 5.1 Centre for Resource and Environmental Studies, Australian National University.



Malcolm Hudson, Macquarie University
Block Fisher scoring optimization of penalized likelihoods in emission tomography

Following the contributions of Osborne (ISI Review, 1992) in generalized linear models, we study properties of Fisher scoring and its block iterative implementation for density estimation utilizing Poisson models for indirect projection counts in emission tomography. Adopting a quadratic penalty, we provide a duality based decomposition of the penalized likelihood which provides a block-iterative solution. This is joint work with Jun Ma.

Alex Smola, NICTA
Nonparametric tests for distributions
Presentation (pdf)

Linear criteria give easy to test necessary conditions for testing independence and identity between two distributions. More to the point, it is a necessary condition for independence that the covariance between random variables vanishes. Equally, it is a necessary condition for identity that the means of random variables agree between distributions. These conditions can be turned into sufficient criteria by extending this requirement to hold for a set of functions of the random variables, provided that this class is sufficiently rich. By using Reproducing Kernel Hilbert Spaces we obtain simple tests which have nonparametric power yet are very easy to compute, essentially by taking sums over entries of matrices without requiring a separate optimization stage.

Sergey Bakin, Suncorp
Regression tree ensemble models
Presentation (ppt)

Tree-based models have gained a good reputation in macnine learning and applied statistics due to their robustness, interpretability and high scalability of the algorithms used to generate such models. A tree ensemble model is a collection of trees. A prediction from such a model is obtained by combining predictions from individual trees. It has been demonstrated that the prediction accuracy of an ensemble of trees is often much better than accuracy of any individual tree.

In this presentation I will concentrate on regression tree ensemble methods for modelling data with non-normal distribution of the target (response) variable.

Berwin Turlach, University of Western Australia
On homotopy algorithms in statistics
Presentation (pdf)

The LASSO, proposed by Tibshirani (1996, JRSS B), is a colourful name for, essentially, minimising the residual sum of squares of a least squares problem while imposing a constraint (or penalty) on the L1 norm of the parameter vector. However, this methodology has some interesting shrinkage and variables selection features. Various authors have since studied (and extended) this methodology; among them Osborne et al. (2000, IMA J. Num. Anal.) who showed that the solution of the LASSO problem, as a function of the constraint, is piecewise linear and continuous. Osborne et al. also describe a homotopy algorithm for calculating the complete solution path of the LASSO. This algorithm was later popularised in the statistical literature by Efron et al. (2004, Ann. Stat.) who develop the least angle regression (LARS) framework which, amongst other results, provides further insight into how the LASSO works and how it is related to other variable selection techniques. In the wake of Efron et al. some research effort in the statistics community was directed towards developing homotopy algorithms for various statistical problems. Time permitting, we will review some of these algorithm before we discuss the latest addition, a homotopy algorithm that calculates the complete solution path for a particular extension of the LASSO to the multivariate regression case proposed by Turlach et al. (2005, Technometrics).

Moshe Haviv, The Hebrew University of Jerusalem
On singularly perturbed Markov chains

Nothing structurally happens in case that transition probabilities in an ergodic Markov chain are slightly perturbed as all measures are continuous with the perturbation. This is not the case when the state space can be decomposed into a number of ergodic classes plus a number of transient states, and when the perturbation coupled them all. For example, the stationary distribution is now uniquely defined and mean passage times between states belonging to different classes are now well defined. The talk will survey some results which are known and some which are under construction. Special emphasis will be given to the issue of time scales.

Kevin Burrage, University of Queensland
Multiscale modelling of cellular processes
Presentation (pdf)

The Stochastic Simulation Algorithm (SSA) of Dan Gillespie is a crucial technique for simulating the interactions of small numbers of molecules in cellular environments and in particular genetic regulation. However, simulation codes based on this approach can be very slow because of the possibly very small time steps involved. In this talk we will discuss some multi-scaled approaches suited to the modelling of the kinetics of cellular systems when there is a mix of small numbers and large numbers of molecules.

Tania Prvan, Macquarie University
Thin film models in a stochastic setting
Presentation (pdf)

Dunn & Tichenor (1988) proposed a class of differential equation models to describe the phenomenon of transient sink behaviour for organic emissions exhibited by interior surface films in state-of-the-art emission test chambers. Data from a particular application is used to exemplify the use of a model selection scheme which embeds the derived models within a class of stochastic differential equations. These embeddings have the property that the quality of model fit varies inversely with the strength of the stochastic forcing term. Results of this modelling application are discussed.

Steve Roberts, Australian National University
An application of sparse grids to the estimation of probability density functions
Presentation (pdf)

In this talk we will provide an overview of the approximation properties of sparse grids and then show how sparse grids can be used to provide an efficient method for the estimation of high dimensional probability density functions (pdfs). We take a simple approach and obtain an estimate of a pdf using an L2+Mixed Sobolov semi-norm projection method. This method in combination with the work of Markus Hegland on adaptive sparse grids has the potential to provide a effcient method for estimating pdfs in very high dimensions.

Fraces Kuo, University of New South Wales
Importance sampling plays a crucial role for quasi-Monte Carlo methods
Presentation (pdf)

Many practical problems in statistics and mathematical finance involve integrals with hundreds or even thousands of dimensions. These integrals often arise from multivariate expected values. In some cases the actual integrals are the goal (e.g. option pricing), while in others only the best parameters for the model were required (e.g. maximum likelihood problems). The common strategy is to use Monte Carlo methods with importance sampling, i.e., an integral is approximated by an average of function values at a number of points sampled randomly from a certain probability density.

In the past decade or so, quasi-Monte Carlo methods have emerged as an efficient alternative to Monte Carlo methods. Instead of generating the sample points randomly, these points are constructed deterministically over the unit cube [0,1]d to be as uniform as possible. However, since all the underlying theory behind quasi-Monte Carlo methods reside on the unit cube while most integrals in practice are defined over the unbounded domain Rd an appropriate transformation is required to bring the problem into the unit cube. This transformation plays a crucial role in the effectiveness of quasi-Monte Carlo methods since it determines the features of the transformed integrand. As it turns out, this is equivalent to importance sampling for Monte Carlo methods.

This talk is devoted to a discussion of various techniques for choosing a good transformation, with a focus on some maximum likelihood problems in statistics. It is based on a joint project with William Dunsmuir, Matt Wand, Rob Womersley, and Ian Sloan.

Roger Sidje, University of Queensland
A QR-Based tridiagonalization algorithm for nonsymmetric matrices
Presentation (pdf)

The tridiagonalization of a general square matrix represents the most compact similarity reduction that can be computed directly. Attempting to reduce further (diagonalization or bidiagonalization) implies the retrieval of eigenvalues, which can only be done iteratively in general, unless the order of the matrix is under five. Consequently, finite tridiagonalization algorithms tend to be prone to numerical instabilities. The tridiagonal representation is too compact for its own good.

The talk will first outline a general framework for tridiagonalization algorithms. We subsequently present a new elimination technique combining orthogonal similarity transformations that are stable. We then discuss recovery techniques when serious breakdowns are encountered. It is established that the new algorithm is nearly optimal in minimizing roundoff errors. These are recent developments that have just appeared in the SIAM Journal on Matrix Analysis and Applications.

Applications of this study include eigenvalue calculation and the approximation of matrix functions.

Garry Newsam, Defence, Science & Technology Org. (DSTO)
Domain decomposition relations as a source of fast algorithms for evaluating integral operators
Presentation (pdf)

For a number of important integral operators it is possible to derive simple equations that define the operator over a given domain in terms of the output of the operator over sub domains: such equations in turn provide a basis for constructing fast, divide-and-conquer algorithms for computing the transform. The FFT can be derived from one such equation satisfied by the Fourier transform, but the talk will instead illustrate this idea through analysis of a simple equation expressing the Radon transform over a square in terms of the transforms over the four component sub squares. Approximation of this equation, e.g. by collocation or Galerkin methods, across a range of scales implicitly discretises the operator into a product of sparse matrices. This makes it possible to compute an approximation to the Radon transform of an N×N image in O(N2)log N operations, where the constant depends on the order of the approximation and the accuracy desired. This is an improvement over most existing algorithms (e.g. the Radon transform routine in Matlab) as these require O(N3) operations. The technique has several bonuses: the Radon transforms over a hierarchy of sub domains are provided for free in the course of the calculation (a big advantage over competing, Fourier-based, fast algorithms for computing the Radon transform); while the matrix product representation of the discretisation also provides fast algorithms for computing adjoints (in this case the back-projection operator) and inverses. Computations illustrating the algorithm and achievable accuracies will be included.

Rob Womersley, University of New South Wales
Distributing points on manifolds
Presentation (pdf)

There are many criteria that can be used to distribute points on manifolds such as spheres and tori. This includes minimizing the Riesz s-energy, maximizing the sum of distances, minimizing the mesh ratio or properties of a basis matrix used for interpolation or to get cubature weights. The choice is naturally influenced by the aim, whether it be to discretize the manifold by a “uniformly distributed” set of points, to find good points for polynomial or radial basis interpolation, or to find points and weights for numerical integration over the manifold.

All of these can be posed as nonlinear optimization problems, which typically have many local optima close to a global optimum. This talk will look at optimizing several of these criteria on sphere and tori. In part the work is motivated by the remarkable result of Hardin and Saff that the minimal Riesz s-energy point sets are asymptotically uniformly distributed for s larger than the dimension of the manifold for a wide class of manifolds.

David Ryan, University of Auckland
The use of severely limited subsequence in set partitioning optimisation

Many real world scheduling applications can be modeled with the Set Partitioning optimisation model. Optimal integer solutions of these models exhibit the property of “unique subsequence”. Set Partitioning Problems possessing a “limited subsequence” property are known to be easier to solve but optimality may be compromised if the limited subsequences do not contain the optimal (unique) subsequences. Fortunately the limited subsequence property occurs naturally in many applications. In this talk we consider the use of limited subsequence formulations in which we severely restrict subsequences to improve the efficiency of the solution algorithm but then identify or generate additional subsequences which have the potential to improve solution quality.

Andreas Griewank, Humboldt-Universität Berlin (member of MATHEON)
On the computational complexity of Jacobians and Hessians
Presentation (pdf)

Most numerical algorithms for the solution of nonlinear problems use local approximations based on first and possibly second derivatives. In many calculations the effort for evaluating or approximating these derivative matrices dominates the overall computational cost.

For functions that are given by evaluation programs or otherwise specified as a composite of elemental functions, the chain rule based technique of algorithmic or automatic differentiation yields truncation error free derivative values at a priori bounded computational costs. However, as had been conjectured for at least a decade and just recently proven by Uwe Naumann, the task of minimizing the number of multiplications in computing Jacobians is NP hard.

Fortunately, the explicit accumulation of Jacobians and Hessians as real arrays can often be avoided, for example in Krylov solvers that require only matrix vector products. The attempt to provide these at minimal cost leads to another combinatorial optimization problem, whose complexity is not yet clear.

Jorge Moré, Argonne National Laboratory
Evaluating optimization
Presentation (pdf)

We introduce the concept of performance profiles to evaluate the behavior of optimization algorithms on model applications. We use case studies to illustrate the difficulties of developing and evaluating optimization software for the solution of large optimization problems. Our presentation centers on codes for solving constrained optimization problems, but our results have applicability to algorithms in most areas of scientific computing. The presentation centers on the results of a series of computational experiments designed to measure the performance of optimization solvers. Our aim is to answer a basic question: Are we making progress in optimization?

Michael Saunders, Stanford University
Least-squares methods in optimization
Presentation (pdf)

Least-squares methods are vital tools for data analysis. They also arise in algorithms for constrained optimization. We discuss the least-squares solvers LSQR, LSSOL, SQOPT, PDCO, and their application to various problems including LASSO, Fused LASSO, and Basis Pursuit De-Noising.

Jerzy Filar, University of South Australia
Doubly stochastic matrices, the Hamiltonian cycle problem & a control theoretic perspective of static optimization problems.
Presentation (pdf)

When seeking a numerical solution to a control problem, it is common to convert it to an approximate, static, optimization problem and then to solve the latter. The perspective promoted in this presentation is that there are situations where a reverse approach can prove fruitful. Namely, a difficult static optimization problem may, sometimes, be approximated by a dynamic (possibly stochastic) control problem and a solution sought with the help of the dynamic theories.

We illustrate the above perspective by considering the classical (NP-hard) problem of determining whether a given graph possesses a Hamiltonian cycle. Exploiting singularly perturbed controlled Markov chains and their fundamental matrices, we reduce the problem to that of minimising the variability of the “first return time” to the home node. The latter is a highly structured, non-convex, optimisation problem whose properties shed light on both the theoretical complexity of the problem and its algorithmic  tractability. A key recent innovation in this approach is the restriction to those controls that induce doubly stochastic Markov chains. This domain is convex, rich enough to include Hamiltonian cycles and structured enough to considerably simplify fundamental matrices of the relevant Markov chains.

Arguably, the problem of determining whether a given graph possesses a Hamiltonian cycle contains the essential difficulty of the famous “Travelling Salesman Problem”.  We claim that it is conceptually insightful, and algorithmically useful, to characterise this difficulty in terms of quantities such as the variability of returns (to the initial state) in a controlled stochastic process, or gaps between the minimal values of suitably constructed mathematical programming problems.

This is joint work with Vivek S Borkar, Tata Institute for Fundamental Research (Mumbai) and Vladimir Ejoy of CIAM, UniSA.

Phil Howlett, University of South Australia
Calculation of optimal driving strategies for freight trains
Presentation (pdf)

The problem of finding a driving strategy that minimises fuel consumption subject to completing the required journey within a given time is now well understood. A key paper by Asnis et al. which minimised mechanical energy subject to uniformly bounded acceleration was published in 1985 and a similar approach was used by Howlett in 1990. A more realistic model which sought to minimise fuel consumption and allowed only discrete throttle settings was formulated by Benjamin et al. at UniSA and subsequently solved in a sequence of papers by Howlett, Pudney and Cheng during the period 1991 – 1999. A corresponding formulation and solution with continuous throttle settings was given by Khmelnitsky in 1994 but was not published until 2000.

Suffice it to say that on level track the solution is a power, speedhold, coast, brake strategy in which the dominant speedhold mode is actually a singular control. This solution is preserved on non-steep track but on a track with steep inclines and declines where the desired holding speed cannot be maintained it is necessary to anticipate the steep grades by switching to power before the steep incline is reached and by switching to coast before the steep decline is reached. Elsewhere the speedhold mode is maintained. In this paper we will consider the problem of calculating the key switching points and will present a new formulation of the problem and a corresponding numerical solution procedure that may also provide a constructive proof that the optimal strategy is uniquely determined by the holding speed.

The rather encouraging feature of this work is that a small Australian firm, TMG International, has used algorithms developed by the Centre for Industrial and Applied Mathematics (CIAM) at UniSA to produce an on-board system, known as FreightMiser, that will provide advice to train drivers about optimal driving strategies. FreightMiser is currently being trialled by Pacific National for use on Australian freight trains. In-service trials indicate that fuel savings of more than 10% will be possible. This will amount to many millions of dollars per annum. This is joint work with Peter Pudney and Xuan Vu.

Dan Gordon, Constraint Technologies International
Airline optimisation problems – an industry perspective
Presentation (ppt)

Constraint Technologies International provides planning and decision support systems to industry. In particular, we are very active in the airline industry. Running a large airline requires the ability to solve large and challenging optimisation problems. One of the greatest challenges arises from the tension between the need for powerful numerical optimisation on the one hand, and the requirement to take into account a complicated and changing set of business rules on the other hand. Another major challenge is the problem of real-time disruption management, which requires a high degree of integration between problems that would normally be dealt with in a disjoint fashion.

We begin by describing some of the classic problems and challenges faced by airlines, such as crew scheduling and rostering, and some standard ways of soving these problems. We then go on to discuss particulars of CTI's approach to these problems; some of the practical challenges that we have been faced with and our response to these challenges; and our strategies for achieving flexibility without adversely affecting performance.

We finish by discussing some of the hardest (and indeed only partially solved) challenges of airline decision support: the increasing need for effective disruption management tools, and more generally the desirability of achieving greater integration between the various stages of planning and operations.

Jonathan Manton, Australian National University
On the theory of optimisation on manifolds

The traditional approach to developing a numerical algorithm for minimising a cost function on a manifold is to endow the manifold with a metric structure and then use judiciously parallel transport and the Riemannian exponential map to generalise Euclidean algorithms (such as the conjugate gradient or the Newton method) to the manifold setting. This talk will present a more general framework and explain the theoretical and practical advantages of working in this greater generality. Specifically, this framework enables an arbitrary algorithm in Euclidean space to be transported to a manifold, and is sufficiently flexible to allow the domains of attraction, the computational complexity and the convergence rate to be altered significantly, while at the same time guaranteeing the local convergence of the transported algorithm to be the same or faster than the original algorithm.

John Butcher, University of Auckland
Towards efficient general linear methods
Presentation (pdf)

The central aim in algorithms to solve initial value problems is the minimization of error+T*cost, where the Lagrange multiplier T is a user-supplied tolerance. We will discuss the derivation of a class of general linear methods which provide the ingredients for this optimisation to be carried out dynamically and which have a number of additional desirable properties.

Helmut Podhaisky, Martin Luther Universität
On variable order for general linear methods in Nordsieck form

Some general linear methods such as peer methods or methods with inherent Runge-Kutta stability (IRKS) are an attractive alternative to well established Runge-Kutta and multistep methods. They combine L-stability and high order and stage order at relatively low cost (i.e. with a diagonally implicit scheme) and allow an implementation as a linearly implicit method (with one Newton iteration). Promising methods up to order p=4 have been developed and tested.

Monitoring the error recursion leads to accurate error estimators. While integrating with a method of order p it is possible not only to estimate the truncation error of this method but also the truncation error of the method of order p+1 asymptotically correct. Numerical results for a variable stepsize and variable order implementation for stiff ODEs are given.

Dirk Nuyens, Katholieke Universiteit Leuven, Belgium
On the contruction of rank-1 lattice rules
Presentation (pdf)

In the last decade, so-called quasi-Monte Carlo rules have gained a lot of interest for the numerical evaluation of high dimensional integrals. These are equal-weight cubature rules where the point set Pn is a low-discrepancy point set over [0,1)s. A popular construction algorithm for one class of low-discrepancy point sets, namely rank-1 lattices, is the component-by-component (CBC) construction algorithm due to Ian H.~Sloan et al. The main benefits of this construction algorithm are its extensibility in dimension s, the usage of a weighted reproducing kernel Hilbert space and(strong) tractability in s under appropriate conditions on the weights.

Recently this algorithm was reformulated into a fast variant (for a prime number of points n) by exploiting its inherent structure which is represented explicitly by a linear algebra formulation, N. & Cools (2004). The algorithm hence shows its inner beauty in a repeated matrix-vector multiplication operating on the spectral properties of the error representer in the sequence of higher and higher dimensional function spaces. In N. & Cools (2005) the multiplicative structure for non-prime n was studied, giving a more complete view on the problem at hand.

From the last work we observe embedding properties of the matrix in the algorithm and in the projections of the point set. This can be put to use in creating lattice rules which are “extensible” in the number of points, which is a defect of current lattice rules (and constructions). These lattice rules are not so much “extensible”, but are actually embedded rules, where the smaller ones can be extended up to the maximum number of points. This is work in progress with Ronald Cools, Frances Y.~Kuo and Ian H.~Sloan.

Mike Osborne, Australian National University
Multiple Shooting Revisited
Presentation (pdf)

This talk will consider some stability questions relating to the use of a boundary value problem solver possessing a family relationship with multiple shooting for solving ODE estimation problems when the basic equations possess a strong dichotomy. Questions relating to the provision of compatible boundary conditions will be considered and evidence presented that stiffly stable numerical procedures may be possible.

Gene Golub, Stanford University
Matrices, moments and quadrature with emphasis on least squares problems
Presentation (pdf)

Given an n×n symmetric, positive definite matrix A and a real vector u, it is desirable to estimate and bound the quadratic form uTF(A)u/uTu where F is a differentiable function. This problem arises in estimating errors of linear systems, computing a parameter in a least squares problem with a quadratic constraint and determining bounds for the perturbations in least squares problems.

We describe a method based on the theory of moments and numerical quadrature for estimating the quadratic form. A basic tool is the Lanczos algorithm which can be used for computing the recursive relationship for orthogonal polynomials.

Georg Bock, University of Heidelberg
From parameter estimation to the design of optimal experiments in differential algebraic equations

The quantitative validation of nonlinear differential algebraic and partial differential equations models is a difficult task that requires mathematical methods for parameter estimation, numerical and statistical sensitivity analysis, and the optimal design of experiments. In the talk, we discuss the advantages of the multiple shooting discretization – or rather parametrization – for parameter estimation, and treatment of the discretized BVP as a specially structured, nonlinear equality constraint by Gauss-Newton methods. Based on an assessment of the statistical accuracy of the parameter estimates in terms of confidence regions, efficient optimization methods for the determination of experiments are presented which maximize the information gain subject to feasibility and cost constraints.

Special emphasis is put on questions of robustness: both of the parameter estimation wrt measurement errors – where polyhedral norms are important, and of optimal experiment designs wrt the uncertainty in the parameters. Applications from satellite dynamics, chemical and biochemical reactions will be given. As an outlook, the incorporation of these concepts into real estimation and control schemes are discussed.

Åke Björck, Linköping University
Recent developments in least squares computations — Bidiagonal decomposition and least squares
Presentation (pdf)

Any rectangular matrix A can be decomposed as A = UBVT, where U and V are orthogonal and B upper (or lower) bidiagonal. In a seminal paper from 1965, Golub and Kahan gave two different, but mathematically equivalent, algorithms for computing this decomposition. The first algorithm uses Householder transformations applied, alternately from left and right, to A. The second algorithm uses two two-term recurrence relations related to a Lanczos process. In this, only matrix-vector products with A and AT are needed.

The Householder algorithm led to the first stable algorithm for computing the SVD of A. This plays a crucial role in direct methods for determining the numerical rank of a matrix A and for solving possibly ill-conditioned, least squares problems. The Lanczos bidiagonalization algorithm forms the core of the Krylov subspace algorithm LSQR. This is the method of choice for solving large sparse least squares problems.

In this talk we review several more recent developments of both algorithms for reduction to bidiagonal form. Paige and Strakoš have shown that reduction to upper bidiagonal form of the matrix (b A) provides an elegant way to extract a minimally dimensioned core problem, not only for the linear least squares but also for the (weighted) total least squares problems.

For some time it has been known that, for a class of least squares problems derived from an ill-posed problem, the truncated bidiagonal reduction is often superior to the use of truncated SVD. This is related to the popular method of Partial Least Squares (PLS), used in statistics for handling colinearities. In some applications, the familiar loss of orthogonality that occurs in a Lanczos-type bidiagonalization algorithm is not acceptable. This can be cured by applying a (selective) Gram-Schmidt reorthogonalization. We discuss some possible alternatives to this. Several modifications to Householder bidiagonalization have also been suggested in order to improve complexity and accuracy in algorithms for computing the SVD.

Gordon Smyth, Walter & Eliza Hall Institute of Medical Research
Random effects, heteroscedastic regression, algorithms and microarray data analysis
Presentation (pdf)

Microarrays are one of the technologies underlying the genomic revolution in modern biology. They provide a cheap and accessible way to measure the activity level of genes in an organism, providing measurements for tens of thousands of genes simultaneously. All the usual statistical considerations for designed experiments apply when analysing microarray experiments, but in a massively parallel fashion. Traditional two-colour microarrays especially embody some of the traditional statistical notions of blocking, random effects and non-constant variance. The high throughput nature of the data gives rise to new issues both statistically and computationally. The fact that microarray software is typically used by non-mathematical scientists emphasises the need for reliability and efficiency of the computational algorithms. This talk will describe an approach to microarray data analysis which takes advantage of algorithmic developments for random effects and heteroscedastic regression models. The approach has been implemented in the R computing environment where it takes advantage of the R support for data structures, classes and elementary object-oriented programming.

Richard Brent, Australian National University
Challenges in solving large sparse linear systems over finite fields
Presentation (pdf)

This talk outlines how very large, sparse linear systems arise in the solution of problems of interest in computational number theory and public-key cryptography, such as the integer factorization and discrete logarithm problems. The linear systems are over finite fields, often the field GF(2) of two elements. We describe some algorithms for solving large sparse linear systems over GF(2), and compare them with algorithms for the real field. In particular, some “iterative” algorithms which are well-known to numerical analysts, such as the conjugate gradient and Lanczos algorithms, can be adapted to work over GF(2), but there are significant differences between algorithms for the real field and for GF(2).

William McLean, University of New South Wales
The conditioning of boundary element equations
Presentation (pdf)

We consider Galerkin discretization with standard nodal basis functions of first-kind boundary integral equations. It is well known that in such cases the condition number of the stiffness matrix grows as the mesh is refined. We begin by explaining this phenomenon in the familiar case when the mesh is quasi-uniform and then describe some results of joint work with Thanh Tran and Mark Ainsworth on the effect of local refinement. The final part of the talk discusses recent work with Ivan Graham in which the local refinement may be anisotropic, leading to meshes that fail to be shape-regular. This type of refinement arises naturally when one seeks to resolve edge singularities in 3D boundary element methods. We also see how a simple diagonal preconditioning can significantly reduce the condition number.


Nick Stokes, CSIRO Mathematical and Information Sciences
Material flows — matching micro-scale computations with macroscopic effects

For complex media in motion which is not governed by a partial differential equation with fixed boundaries, progress can be made by numerically modelling the interactions of smaller parts with more predictable behaviour. Examples include the various kinds of meshless fluid modelling, and the modelling of granular flows as individual collisions of rigid bodies. The direct numerical simulation of high Reynolds number flows is another. A challenge is that the micro-scale dictates a necessary resolution, which may be very different to the scale of the phenomenon being modelled. It is, for example, feasible to model a mill full of rocks, but not one full of sand, and yet this may represent the progress of the milling process.

What is badly needed is a kind of intermediate scale, in which the micro-scale modelling determines a larger scale quasi-continuum behaviour, at which scale further modelling can be performed.

This talk will unfortunately not fulfil this need, but will give some indications of the advances that are being made at the micro level.


Ian Sloan, University of New South Wales
The curse of dimensionality

Richard Bellman coined the phrase “the curse of dimensionality” to describe the extraordinarily rapid increase in the difficulty of most problems as the number of variables increases. A typical problem is numerical multiple integration, where the cost of any integration formula of product type obviously rises exponentially with the number of variables. Nevertheless, problems with hundreds or even thousands of variables do arise, and are now being tackled successfully. In this talk I will describe recent strategies, mathematical settings and constructions with which suitable integration problems (from mathematical finance, for example) are being successfully handled.

Markus Hegland, Australian National University
Resolution enhancement of spectra using differentiation
Presentation (pdf)

The analysis of measured spectra, when it reduces to finding the positions and heights of spectral lines, becomes an extremely difficult problem when the lines are closely spaced and broadened to the point where overlapping occurs, and explicit models for the shape of the peaks are unknown. In this talk, I will cover some recent results we obtained by applying a Hilbert space framework to analyse several methods and I will present new enhancement techniques based on numerical differentiation.

This is joint work with Bob Anderssen, of CSIRO Mathematical and Information Sciences.

Reference:
Resolution enhancement of spectra using differentiation **Hegland, Markus; Anderssen, Robert S. **Inverse Problems <http://www.ingentaconnect.com/content/iop/ip>, June 2005, vol. 21, no. 3, pp. 915-934(20)*