Past Computational Mathematics and Applications Seminar

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.

The "secular equation" is a special way of expressing eigenvalue

problems in a variety of applications. We describe the secular

equation for several problems, viz eigenvector problems with a linear

constraint on the eigenvector and the solution of eigenvalue problems

where the given matrix has been modified by a rank one matrix. Next we

show how the secular equation can be approximated by use of the

Lanczos algorithm. Finally, we discuss numerical methods for solving

the approximate secular equation.

The numerical study of lattice quantum chromodynamics (QCD) is an attempt to extract predictions about the world around us from the standard model of physics. Worldwide, there are several large collaborations on lattice QCD methods, with terascale computing power dedicated to these problems. Central to the computation in lattice QCD is the inversion of a series of fermion matrices, representing the interaction of quarks on a four-dimensional space-time lattice. In practical computation, this inversion may be approximated based on the solution of a set of linear systems.

In this talk, I will present a basic description of the linear algebra problems in lattice QCD and why we believe that multigrid methods are well-suited to effectively solving them. While multigrid methods are known to be efficient solution techniques for many operators, those arising in lattice QCD offer new challenges, not easily handled by classical multigrid and algebraic multigrid approaches. The role of adaptive multigrid techniques in addressing the fermion matrices will be highlighted, along with preliminary results for several model problems.

In this talk we shall be looking at recent and forthcoming developments in the widely used LAPACK and ScaLAPACK numerical linear algebra libraries.

Improvements include the following: Faster algorithms, better numerical methods, memory hierarchy optimizations, parallelism, and automatic performance tuning to accommodate new architectures; more accurate algorithms, and the use of extra precision; expanded functionality, including updating and downdating and new eigenproblems; putting more of LAPACK into ScaLAPACK; and improved ease of use with friendlier interfaces in multiple languages. To accomplish these goals we are also relying on better software engineering techniques and contributions from collaborators at many institutions.

After an overview, this talk will highlight new more accurate algorithms; faster algorithms, including those for pivoted Cholesky and updating of factorizations; and hybrid data formats.

This is joint work with Jim Demmel, Jack Dongarra and the LAPACK/ScaLAPACK team.

Several formulations of structural optimization problems based on linear and nonlinear semidefinite programming will be presented. SDP allows us to formulate and solve problems with difficult constraints that could hardly be solved before. We will show that sometimes it is advantageous to prefer a nonlinear formulation to a linear one. All the presented formulations result in large-scale sparse (nonlinear) SDPs. In the second part of the talk we will show how these problems can be solved by our augmented Lagrangian code PENNON. Numerical examples will illustrate the talk.

Joint work with Michael Stingl.

We consider the parallel solution of sparse linear systems of equations in a limited memory environment. A preliminary out-of core version of a sparse multifrontal code called MUMPS (MUltifrontal Massively Parallel Solver) has been developed as part of a collaboration between CERFACS, ENSEEIHT and INRIA (ENS-Lyon and Bordeaux).

We first briefly describe the current status of the out-of-core factorization phase. We then assume that the factors have been written on the hard disk during the factorization phase and we discuss the design of an efficient solution phase.Two different approaches are presented to read data from the disk, with a discussion on the advantages and the drawbacks of each one.

Our work differs and extends the work of Rothberg and Schreiber (1999) and of Rotkin and Toledo (2004) because firstly we consider a parallel out-of-core context, and secondly we also study the performance of the solve phase.

This is work on collaboration with E. Agullo, I.S Duff, A. Guermouche, J.-Y. L'Excellent, T. Slavova

A strict adherence to threshold pivoting in the direct solution of symmetric indefinite problems can result in substantially more work and storage than forecast by an sparse analysis of the symmetric problem. One way of avoiding this is to use static pivoting where the data structures and pivoting sequence generated by the analysis are respected and pivots that would otherwise be very small are replaced by a user defined quantity. This can give a stable factorization but of a perturbed matrix.

The conventional way of solving the sparse linear system is then to use iterative refinement (IR) but there are cases where this fails to converge. We will discuss the use of more robust iterative methods, namely GMRES and its variant FGMRES and their backward stability when the preconditioning is performed by HSL_M57 with a static pivot option.

Several examples under Matlab will be presented.

Radial basis functions have been used for decades for the interpolation of scattered,

high-dimensional data. Recently they have attracted interest as methods for simulating

partial differential equations as well. RBFs do not require a grid or triangulation, they

offer the possibility of spectral accuracy with local refinement, and their implementation

is very straightforward. A number of theoretical and practical breakthroughs in recent years

has improved our understanding and application of these methods, and they are currently being

tested on real-world applications in shallow water flow on the sphere and tear film evolution

in the human eye.

Joint work with Yogi Erlangga and Kees Vuik.

Shifted Laplace preconditioners have attracted considerable attention as a technique to speed up convergence of iterative solution methods for the Helmholtz equation. In this paper we present a comprehensive spectral analysis of the Helmholtz operator preconditioned with a shifted Laplacian. Our analysis is valid under general conditions. The propagating medium can be heterogeneous, and the analysis also holds for different types of damping, including a radiation condition for the boundary of the computational domain. By combining the results of the spectral analysis of the preconditioned Helmholtz operator with an upper bound on the GMRES-residual norm we are able to provide an optimal value for the shift, and to explain the mesh-depency of the convergence of GMRES preconditioned with a shifted Laplacian. We illustrate our results with a seismic test problem.

Many large-scale optimization problems arise in the context of the discretization of infinite dimensional applications. In such cases, the description of the finite-dimensional problem is not unique, but depends on the discretization used, resulting in a natural multi-level description. How can such a problem structure be exploited, in discretized problems or more generally? The talk will focus on discussing this issue in the context of unconstrained optimization and in relation with the classical multigrid approach to elliptic systems of partial differential equations. Both theoretical convergence properties of special purpose algorithms and their numerical performances will be discussed. Perspectives will also be given.

Collaboration with S. Gratton, A. Sartenaer and M. Weber.