Past Computational Mathematics and Applications Seminar

The advent of multicore processors represents a disruptive event in the history of computer science as conventional parallel programming paradigms are proving incapable of fully exploiting their potential for concurrent computations. The need for different or new programming models clearly arises from recent studies which identify fine-granularity and dynamic execution as the keys to achieve high efficiency on multicore systems. This talk shows how these models can be effectively applied to the multifrontal method for the QR factorization of sparse matrices providing a very high efficiency achieved through a fine-grained partitioning of data and a dynamic scheduling of computational tasks relying on a dataflow parallel programming model. Moreover, preliminary results will be discussed showing how the multifrontal QR factorization can be accelerated by using low-rank approximation techniques.

The talk will survey recent developments concerning the existence and the approximation of global weak solutions to a general class of coupled microscopic-macroscopic bead-spring chain models that arise in the kinetic theory of dilute solutions of polymeric liquids with noninteracting polymer chains. The class of models involves the unsteady incompressible Navier-Stokes equations in a bounded domain for the velocity and the pressure of the fluid, with an elastic extra-stress tensor appearing on the right-hand side of the momentum equation. The extra-stress tensor stems from the random movement of the polymer chains and is defined by the Kramers expression through the associated probability density function that satisfies a Fokker-Planck type parabolic equation. Models of this kind were proposed in work of Hans Kramers in the early 1940's, and the existence of global weak solutions to the model has been a long-standing question in the mathematical analysis of kinetic models of dilute polymers.

\\

\\

We also discuss computational challenges associated with the numerical approximation of the high-dimensional Fokker-Planck equation featuring in the model.

In this talk I review the use of the spectral decomposition for understanding the solution of ill-posed inverse problems. It is immediate to see that regularization is needed in order to find stable solutions. These solutions, however, do not typically allow reconstruction of signal features such as edges. Generalized regularization assists but is still insufficient and methods of total variation are commonly suggested as an alternative. In the talk I consider application of standard approaches from Tikhonov regularization for finding appropriate regularization parameters in the total variation augmented Lagrangian implementations. Areas for future research will be considered.

In this talk we present an overview of some recent developments concerning the a posteriori error analysis and adaptive mesh design of $h$- and $hp$-version discontinuous Galerkin finite element methods for the numerical approximation of second-order quasilinear elliptic boundary value problems. In particular, we consider the derivation of computable bounds on the error measured in terms of an appropriate (mesh-dependent) energy norm in the case when a two-grid approximation is employed. In this setting, the fully nonlinear problem is first computed on a coarse finite element space $V_{H,P}$. The resulting 'coarse' numerical solution is then exploited to provide the necessary data needed to linearise the underlying discretization on the finer space $V_{h,p}$; thereby, only a linear system of equations is solved on the richer space $V_{h,p}$. Here, an adaptive $hp$-refinement algorithm is proposed which automatically selects the local mesh size and local polynomial degrees on both the coarse and fine spaces $V_{H,P}$ and $V_{h,p}$, respectively. Numerical experiments confirming the reliability and efficiency of the proposed mesh refinement algorithm are presented.

Standard numerical schemes for acoustic scattering problems suffer from the restriction that the number of degrees of freedom required to achieve a prescribed level of accuracy must grow at least linearly with respect to frequency in order to maintain accuracy as frequency increases. In this talk, we review recent progress on the development and analysis of hybrid numerical-asymptotic boundary integral equation methods for these problems. The key idea of this approach is to form an ansatz for the solution based on knowledge of the high frequency asymptotics, allowing one to achieve any required accuracy via the approximation of only (in many cases provably) non-oscillatory functions. In particular, we discuss very recent work extending these ideas for the first time to non-convex scatterers.

Computing the eigenvalue decomposition of a symmetric matrix and the singular value decomposition of a general matrix are two of the central tasks in numerical linear algebra. There has been much recent work in the development of linear algebra algorithms that minimize communication cost. However, the reduction in communication cost sometimes comes at the expense of significantly more arithmetic and potential instability.

\\

\\

In this talk I will describe algorithms for the two decompositions that have optimal communication cost and arithmetic cost within a small factor of those for the best known algorithms. The key idea is to use the best rational approximation of the sign function, which lets the algorithm converge in just two steps. The algorithms are backward stable and easily parallelizable. Preliminary numerical experiments demonstrate their efficiency.

We show that the steepest-descent and Newton's methods for unconstrained nonconvex optimization

under standard assumptions may both require a number of iterations and function evaluations

arbitrarily close to the steepest-descent's global worst-case complexity bound. This implies that

the latter upper bound is essentially tight for steepest descent and that Newton's method may be as

slow as the steepest-descent method in the worst case. Then the cubic regularization of Newton's

method (Griewank (1981), Nesterov & Polyak (2006)) is considered and extended to large-scale

problems, while preserving the same order of its improved worst-case complexity (by comparison to

that of steepest-descent); this improved worst-case bound is also shown to be tight. We further

show that the cubic regularization approach is, in fact, optimal from a worst-case complexity point

of view amongst a wide class of second-order methods. The worst-case problem-evaluation complexity

of constrained optimization will also be discussed. This is joint work with Nick Gould (Rutherford

Appleton Laboratory, UK) and Philippe Toint (University of Namur, Belgium).

We present progress on an Interior Point based multi-step solution approach for stochastic programming problems. Our approach works with a series of scenario trees that can be seen as successively more accurate discretizations of an underlying probability distribution and employs IPM warmstarts to "lift" approximate solutions from one tree to the next larger tree.

Nonnormality is a well studied subject in the context of partial differential operators. Yet, only little is known for boundary integral operators. The only well studied case is the unit ball, where the standard single layer, double layer and conjugate double layer potential operators in acoustic scattering diagonalise in a unitary basis. In this talk we present recent results for the analysis of spectral decompositions and nonnormality of boundary integral operators on more general domains. One particular application is the analysis of stability constants for boundary element discretisations. We demonstrate how these are effected by nonnormality and give several numerical examples, illustrating these issues on various domains.

Many iterative algorithms for large sparse matrix problems are based on orthogonality (or $A$-orthogonality, bi-orthogonality, etc.), but these properties can be lost very rapidly using vector orthogonalization (subtracting multiples of earlier supposedly orthogonal vectors from the latest vector to produce the next orthogonal vector). Yet many of these algorithms are some of the best we have for very large sparse problems, such as Conjugate Gradients, Lanczos' method for the eigenproblem, Golub and Kahan bidiagonalization, and MGS-GMRES.

\\

\\

Here we describe an ideal form of orthogonal matrix that arises from any sequence of supposedly orthogonal vectors. We illustrate some of its fascinating properties, including a beautiful measure of orthogonality of the original set of vectors. We will indicate how the ideal orthogonal matrix leads to expressions for new concepts of stability of such iterative algorithm. These are expansions of the concept of backward stability for matrix transformation algorithms that was so effectively developed and applied by J. H. Wilkinson (FRS). The resulting new expressions can be used to understand the subtle and effective performance of some (and hopefully eventually all) of these iterative algorithms.

The fundamental task in climate variability research is to eke

out structure from climate signals. Ideally we want a causal

connection between a physical process and the structure of the

signal. Sometimes we have to settle for a correlation between

these. The challenge is that the data is often poorly

constrained and/or sparse. Even though many data gathering

campaigns are taking place or are being planned, the very high

dimensional state space of the system makes the prospects of

climate variability analysis from data alone impractical.

Progress in the analysis is possible by the use of models and

data. Data assimilation is one such strategy. In this talk we

will describe the methodology, illustrate some of its

challenges, and highlight some of the ways our group has

proposed to improving the methodology.

The liquid crystal (LC) flow model is a coupling between

orientation (director field) of LC molecules and a flow field.

The model may probably be one of simplest complex fluids and

is very similar to a Allen-Cahn phase field model for

multiphase flows if the orientation variable is replaced by a

phase function. There are a few large or small parameters

involved in the model (e.g. the small penalty parameter for

the unit length LC molecule or the small phase-change

parameter, possibly large Reynolds number of the flow field,

etc.). We propose a C^0 finite element formulation in space

and a modified midpoint scheme in time which accurately

preserves the inherent energy law of the model. We use C^0

elements because they are simpler than existing C^1 element

and mixed element methods. We emphasise the energy law

preservation because from the PDE analysis point of view the

energy law is very important to correctly catch the evolution

of singularities in the LC molecule orientation. In addition

we will see numerical examples that the energy law preserving

scheme performs better under some choices of parameters. We

shall apply the same idea to a Cahn-Hilliard phase field model

where the biharmonic operator is decomposed into two Laplacian

operators. But we find that under our scheme non-physical

oscillation near the interface occurs. We figure out the

reason from the viewpoint of differential algebraic equations

and then remove the non-physical oscillation by doing only one

step of a modified backward Euler scheme at the initial time.

A number of numerical examples demonstrate the good

performance of the method. At the end of the talk we will show

how to apply the method to compute a superconductivity model,

especially at the regime of Hc2 or beyond. The talk is based

on a few joint papers with Chun Liu, Qi Wang, Xingbin Pan and

Roland Glowinski, etc.

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.