![]() |
Mathematical Sciences Institute (MSI)
Events - Special Year 2005
|
|
2005 Special Year on Computational Mathematics
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 |
| 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.
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: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.
Nonparametric tests for distributionsLinear 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.
Regression tree ensemble modelsTree-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.
On homotopy algorithms in statisticsThe 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).
On singularly perturbed Markov chainsNothing 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.
Multiscale modelling of cellular processesThe 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.
Thin film models in a stochastic settingDunn & 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.
An application of sparse grids to the estimation of probability density functionsIn 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.
Importance sampling plays a crucial role for quasi-Monte Carlo methodsMany 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.
A QR-Based tridiagonalization algorithm for nonsymmetric matricesThe 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.
Domain decomposition relations as a source of fast algorithms for evaluating integral operatorsFor 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.
Distributing points on manifoldsThere 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.
The use of severely limited subsequence in set partitioning optimisationMany 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.
On the computational complexity of Jacobians and HessiansMost 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.
Evaluating optimizationWe 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?
Least-squares methods in optimizationLeast-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.
Doubly stochastic matrices, the Hamiltonian cycle problem & a control theoretic perspective of static optimization problems.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.
Calculation of optimal driving strategies for freight trainsThe 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.
Airline optimisation problems – an industry perspectiveConstraint 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.
On the theory of optimisation on manifoldsThe 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.
Towards efficient general linear methodsThe 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.
On variable order for general linear methods in Nordsieck formSome 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.
On the contruction of rank-1 lattice rulesIn 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.
Multiple Shooting RevisitedThis 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.
Matrices, moments and quadrature with emphasis on least squares problemsGiven 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.
From parameter estimation to the design of optimal experiments in differential algebraic equationsThe 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.
Recent developments in least squares computations — Bidiagonal decomposition and least squaresAny 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.
Random effects, heteroscedastic regression, algorithms and microarray data analysisMicroarrays 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.
Challenges in solving large sparse linear systems over finite fieldsThis 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).
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.
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.
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 UniversityThe 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)*
|
Page last updated: 3 November, 2006 Please direct all enquiries to: MSI webmaster Page authorised by: Director, MSI |
| The Australian National University - CRICOS Provider Number 00120C |