Past Computational Mathematics and Applications Seminar

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.

We present algorithms for dense LU and QR factorizations that minimize the cost of communication. One of today's challenging technology trends is the increased communication cost. This trend predicts that arithmetic will continue to improve exponentially faster than bandwidth, and bandwidth exponentially faster than latency. The new algorithms for dense QR and LU factorizations greatly reduce the amount of time spent communicating, relative to conventional algorithms.

This is joint work with James Demmel, Mark Hoemmen, Julien Langou, and Hua Xiang.

A new primal-dual augmented Lagrangian merit function is proposed that may be minimized with respect to both the primal and dual variables. A benefit of this approach is that each subproblem may be regularized by imposing explicit bounds on the dual variables. Two primal-dual variants of classical primal methods are given: a primal-dual bound constrained Lagrangian (pdBCL) method and a primal-dual l1 linearly constrained Lagrangian (pdl1-LCL) method.

Model reduction (also called system reduction, order reduction) is an ubiquitous tool in the analysis and simulation of dynamical systems, control design, circuit simulation, structural dynamics, CFD, etc. In the past decades many approaches have been developed for reducing the complexity of a given model. In this introductory talk, we will survey some of the most prominent methods used for linear systems, compare their properties and highlight similarities. In particular, we will emphasize the role of recent developments in numerical linear algebra in the different approaches. Efficiently using these techniques, the range of applicability of some of the methods has considerably widened.

The performance of several approaches will be demonstrated using real-world examples from a variety of engineering disciplines.

I will address the usage of the elliptic reconstruction technique (ERT) in a posteriori error analysis for fully discrete schemes for parabolic partial differential equations. A posteriori error estimates are effective tools in error control and adaptivity and a mathematical rigorous derivation justifies and improves their use in practical implementations.

The flexibility of the ERT allows a virtually indiscriminate use of various parabolic PDE techniques such as energy methods, duality methods and heat-kernel estimates, as opposed to direct approaches which leave less maneuver room. Thanks to ERT parabolic stability techniques can be combined with different elliptic a posteriori error analysis techniques, such as residual or recovery estimators, to derive a posteriori error bounds. The method has the merit of unifying previously known approaches, as well as providing new ones and providing us with novel error bounds (e.g., pointwise norm error bounds for the heat equation). [These results are based on joint work with Ch. Makridakis and A. Demlow.]

Another feature, which I would like to highlight, of the ERT is its simplifying power. It allows us to derive estimates where the analysis would be very complicated otherwise. As an example, I will illustrate its use in the context of non-conforming methods, with a special eye on discontinuous Galerkin methods. [These are recent results obtained jointly with E. Georgoulis.]

We consider the computational complexity of optimizing various classes

of continuous functions over a simplex, hypercube or sphere. These

relatively simple optimization problems arise naturally from diverse

applications. We review known approximation results as well as negative

(inapproximability) results from the recent literature.

Whether the 3D incompressible Euler or Navier-Stokes equations

can develop a finite time singularity from smooth initial data has been

an outstanding open problem. Here we review some existing computational

and theoretical work on possible finite blow-up of the 3D Euler equations.

We show that the geometric regularity of vortex filaments, even in an

extremely localized region, can lead to dynamic depletion of vortex

stretching, thus avoid finite time blowup of the 3D Euler equations.

Further, we perform large scale computations of the 3D Euler equations

to re-examine the two slightly perturbed anti-parallel vortex tubes which

is considered as one of the most attractive candidates for a finite time

blowup of the 3D Euler equations. We found that there is tremendous dynamic

depletion of vortex stretching and the maximum vorticity does not grow faster

than double exponential in time. Finally, we present a new class of solutions

for the 3D Euler and Navier-Stokes equations, which exhibit very interesting

dynamic growth property. By exploiting the special nonlinear structure of the

equations, we prove nonlinear stability and the global regularity of this class of solutions.

Many recent algorithmic approaches involve the construction of a differential equation model for computational purposes, typically by introducing an artificial time variable. The actual computational model involves a discretization of the now time-dependent differential system, usually employing forward Euler. The resulting dynamics of such an algorithm is then a discrete dynamics, and it is expected to be ''close enough'' to the dynamics of the continuous system (which is typically easier to analyze) provided that small -- hence many -- time steps, or iterations, are taken. Indeed, recent papers in inverse problems and image processing routinely report results requiring thousands of iterations to converge. This makes one wonder if and how the computational modeling process can be improved to better reflect the actual properties sought.

In this talk we elaborate on several problem instances that illustrate the above observations. Algorithms may often lend themselves to a dual interpretation, in terms of a simply discretized differential equation with artificial time and in terms of a simple optimization algorithm; such a dual interpretation can be advantageous. We show how a broader computational modeling approach may possibly lead to algorithms with improved efficiency.

The development and quantitative validation of complex nonlinear differential equation models is a difficult task that requires the support by numerical methods for sensitivity analysis, parameter estimation, and the optimal design of experiments. The talk first presents particularly efficient "simultaneous" boundary value problems methods for parameter estimation in nonlinear differential algebraic equations, which are based on constrained Gauss-Newton-type methods and a time domain decomposition by multiple shooting. They include a numerical analysis of the well-posedness of the problem and an assessment of the error of the resulting parameter estimates. Based on these approaches, efficient optimal control methods for the determination of one, or several complementary, optimal experiments are developed, which maximize the information gain subject to constraints such as experimental costs and feasibility, the range of model validity, or further technical constraints.

Special emphasis is placed on issues of robustness, i.e. how to reduce the sensitivity of the problem solutions with respect to uncertainties - such as outliers in the measurements for parameter estimation, and in particular the dependence of optimum experimental designs on the largely unknown values of the model parameters. New numerical methods will be presented, and applications will be discussed that arise in satellite orbit determination, chemical reaction kinetics, enzyme kinetics and robotics. They indicate a wide scope of applicability of the methods, and an enormous potential for reducing the experimental effort and improving the statistical quality of the models.

(Based on joint work with H. G. Bock, S. Koerkel, and J. P. Schloeder.)

Spectral methods are a class of methods for solving PDEs numerically.

If the solution is analytic, it is known that these methods converge

exponentially quickly as a function of the number of terms used.

The basic spectral method only works in regular geometry (rectangles/disks).

A huge amount of effort has gone into extending it to

domains with a complicated geometry. Domain decomposition/spectral

element methods partition the domain into subdomains on which the PDE

can be solved (after transforming each subdomain into a

regular one). We take the dual approach - embedding the domain into

a larger regular domain - known as the fictitious domain method or

domain embedding. This method is extremely simple to implement and

the time complexity is almost the same as that for solving the PDE

on the larger regular domain. We demonstrate exponential convergence

for Dirichlet, Neumann and nonlinear problems. Time permitting, we

shall discuss extension of this technique to PDEs with discontinuous

coefficients.

In this talk we will mainly analyze the vibrations of a simplified 1-d model for a multi-body structure consisting of a finite number of flexible strings distributed along a planar graph. In particular we shall analyze how solutions propagate along the graph as time evolves. The problem of the observation of waves is a natural framework to analyze this issue. Roughly, the question can be formulated as follows: Can we obtain complete information on the vibrations by making measurements in one single extreme of the network? This formulation is relevant both in the context of control and inverse problems.

Using the Fourier development of solutions and techniques of Nonharmonic Fourier Analysis, we give spectral conditions that guarantee the observability property to hold in any time larger than twice the total lengths of the network in a suitable Hilbert that can be characterized in terms of Fourier series by means of properly chosen weights. When the network graph is a tree these weights can be identified.

Once this is done these results can be transferred to other models as the Schroedinger, heat or beam-type equations.

This lecture is based on results obtained in collaboration with Rene Dager.