Past Computational Mathematics and Applications Seminar

Variational data assimilation techniques for optimal state estimation in very large environmental systems currently use approximate Gauss-Newton (GN) methods. The GN method solves a sequence of linear least squares problems subject to linearized system constraints. For very large systems, low resolution linear approximations to the model dynamics are used to improve the efficiency of the algorithm. We propose a new approach for deriving low order system approximations based on model reduction techniques from control theory which can be applied to unstable stochastic systems. We show how this technique can be combined with the GN method to retain the response of the dynamical system more accurately and improve the performance of the approximate GN method.

Computational Fluid Dynamics (CFD) has become an

indispensable tool in designing turbomachinery components in all sectors of

Rolls-Royce business units namely, Aerospace, Industrial, Marine and Nuclear.

Increasingly sophisticated search and optimisation techniques are used based on

both traditional optimisation methods as well as, design of computer experiment

techniques, advanced surrogate methods, and evolutionary optimisation

techniques. Geometry and data representation as well as access, queuing and

loading control of large high performance computing clusters are areas of

research to establish the most efficient techniques for improving the

performance of an already highly efficient modern jet engine.

\\

\\

This presentation focuses on a high fidelity design

optimisation framework called SOPHY that is used in Rolls-Royce to provide

parametric geometry, automatic meshing, advanced design-space search

algorithms, accurate and robust CFD methodology and post-processing. The

significance of including the so-called real geometry features and interaction

of turbomachinery components in the optimisation cycle are discussed. Examples are drawn from real world

applications of the SOPHY design systems in an engine project.

We will discuss the use of hypergraph-based methods for orderings of sparse matrices in Cholesky, LU and QR factorizations. For the Cholesky factorization case, we will investigate a recent result on pattern-wise decomposition of sparse matrices, generalize the result and develop algorithmic tools to obtain effective ordering methods. We will also see that the generalized results help us formulate the ordering problem in LU much like we do for the Cholesky case, without ever symmetrizing the given matrix $A$ as $A+A^{T}$ or $A^{T}A$. For the QR factorization case, the use of hypergraph models is fairly standard. We will nonetheless highlight the fact that the method again does not form the possibly much denser matrix $A^{T}A$. We will see comparisons of the hypergraph-based methods with the most common alternatives in all three cases.

\\

\\

This is joint work with Iain S. Duff.

We will study the question of whether the adjoint of a given matrix can be written as a rational function in the matrix. After showing necessary and sufficient conditions, rational interpolation theory will help to characterize the most important existing cases. Several topics related to our question will be explored. They range from short recurrence Krylov subspace methods to the roots of harmonic polynomials and harmonic rational functions. The latter have recently found interesting applications in astrophysics, which will briefly be discussed as well.

We develop a randomized block-coordinate descent method for minimizing the sum of a smooth and a simple nonsmooth block-separable convex function and prove that it obtains an $\epsilon$-accurate solution with probability at least $1-\rho$ in at most $O[(2n/\epsilon)\log(1/\rho)]$ iterations, where $n$ is the dimension of the problem. For strongly convex functions the method converges linearly. This extends recent results of Nesterov [Efficiency of coordinate descent methods on huge-scale optimization problems, CORE Discussion Paper \#2010/2], which cover the smooth case, to composite minimization, while at the same time improving the complexity by the factor of 4 and removing $\epsilon$ from the logarithm. More importantly, in contrast with the aforementioned work in which the authors achieve the results by applying their method to a regularized version of the objective function with an unknown scaling factor, we show that this is not necessary, thus achieving true iteration complexity bounds. In the smooth case we also allow for arbitrary probability vectors and non-Euclidean norms. Time permitting, we will also mention new iteration complexity results for a parallel version of the method.

\\

\\

In the second part of the talk we demonstrate numerically that the algorithm is able to solve huge-scale $\ell_1$-regularized support vector machine and least squares problems with billion variables. Finally, we present preliminary computational results with a GPU-accelerated parallel version of the method, on truss topology design instances, achieving speedups of up to two orders of magnitude when compared to a single-core implementation in C.

\\

\\

References:

\\

P. Richtarik and M. Takac, Iteration complexity of randomized block-coordinate descent methods for minimizing a composite function, 2011; http://arxiv.org/abs/1107.2848

\\

P. Richtarik and M. Takac, Efficient serial and parallel coordinate descent methods for huge-scale truss topology design, 2011; http://www.optimization-online.org/DB_HTML/2011/08/3118.html

\\

P. Richtarik and M. Takac, Efficiency of randomized coordinate descent methods on minimization problems with a composite objective function, Proceedings of SPARS11

Liquid impact is a key issue in various industrial applications (seawalls, offshore structures, breakwaters, sloshing in tanks of liquefied natural gas vessels, wave energy converters, offshore wind turbines, etc). Numerical simulations dealing with these applications have been performed by many groups, using various types of numerical methods. In terms of the numerical results, the outcome is often impressive, but the question remains of how relevant these results are when it comes to determining impact pressures. The numerical models are too simplified to reproduce the high variability of the measured pressures. In fact, for the time being, it is not possible to simulate accurately both global and local effects. Unfortunately it appears that local effects predominate over global effects when the behaviour of pressures is considered.

\\

\\

Having said this, it is important to point out that numerical studies can be quite useful to perform sensitivity analyses in idealized conditions such as a liquid mass falling under gravity on top of a horizontal wall and then spreading along the lateral sides. Simple analytical models inspired by numerical results on idealized problems can also be useful to predict trends.

\\

\\

The talk is organized as follows: After an introduction on some of the industrial applications, it will be explained to what extent numerical studies can be used to improve our understanding of impact pressures. Results on a liquid mass hitting a wall obtained by various numerical codes will be shown.

In this talk, I will discuss various aspects of approximation by radial basis functions on spheres. After a short introduction to the subject of scattered data approximation on spheres and optimal recovery, I will particularly talk about error analysis, a hybrid approximation scheme involving polynomials and radial basis functions and, if time permits, solving nonlinear parabolic equations on spheres.

there will be no seminar in this week.

Chebyshev polynomials are arguably the most useful orthogonal polynomials for computational purposes. In one dimension they arise from the close relationship that exists between Fourier series and polynomials. We describe how this relationship generalizes to Fourier series on certain symmetric lattices, that exist in all dimensions. The associated polynomials can not be seen as tensor-product generalizations of the one-dimensional case. Yet, they still enjoy excellent properties for interpolation, integration, and spectral approximation in general, with fast FFT-based algorithms, on a variety of domains. The first interesting case is the equilateral triangle in two dimensions (almost). We further describe the generalization of Chebyshev polynomials of the second kind, and many new kinds are found when the theory is completed. Connections are made to Laplacian eigenfunctions, representation theory of finite groups, and the Gibbs phenomenon in higher dimensions.

Following the framework of the heterogeneous multiscale method, we present a numerical method for nonlinear elliptic homogenization problems. We briefly review the numerical, relying on an efficient coupling of macro and micro solvers, for linear problems. A fully discrete analysis is then given for nonlinear (nonmonotone) problems, optimal convergence rates in the H¹ and L² norms are derived and the uniqueness of the method is shown on sufficiently fine macro and micro meshes.

Numerical examples confirm the theoretical convergence rates and illustrate the performance and versatility of our approach.

This talk is about the Induced Dimension Reduction (IDR) methods developed by Peter Sonneveld and, more recently, Martin van Gijzen. We sketch the history, outline the underlying principle, and give a few details about different points of view on this class of Krylov subspace methods. If time permits, we briefly outline some recent developments in this field and the benefits and drawbacks of these and IDR methods in general.

Many continuum models have been derived in recent years which describe the self-assembly of industrially utilisable crystalline films to a level of detail that allows qualitative comparisons with experiments. For thin-film problems, where the characteristic length scales in vertical and horizontal directions differ significantly, the governing surface diffusion equations can be reduced to simpler PDEs by making use of asymptotic expansions. Many mathematical problems and solutions emerge from such new evolution equations and many of them remind of Cahn-Hilliard type equations. The surface diffusion models are of high, of fourth or even sixth, order.

We present the modeling, model reduction and simulation results for heteroepitaxial growth as for Ge/Si quantum dot self-assembly. The numerical methods we are using are based on trigonometric interpolation. These kind of pseudospectral methods seem very well suited for simulating the coarsening of large quantum dot arrays. When the anisotropy of the growing crystalline film is strong, it might become necessary to add a corner regularisation to the model. Then the transition region between neighboring facets is still smooth, but its scale is rather small. In this case it might be useful to think about an adaptive extension of the existing method.

Figure 1: Ostwald ripening process of quantum dots depicted at consecutive time points. One fourth of the whole, periodic, simulated domain is shown.

Joint work with Peter Evans and Barbara Wagner

This seminar will be held at the Rutherford Appleton Laboratory near Didcot.

Abstract:

Numerical calculations of laminar flow in a two-dimensional channel with a sudden expansion exhibit a symmetry-breaking bifurcation at Reynolds number 40.45 when the expansion ratio is 3:1. In the experiments reported by Fearn, Mullin and Cliffe [1] there is a large perturbation to this bifurcation and the agreement with the numerical calculations is surprisingly poor. Possible reasons for this discrepancy are explored using modern techniques for uncertainty quantification.

When experimental equipment is constructed there are, inevitably, small manufacturing imperfections that can break the symmetry in the apparatus. In this work we considered a simple model for these imperfections. It was assumed that the inlet section of the channel was displaced by a small amount and that the centre line of the inlet section was not parallel to the centre line of the outlet section. Both imperfections were modelled as normal random variables with variance equal to the manufacturing tolerance. Thus the problem to be solved is the Navier-Stokes equations in a geometry with small random perturbations. A co-ordinate transformation technique was used to transform the problem to a fixed deterministic domain but with random coefficient appearing in the transformed Navier-Stokes equations. The resulting equations were solved using a stochastic collocation technique that took into account the fact that the problem has a discontinuity in parameter space arising from the bifurcation structure in the problem.

The numerical results are in the form of an approximation to a probability measure on the set of bifurcation diagrams. The experimental data of Fearn, Mullin and Cliffe are consistent with the computed solutions, so it appears that a satisfactory explanation for the large perturbation can be provided by manufacturing imperfections in the experimental apparatus.

The work demonstrates that modern methods for uncertainty quantification can be applied successfully to a bifurcation problem arising in fluid mechanics. It should be possible to apply similar techniques to a wide range of bifurcation problems in fluid mechanics in the future.

References:

[1] R M Fearn, T Mullin and K A Cliffe Nonlinear flow phenomena in a symmetric sudden expansion, J. Fluid Mech. 211, 595-608, 1990.

Please note that this is a short notice change from the originally advertised talk by Dr Shahrokh Shahpar (Rolls-Royce plc.)

The new talk "Multilevel Monte Carlo method" is given by Mike Giles, Oxford-Man Institute of Quantitative Finance, Mathematical Institute, University of Oxford.

Joint work with Rob Scheichl, Aretha Teckentrup (Bath) and Andrew Cliffe (Nottingham)

Over the last several decades, many mesh generation methods and a plethora of adaptive methods for solving differential equations have been developed.  In this talk, we take a general approach for describing the mesh generation problem, which can be considered as being in some sense equivalent to determining a coordinate transformation between physical space and a computational space.  Our description provides some new theoretical insights into precisely what is accomplished from mesh equidistribution (which is a standard adaptivity tool used in practice) and mesh alignment.  We show how variational mesh generation algorithms, which have historically been the most common and important ones, can generally be compared using these mesh generation principles.  Lastly, we relate these to a variety of moving mesh methods for solving time-dependent PDEs.

This is joint work with Weizhang Huang, Kansas University

In this talk we discuss the design of efficient numerical methods for solving symmetric indefinite linear systems arising from mixed approximation of elliptic PDE problems with associated constraints. Examples include linear elasticity (Navier-Lame equations), steady fluid flow (Stokes' equations) and electromagnetism (Maxwell's equations).

The novel feature of our iterative solution approach is the incorporation of error control in the natural "energy" norm in combination with an a posteriori estimator for the PDE approximation error. This leads to a robust and optimally efficient stopping criterion: the iteration is terminated as soon as the algebraic error is insignificant compared to the approximation error. We describe a "proof of concept" MATLAB implementation of this algorithm, which we call EST_MINRES, and we illustrate its effectiveness when integrated into our Incompressible Flow Iterative Solution Software (IFISS) package (http://www.manchester.ac.uk/ifiss/).

In this talk, we focus on the analytical properties of a decentralized auction, namely the PAUSE Auction Procedure. We prove that the revenue of the auctioneer from PAUSE is greater than or equal to the profit from the well-known VCG auction when there are only two bidders and provide lower bounds on the profit for arbitrary number of bidders. Based on these bounds and observations from auctions with few items, we propose a modification of the procedure that increases the profit. We believe that this study, which is still in progress, will be a milestone in designing better decentralized auctions since it is the first analytical study on such auctions with promising results.

Polynomial Programs are ussually solved by using hierarchies of convex relaxations. This scheme rapidly becomes computationally expensive and is often tractable only for problems of small sizes. We propose an iterative scheme that improves an initial relaxation without incurring exponential growth in size. The key ingredient is a dynamic scheme for generating valid polynomial inequalities for general polynomial programs. These valid inequalities are then used to construct better approximations of the original problem. As a result, the proposed scheme is in principle scalable to large general combinatorial optimization problems.

Joint work with Bissan Ghaddar and Miguel Anjos

Based on an MPI library written over 10 years ago, OP2 is a new open-source library which is aimed at application developers using unstructured grids. Using a single API, it targets a variety of HPC architectures, including both manycore GPUs and multicore CPUs with vector units. The talk will cover the API design, key aspects of the parallel implementation on the different platforms, and preliminary performance results on a small but representative CFD test code.

Project homepage: http://people.maths.ox.ac.uk/gilesm/op2/

We consider the iterative solution of large sparse linear least squares (LS) problems. Specifically, we focus on the design and implementation of reliable stopping criteria for the widely-used algorithm LSQR of Paige and Saunders. First we perform a backward perturbation analysis of the LS problem. We show why certain projections of the residual vector are good measures of convergence, and we propose stopping criteria that use these quantities. These projections are too expensive to compute to be used directly in practice. We show how to estimate them efficiently at every iteration of the algorithm LSQR. Our proposed stopping criteria can therefore be used in practice.

This talk is based on joint work with Xiao-Wen Chang, Chris Paige, Pavel Jiranek, and Serge Gratton.

In various fields of application, one must solve very large scale boundary value problems using parallel solvers and supercomputers. The domain decomposition approach partitions the large computational domain into smaller computational subdomains. In order to speed up the convergence, we have several ``optimized'' algorithm that use Robin transmission conditions across the artificial interfaces (FETI-2LM). It is known that this approach alone is not sufficient: as the number of subdomains increases, the number of iterations required for convergence also increases and hence the parallel speedup is lost. A known solution for classical Schwarz methods as well as FETI algorithms is to incorporate a ``coarse grid correction'', which is able to transmit low-frequency information more quickly across the whole domain. Such algorithms are known to ``scale weakly'' to large supercomputers. A coarse grid correction is also necessary for FETI-2LM methods. In this talk, we will introduce and analyze coarse grid correction algorithms for FETI-2LM domain decomposition methods.