Past Computational Mathematics and Applications Seminar

Current community models in the geosciences employ a variety of numerical methods from finite-difference, finite-volume, finite- or spectral elements, to pseudospectral methods. All have specialized strengths but also serious weaknesses. The first three methods are generally considered low-order and can involve high algorithmic complexity (as in triangular elements or unstructured meshes). Global spectral methods do not practically allow for local mesh refinement and often involve cumbersome algebra. Radial basis functions have the advantage of being spectrally accurate for irregular node layouts in multi-dimensions with extreme algorithmic simplicity, and naturally permit local node refinement on arbitrary domains. We will show test examples ranging from vortex roll-ups, modeling idealized cyclogenesis, to the unsteady nonlinear flows posed by the shallow water equations to 3-D mantle convection in the earth’s interior. The results will be evaluated based on numerical accuracy, stability and computational performance.

This is a special talk outside the normal Computational Mathematics and Application seminar series. Please note it takes place on a Wednesday.

An iterative substructuring method with balancing Neumann-Neumann preconditioner is known as an efficient parallel algorithm for the elasticity equations. This method was extended to the Stokes equations by Pavarino and Widlund [2002]. In their extension, Q2/P0-discontinuous elements are used for velocity/pressure and a Schur complement system within "benign space", where incompressibility satisfied, is solved by CG method.

For the construction of the coarse space for the balancing preconditioner, some supplementary solvability conditions are considered. In our algorithm for 3-D computation, P1/P1 elements for velocity/pressure with pressure stabilization are used to save computational cost in the stiffness matrix. We introduce a simple coarse space similar to the one of elasticity equations. Owing to the stability term, solvabilities of local Dirichlet problem for a Schur complement system, of Neumann problem for the preconditioner, and of the coarse space problem are ensured. In our implementation, local Dirichlet and Neumann problems are solved by a 4x4-block modified Cholesky factorization procedure with an envelope method, which leads to fast computation with small memory requirement. Numerical result on parallel efficiency with a shared memory computer will be shown. Direct use of the Stokes solver in an application of Earth's mantle convection problem will be also shown.

Work with N.K. Nichols (Reading), C. Boess & A. Bunse-Gerstner (Bremen)

The Gauss-Newton (GN) method is a well known iterative technique for solving nonlinear least squares problems subject to dynamical system constraints. Such problems arise commonly from applications in optimal control and state estimation. Variational data assimilation systems for weather, ocean and climate prediction currently use approximate 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 method for deriving low order system approximations based on model reduction techniques from control theory. We show how this technique can be combined with the GN method to give a state estimation technique that retains more of the dynamical information of the full system. Numerical experiments using a shallow-water model illustrate the superior performance of model reduction to standard truncation techniques.

For the task of solving PDEs, finite difference (FD) methods are particularly easy to implement. Finite element (FE) methods are more flexible geometrically, but tend to be difficult to make very accurate. Pseudospectral (PS) methods can be seen as a limit of FD methods if one keeps on increasing their order of accuracy. They are extremely effective in many situations, but this strength comes at the price of very severe geometric restrictions. A more standard introduction to PS methods (rather than via FD methods of increasing orders of accuracy) is in terms of expansions in orthogonal functions (such as Fourier, Chebyshev, etc.).

Radial basis functions (RBFs) were first proposed around 1970 as a tool for interpolating scattered data. Since then, both our knowledge about them and their range of applications have grown tremendously. In the context of solving PDEs, we can see the RBF approach as a major generalization of PS methods, abandoning the orthogonality of the basis functions and in return obtaining much improved simplicity and flexibility. Spectral accuracy becomes now easily available also when using completely unstructured meshes, permitting local node refinements in critical areas. A very counterintuitive parameter range (making all the RBFs very flat) turns out to be of special interest. Computational cost and numerical stability were initially seen as serious difficulties, but major progress have recently been made also in these areas.

Quasicontinuum methods are a prototypical class of atomistic-to-continuum coupling methods. For example, we may wish to model a lattice defect (a vacancy or a dislocation) by an atomistic model, but the elastic far field by a continuum model. If the continuum model is consistent with the atomistic model (e.g., the Cauchy--Born model) then the main question is how the interface treatment affects the method.

In this talk I will introduce three of the main ideas how to treat the interface. I will explain their strengths and weaknesses by formulating the simplest possible non-trivial model problem and then simply analyzing the two classical concerns of numerical analysis: consistency and stability.

Inverse problems are often ill-posed, with solutions that depend sensitively on data. Regularization of some form is often used to counteract this. I will describe an approach to regularization, based on a Bayesian formulation of the problem, which leads to a notion of well-posedness for inverse problems, at the level of probability measures.

The stability which results from this well-posedness may be used as the basis for understanding approximation of inverse problems in finite dimensional spaces. I will describe a theory which carries out this program.

The ideas will be illustrated with the classical inverse problem for the heat equation, and then applied to so more complicated inverse problems arising in data assimilation, such as determining the initial condition for the Navier-Stokes equation from observations.

In recent years, the interdisciplinary field of imaging science has been experiencing an explosive growth in research activities including more models being developed, more publications generated, and above all wider applications attempted.
In this talk I shall first give an overview of the various imaging work carried out in our Liverpool group, some with collaborations with UCLA (T F Chan), CUHK (R H Chan) and Bergen (X C Tai) and several colleagues from other departments in Liverpool. Then I shall focus on two pieces of recent work, denoising and segmentation respectively:
(i) Image denoising has been a research topic deeply investigated within the last two decades. Even algorithmically the well-known ROF model (1992) can be solved efficiently. However less work has been done on models using high order regularization. I shall describe our first and successful attempt to develop a working multilevel algorithm for a 4th order nonlinear denoising model, and our work on solving the combined denoising and deblurring problem, different from the reformulation approach by M N Ng and W T Yin (2008) et al.
(ii) the image active contour model by Chan-Vese (2001) can be solved efficiently both by a geometric multigrid method and by an optimization based multilevel method. Surprisingly the new multilevel methods can find a solution closer to the global minimize than the existing unilevel methods. Also discussed are some recent work (jointly with N Badshah) on selective segmentation that has useful medical applications.

Geometric integration is the numerical integration of a differential equation, while preserving one or more of its geometric/physical properties exactly, i.e. to within round-off error.

Many of these geometric properties are of crucial importance in physical applications: preservation of energy, momentum, angular momentum, phase-space volume, symmetries, time-reversal symmetry, symplectic structure and dissipation are examples. The field has tantalizing connections to dynamical systems, as well as to Lie groups.

In this talk we first present a survey of geometric numerical integration methods for differential equations, and then exemplify this by discussing symplectic vs energy-preserving integrators for ODEs as well as for PDEs.

We discuss numerical methods for the solution of the palindromic eigenvalue problem Ax=λ A^Tx, where A is a complex matrix. Such eigenvalue problems occur, for example, in the vibration analysis of rail tracks.

The structure of palindromic eigenvalue problems leads to a symmetry in the spectrum: all eigenvalues occur in reciprocal pairs. The need for preservation of this symmetry in finite precision arithmetic requires the use of structure-preserving numerical methods. In this talk, we explain how such methods can be derived.

The aim of this talk is to render the power of (short-step) interior-point methods for linear programming (and by extension, convex programming) intuitively understandable to those who have a basic training in numerical methods for dynamical systems solving. The connection between the two areas is made by interpreting line-search methods in a forward Euler framework, and by analysing the algorithmic complexity in terms of the stiffness of the vector field of search directions. Our analysis cannot replicate the best complexity bounds, but due to its weak assumptions it also applies to inexactly computed search directions and has explanatory power for a wide class of algorithms.

Co-Author: Coralia Cartis, Edinburgh University School of Mathematics.

This talk is concerned with the probabilistic behaviour of a condition

number C(A) for the problem of deciding whether a system of n

homogeneous linear inequalities in m unknowns has a non-zero solution.

In the case where the input system is feasible, the exact probability

distribution of the condition number for random inputs is determined,

and a sharp bound for the general case. In particular, for the

expected value of the logarithm log C(A), an upper bound of order

O(log m) (m the number of variables) is presented which does not

depend on the number of inequalities.

The probability distribution of the condition number C(A) is closely

related to the probability of covering the m-sphere with n spherical

caps of a given radius. As a corollary, we obtain bounds on the

probability of covering the sphere with random caps.

The climate is largely determined by the ocean flow, which in itself is driven by wind and by gradients in temperature and salinity. Nowadays numerical models exist that are able to describe the occurring phenomena not only qualitatively but also quantitatively. At the Institute for Marine and Atmospheric research Utrecht (IMAU) a so-called thermohaline circulation model is developed in which methods of dynamical systems theory are used to study the stability of ocean flows. Here bifurcation diagrams are constructed by varying the strength of the forcing, for instance the amount of fresh water coming in from the north due to melting. For every value of the strength we have to solve a nonlinear system, which is handled by a Newton-type method. This produces many linear systems to be solved. 

In the talk the following will be addressed: the form of the system of equations, a special purpose method which uses Trilinos and MRILU. The latter is a multilevel ILU preconditioner developed at Groningen University. Results of the approach obtained on the Dutch national supercomputer will be shown.

For large--scale saddle point problems, the application of exact iterative schemes and preconditioners may be computationally expensive. In practical situations, only approximations to the inverses of the diagonal block or the related cross-product matrices are considered, giving rise to inexact versions of various solvers. Therefore, the approximation effects must be carefully studied. In this talk we study numerical behavior of several iterative Krylov subspace solvers applied to the solution of large-scale saddle point problems. Two main representatives of the segregated solution approach are analyzed: the Schur complement reduction method, based on an (iterative) elimination of primary variables and the null-space projection method which relies on a basis for the null-space for the constraints. We concentrate on the question what is the best accuracy we can get from inexact schemes solving either Schur complement system or the null-space projected system when implemented in finite precision arithmetic. The fact that the inner solution tolerance strongly influences the accuracy of computed iterates is known and was studied in several contexts.

In particular, for several mathematically equivalent implementations we study the influence of inexact solving the inner systems and estimate their maximum attainable accuracy. When considering the outer iteration process our rounding error analysis leads to results similar to ones which can be obtained assuming exact arithmetic. The situation is different when we look at the residuals in the original saddle point system. We can show that some implementations lead ultimately to residuals on the the roundoff unit level independently of the fact that the inner systems were solved inexactly on a much higher level than their level of limiting accuracy. Indeed, our results confirm that the generic and actually the cheapest implementations deliver the approximate solutions which satisfy either the second or the first block equation to the working accuracy. In addition, the schemes with a corrected direct substitution are also very attractive. We give a theoretical explanation for the behavior which was probably observed or it is already tacitly known. The implementations that we pointed out as optimal are actually those which are widely used and suggested in applications.

Multicore chips are nearly ubiquitous in modern machines, and to fully exploit this continuation of Moore's Law, numerical algorithms need to be able to exploit parallelism. We describe recent approaches to both dense and sparse parallel Cholesky factorization on shared memory multicore systems and present results from our new codes for problems arising from large real-world applications. In particular we describe our experiences using directed acyclic graph based scheduling in the dense case and retrofitting parallelism to a

sparse serial solver.

As research in topology optimisation has reached a level of maturity, two main classes of methods have emerged and their applications to real engineering design in industry are increasing. It has therefore become important to identify the limitations and challenges in order to ensure that topology optimisation is appropriately employed during the design process whilst research may continue to offer a more reliable and fast design tool to engineers.

The seminar will begin by introducing the topology optimisation problem and the two popular finite element based approaches. A range of numerical methods used in the typical implementations will be outlined. This will form the basis for the discussion on the short-comings and challenges as an easy-to-use design tool for engineers, particularly in the context of reliably providing the consistent optimum solutions to given problems with minimum a priori information. Another industrial requirement is a fast solution time to easy-to-set-up problems. The seminar will present the recent efforts in addressing some of these issues and the remaining challenges for the future.

Partial differential equations with a nonlinear pointwise constraint defined through a manifold occur in a variety of applications: The magnetization of a ferromagnet can be described by a unit length vector field and the orientation of the rod-like molecules that constitute a liquid crystal is often modeled by a vector field that attains its values in the real projective plane thus respecting the head-to-tail symmetry of the molecules. Other applications arise in geometric

modeling, quantum mechanics, and general relativity. Simple examples reveal that it is impossible to satisfy pointwise constraints exactly by lowest order finite elements. For two model problems we discuss the practical realization of the constraint, the efficient solution of the resulting nonlinear systems of equations, and weak accumulation of approximations at exact solutions.

We focus on domain decomposition methods for systems of PDEs (versus scalar PDEs). The Smith factorization (a "pure" algebra tool) is used systematically to derive new domain decompositions methods for symmetric and unsymmetric systems of PDEs: the compressible Euler equations, the Stokes and Oseen (linearized Navier-Stokes) problem. We will focus on the Stokes system. In two dimensions the key idea is the transformation of the Stokes problem into a scalar bi-harmonic problem. We show, how a proposed domain decomposition method for the bi-harmonic problem leads to a domain decomposition method for the Stokes equations which inherits the convergence behavior of the scalar problem. Thus, it is sufficient to study the convergence of the scalar algorithm. The same procedure can also be applied to the three-dimensional Stokes problem.

The butterfly-shaped Lorenz attractor is a fractal set made up of infinitely many periodic orbits. Ever since Lorenz (1963) introduced a system of three simple ordinary differential equations, much of the discussion of his system and its strange attractor has adopted a dynamical point of view. In contrast, we allow time to be a complex variable and look upon such solutions of the Lorenz system as analytic functions. Formal analysis gives the form and coefficients of the complex singularities of the Lorenz system. Very precise (> 500 digits) numerical computations show that the periodic orbits of the Lorenz system have singularities which obey that form exactly or very nearly so. Both formal analysis and numerical computation suggest that the mathematical analysis of the Lorenz system is a problem in analytic function theory. (Joint work with S. Sahutoglu).

The advent of multicore processors and other technologies like Graphical Processing Units (GPU) will considerably influence future research in High Performance Computing.

To take advantage of these architectures in dense linear algebra operations, new algorithms are

proposed that use finer granularity and minimize synchronization points.

After presenting some of these algorithms, we address the issue of pivoting and investigate randomization techniques to avoid pivoting in some cases.

In the particular case of GPUs, we show how linear algebra operations can be enhanced using

hybrid CPU-GPU calculations and mixed precision algorithms.

This is not intended to be a systematic History, but a selection of highlights, with some digressions, including:

The early days of the Computing Lab;

How the coming of the Computer changed some of the ways we do Computation;

A problem from the Study Groups;

Influence of the computing environment (hardware and software);

Convergence analysis for the heat equation, then and now.

Dirichlet to Neumann maps and their generalizations are exceptionally useful tools in the study of eigenvalue problems for ODEs and PDEs. They also have real physical significance through their occurrence in electrical impedance tomography, with applications to medical imagine, landmine detection and non-destructive testing. This talk will review some of the basic properties of Dirichlet to Neumann maps, some new abstract results which make it easier to use them for a wide variety of models, and some analytical/numerical results which depend on them, including detection and elimination of spectral pollution.

The Jacobi-Davidson method, proposed by Sleijpen and Van der Vorst more than a decade ago, has been successfully used to numerically solve large matrix eigenvalue problems. In this talk we will give an introduction to and an overview of this method, and also point out some recent developments.

The task is to compute orthogonal eigenvectors (without Gram-Schmidt) of symmetric tridiagonals for isolated clusters of close eigenvalues. We review an "old" method, the Submatrix method, and describe an extension which significantly enlarges the scope to include several mini-clusters within the given cluster. An essential feature is to find the envelope of the associated invariant subspace.

"Eigenvalue avoidance" or "level repulsion" refers to the tendency of eigenvalues of matrices or operators to be distinct rather than degenerate.

The mathematics goes back to von Neumann and Wigner in 1929 and touches many subjects including numerical linear algebra, random matrix theory, chaotic dynamics, and number theory.

This talk will be an informal illustrated discussion of various aspects of this phenomenon.

We discuss general and structured matrix polynomials which may be singular and may have eigenvalues at infinity. We discuss several real industrial applications ranging from acoustic field computations to optimal control problems.

We discuss linearization and first order formulations and their relationship to the corresponding techniques used in the treatment of systems of higher order differential equations.

In order to deal with structure preservation, we derive condensed/canonical forms that allow (partial) deflation of critical eigenvalues and the singular structure of the matrix polynomial. The remaining reduced order staircase form leads to new types of linearizations which determine the finite eigenvalues and corresponding eigenvectors.

Based on these new linearization techniques we discuss new structure preserving eigenvalue methods and present several real world numerical examples.

Meshfree methods become more and more important for the numerical simulation of complex real-world processes. Compared to classical, mesh-based methods they have the advantage of being more flexible, in particular for higher dimensional problems and for problems, where the underlying geometry is changing. However, often, they are also combined with classical methods to form hybrid methods.

In this talk, I will discuss meshfree, kernel based methods. After a short introduction along the lines of optimal recovery, I will concentrate on results concerning convergence orders and stability. After that I will address efficient numerical algorithms. Finally, I will present some examples, including one from fluid-structure-interaction, which will demonstrate why these methods are currently becoming Airbus's preferred solution in Aeroelasticity.

A well-known problem in protein modeling is the determination of the structure of a protein with a given set of inter-atomic or inter-residue distances obtained from either physical experiments or theoretical estimates. A general form of the problem is known as the distance geometry problem in mathematics, the graph embedding problem in computer science, and the multidimensional scaling problem in statistics. The problem has applications in many other scientific and engineering fields as well such as sensor network localization, image recognition, and protein classification. We describe the formulations and complexities of the problem in its various forms, and introduce a so-called geometric buildup approach to the problem. We present the general algorithm and discuss related computational issues including control of numerical errors, determination of rigid vs. unique structures, and tolerance of distance errors. The theoretical basis of the approach is established based on the theory of distance geometry. A group of necessary and sufficient conditions for the determination of a structure with a given set of distances using a geometric buildup algorithm are justified. The applications of the algorithm to model protein problems are demonstrated.

Graph-theoretic ideas have become very useful in understanding modern large-scale data-mining techniques. We show in this talk that ideas from optimization are also quite useful to better understand the numerical behavior of the corresponding algorithms. We illustrate this claim by looking at two specific graph theoretic problems and their application in data-mining.

The first problem is that of reputation systems where the reputation of objects and voters on the web are estimated; the second problem is that of estimating the similarity of nodes of large graphs. These two problems are also illustrated using concrete applications in data-mining.

Radial basis function (RBF) methods have been successfully used to approximate functions in multidimensional complex domains and are increasingly being used in the numerical solution of partial differential equations. These methods are often called meshfree numerical schemes since, in some cases, they are implemented without an underlying grid or mesh.

The focus of this talk is on the class of RBFs that allow exponential convergence for smooth problems. We will explore the dependence of accuracy and stability on node locations of RBF interpolants. Because Gaussian RBFs with equally spaced centers are related to polynomials through a change of variable, a number of precise conclusions about convergence rates based on the smoothness of the target function will be presented. Collocation methods for PDEs will also be considered.

Adaptive discretizations and iterative multilevel solvers are nowadays well established techniques for the numerical solution of PDEs.

The development of efficient multilevel techniques in the context of PDE-constrained optimization methods is an active research area that offers the potential of reducing the computational costs of the optimization process to an equivalent of only a few PDE solves.

We present a general class of inexact adaptive multilevel SQP-methods for PDE-constrained optimization problems. The algorithm starts with a coarse discretization of the underlying optimization problem and provides

1. implementable criteria for an adaptive refinement strategy of the current discretization based on local error estimators and

2. implementable accuracy requirements for iterative solvers of the PDE and adjoint PDE on the current grid

such that global convergence to the solution of the infinite-dimensional problem is ensured.

We illustrate how the adaptive refinement strategy of the multilevel SQP-method can be implemented by using existing reliable a posteriori error estimators for the state and the adjoint equation. Moreover, we discuss the efficient handling of control constraints and describe how efficent multilevel preconditioners can be constructed for the solution of the arising linear systems.

Numerical results are presented that illustrate the potential of the approach.

This is joint work with Jan Carsten Ziems.

The Snaer program calculates the posterior mean and variance of variables on some of which we have data (with precisions), on some we have prior information (with precisions), and on some prior indicator ratios (with precisions) are available. The variables must satisfy a number of exact restrictions. The system is both large and sparse. Two aspects of the statistical and computational development are a practical procedure for solving a linear integer system, and a stable linearization routine for ratios. We test our numerical method for solving large sparse linear least-squares estimation problems, and find that it performs well, even when the $n \times k$ design matrix is large ( $nk = O (10^{8})$ ).

The evaluation of oscillatory integrals is often considered to be a computationally challenging problem. However, in many cases, the exact opposite is true. Oscillatory integrals are cheaper to evaluate than non-oscillatory ones, even more so in higher dimensions. The simplest strategy that illustrates the general approach is to truncate an asymptotic expansion, where available. We show that an optimal strategy leads to the construction of certain unconventional Gaussian quadrature rules, that converge at twice the rate of asymptotic expansions. We examine a range of one-dimensional and higher dimensional, singular and highly oscillatory integrals.