Past Computational Mathematics and Applications Seminar

We consider the problem of computing a nonnegative low rank factorization to a given nonnegative input matrix under the so-called "separabilty condition".  This assumption makes this otherwise NP hard problem polynomial time solvable, and we will use first order optimization techniques to compute such a factorization. The optimization model use is based on sparse regression with a self-dictionary, in which the low rank constraint is relaxed to the minimization of an l1-norm objective function.  We apply these techniques to endmember detection and classification in hyperspecral imaging data.

 The matrix logarithm, when applied to symmetric positive definite matrices, is known to satisfy a notable concavity property in the positive semidefinite (Loewner) order. This concavity property is a cornerstone result in the study of operator convex functions and has important applications in matrix concentration inequalities and quantum information theory.
In this talk I will show that certain rational approximations of the matrix logarithm remarkably preserve this concavity property and moreover, are amenable to semidefinite programming. Such approximations allow us to use off-the-shelf semidefinite programming solvers for convex optimization problems involving the matrix logarithm. These approximations are also useful in the scalar case and provide a much faster alternative to existing methods based on successive approximation for problems involving the exponential/relative entropy cone. I will conclude by showing some applications to problems arising in quantum information theory.

This is joint work with James Saunderson (Monash University) and Pablo Parrilo (MIT)

Rational Krylov methods are applicable to a wide range of scientific computing problems, and ​the rational Arnoldi algorithm is a commonly used procedure for computing an ​orthonormal basis of a rational Krylov space. Typically, the computationally most expensive component of this​ ​algorithm is the solution of a large linear system of equations at each iteration. We explore the​ ​option of solving several linear systems simultaneously, thus constructing the rational Krylov​ ​basis in parallel. If this is not done carefully, the basis being orthogonalized may become badly​ ​conditioned, leading to numerical instabilities in the orthogonalization process. We introduce the​ ​new concept of continuation pairs which gives rise to a near-optimal parallelization strategy that ​allows to control the growth of the condition number of this nonorthogonal basis. As a consequence we obtain a significantly more accurate and reliable parallel rational Arnoldi algorithm.
​ ​
The computational benefits are illustrated using several numerical examples from different application areas.
​ ​
This ​talk is based on joint work with Mario Berljafa  available as an Eprint at http://eprints.ma.man.ac.uk/2503/
 

We present global rates of convergence for a general class of methods for nonconvex smooth optimization that include linesearch, trust-region and regularisation strategies, but that allow inaccurate problem information. Namely, we assume the local (first- or second-order) models of our function are only sufficiently accurate with a certain probability, and they can be arbitrarily poor otherwise. This framework subsumes certain stochastic gradient analyses and derivative-free techniques based on random sampling of function values. It can also be viewed as a robustness
assessment of deterministic methods and their resilience to inaccurate derivative computation such as due to processor failure in a distribute framework. We show that in terms of the order of the accuracy, the evaluation complexity of such methods is the same as their counterparts that use deterministic accurate models; the use of probabilistic models only increases the complexity by a constant, which depends on the probability of the models being good. Time permitting, we also discuss the case of inaccurate, probabilistic function value information, that arises in stochastic optimization. This work is joint with Katya Scheinberg (Lehigh University, USA).
 

For linear dynamical systems, model reduction has achieved great success. In the case of linear dynamics,  we know how to construct, at a modest cost, (locally) optimal, input-independent reduced models; that is, reduced models that are uniformly good over all inputs having bounded energy. In addition, in some cases we can achieve this goal using only input/output data without a priori knowledge of internal  dynamics.  Even though model reduction has been successfully and effectively applied to nonlinear dynamical systems as well, in this setting,  bot the reduction process and the reduced models are input dependent and the high fidelity of the resulting approximation is generically restricted to the training input/data. In this talk, we will offer remedies to this situation.

To predict the behaviour of a dynamical system using a mathematical model, an accurate estimate of the current state of the system is needed in order to initialize the model. Complete information on the current state is, however, seldom available. The aim of optimal state estimation, known in the geophysical sciences as ‘data assimilation’, is to determine a best estimate of the current state using measured observations of the real system over time, together with the model equations. The problem is commonly formulated in variational terms as a very large nonlinear least-squares optimization problem. The lack of complete data, coupled with errors in the observations and in the model, leads to a highly ill-conditioned inverse problem that is difficult to solve.

To understand the nature of the inverse problem, we examine how different components of the assimilation system influence the conditioning of the optimization problem. First we consider the case where the dynamical equations are assumed to model the real system exactly. We show, against intuition, that with increasingly dense and precise observations, the problem becomes harder to solve accurately. We then extend these results to a 'weak-constraint' form of the problem, where the model equations are assumed not to be exact, but to contain random errors. Two different, but mathematically equivalent, forms of the problem are derived. We investigate the conditioning of these two forms and find, surprisingly, that these have quite different behaviour.

The requirement to compute Jordan blocks for multiple eigenvalues arises in a number of physical problems, for example panel flutter problems in aerodynamical stability, the stability of electrical power systems, and in quantum mechanics. We introduce a general method for computing a 2-dimensional Jordan block in a parameter-dependent matrix eigenvalue problem based on the so called Implicit Determinant Method. This is joint work with Alastair Spence (Bath).

In many applications we need to find or estimate the $p \ge 1$ largest elements of a matrix, along with their locations. This is required for recommender systems used by Amazon and Netflix, link prediction in graphs, and in finding the most important links in a complex network, for example. 

Our algorithm uses only matrix vector products and is based upon a power method for mixed subordinate norms. We have obtained theoretical results on the convergence of this algorithm via a comparison with rook pivoting for the LU  decomposition. We have also improved the practicality of the algorithm by producing a blocked version iterating on $n \times t$ matrices, as opposed to vectors, where $t$ is a tunable parameter. For $p > 1$ we show how deflation can be used to improve the convergence of the algorithm. 

Finally, numerical experiments on both randomly generated matrices and real-life datasets (the latter for $A^TA$ and $e^A$) show how our algorithms can reliably estimate the largest elements of a matrix whilst obtaining considerable speedups when compared to forming the matrix explicitly: over 1000x in some cases.

For many years, sum-factored algorithms for finite elements in rectangular reference geometry have combined low complexity with the mathematical power of high-order approximation.  However, such algorithms rely heavily on the tensor product structure inherent in the geometry and basis functions, and similar algorithms for simplicial geometry have proven elusive.

Bernstein polynomials are totally nonnegative, rotationally symmetric, and geometrically decomposed bases with many other remarkable properties that lead to optimal-complexity algorithms for element wise finite element computations.  The also form natural building blocks for the finite element exterior calculus bases for the de Rham complex so that H(div) and H(curl) bases have efficient representations as well.  We will also their relevance for explicit discontinuous Galerkin methods, where the element mass matrix requires special attention.

We will review the basic building blocks of iterative solvers, i.e. sparse matrix-vector multiplication, in the context of GPU devices such 
as the cards by NVIDIA; we will then discuss some techniques in preconditioning by approximate inverses, and we will conclude with an 
application to an image processing problem from the biomedical field.

The aim of this talk is to design an efficient multigrid method for constrained convex optimization problems arising from discretization  of  some  underlying  infinite  dimensional  problems. Due  to problem  dependency  of this approach, we only consider bound constraints with (possibly) a linear equality constraint. As our aim is to target large-scale problems, we want to avoid computation of second 
derivatives of the objective function, thus excluding Newton like methods. We propose a smoothing operator that only uses first-order information and study the computational efficiency of the resulting method. In the second part, we consider application of multigrid techniques to more general optimization problems, in particular, the topology design problem.

Quadrature is the term for the numerical evaluation of integrals.  It's a beautiful subject because it's so accessible, yet full of conceptual surprises and challenges.  This talk will review ten of these, with plenty of history and numerical demonstrations.  Some are old if not well known, some are new, and two are subjects of my current research.

Partial differential equations with more than three coordinates arise naturally if the model features certain kinds of stochasticity. Typical examples are the Schroedinger, Fokker-Planck and Master equations in quantum mechanics or cell biology, as well as quantification of uncertainty.
The principal difficulty of a straightforward numerical solution of such equations is the `curse of dimensionality': the storage cost of the discrete solution grows exponentially with the number of coordinates (dimensions).

One way to reduce the complexity is the low-rank separation of variables. One can see all discrete data (such as the solution) as multi-index arrays, or tensors. These large tensors are never stored directly.
We approximate them by a sum of products of smaller factors, each carrying only one of the original variables. I will present one of the simplest but powerful of such representations, the Tensor Train (TT) decomposition. The TT decomposition generalizes the approximation of a given matrix by a low-rank matrix to the tensor case. It was found that many interesting models allow such approximations with a significant reduction of storage demands.

A workhorse approach to computations with the TT and other tensor product decompositions is the alternating optimization of factors. The simple realization is however prone to convergence issues.
I will show some of the recent improvements that are indispensable for really many dimensions, or solution of linear systems with non-symmetric or indefinite matrices.

To face the advent of multicore processors and the ever increasing complexity of hardware architectures, programming
models based on DAG parallelism regained popularity in the high performance, scientific computing community. Modern runtime systems offer a programming interface that complies with this paradigm and powerful engines for scheduling the tasks into which the application is decomposed. These tools have already proved their effectiveness on a number of dense linear algebra applications. 

In this talk we present the design of task-based sparse direct solvers on top of runtime systems. In the context of the
qr_mumps solver, we prove the usability and effectiveness of our approach with the implementation of a sparse matrix multifrontal factorization based on a Sequential Task flow parallel programming model. Using this programming model, we developed features such as the integration of dense 2D Communication Avoiding algorithms in the multifrontal method allowing for better scalability compared to the original approach used in qr_mumps.

Following this approach, we move to heterogeneous architectures where task granularity and scheduling strategies are critical to achieve performance. We present, for the multifrontal method, a hierarchical strategy for data partitioning and a scheduling algorithm capable of handling the heterogeneity of resources.   Finally we introduce a memory-aware algorithm to control the memory behavior of our solver and show, in the context of multicore architectures, an important reduction of the memory footprint for the multifrontal QR factorization with a small impact on performance.

Functions are usually approximated numerically in a basis, a non-redundant and complete set of functions that span a certain space. In this talk we highlight a number of benefits of using overcomplete sets, in particular using the more general notion of a "frame". The main 

benefit is that frames are easily constructed even for functions of several variables on domains with irregular shapes. On the other hand, allowing for possible linear depencies naturally leads to ill-conditioning of approximation algorithms. The ill-conditioning is 

potentially severe. We give some useful examples of frames and we first address the numerical stability of best approximations in a frame. Next, we briefly describe special point sets in which interpolation turns out to be stable. Finally, we review so-called Fourier extensions and an efficient algorithm to approximate functions with spectral accuracy on domains without structure.

When assigned with the task of extracting information from given image data the first challenge one faces is the derivation of a truthful model for both the information and the data. Such a model can be determined by the a-priori knowledge about the image (information), the data and their relation to each other. The source of this knowledge is either our understanding of the type of images we want to reconstruct and of the physics behind the acquisition of the data or we can thrive to learn parametric models from the data itself. The common question arises: how can we customise our model choice to a particular application? Or better how can we make our model adaptive to the given data?

Starting from the first modelling strategy this talk will lead us from nonlinear diffusion equations and subdifferential inclusions of total variation type functionals as the most successful image modeltoday to non-smooth second- and third-order variational models, with data models for Gaussian and Poisson distributed data as well as impulse noise. These models exhibit solution-dependent adaptivities in form of nonlinearities or non-smooth terms in the PDE or the variational problem, respectively. Applications for image denoising, inpainting and surface reconstruction are given. After a critical discussion of these different image and data models we will turn towards the second modelling strategy and propose to combine it with the first one using a PDE constrained optimisation method that customises a parametrised form of the model by learning from examples. In particular, we will consider optimal parameter derivation for total variation denoising with multiple noise distributions and optimising total generalised variation regularisation for its application in photography.

Equations of quantum mechanics in the semiclassical regime present an enduring challenge for numerical analysts, because their solution is highly oscillatory and evolves on two scales. Standard computational approaches to the semiclassical Schrödinger equation do not allow for long time integration as required, for example, in quantum control of atoms by short laser bursts. This has motivated our approach of asymptotic splittings. Combining techniques from Lie-algebra theory and numerical algebra, we present a new computational paradigm of symmetric Zassenhaus splittings, which lends itself to a very precise discretisation in long time intervals, at very little cost. We will illustrate our talk by examples of quantum phenomena – quantum tunnelling and quantum scattering – and their computation and, time allowing, discuss an extension of this methodology to time-dependent semiclassical systems using Magnus expansions

This talk focuses on the direct search method, arguably one of the simplest optimization algorithms. The algorithm minimizes an objective function by iteratively evaluating it along a number of (polling) directions, which are typically taken from so-called positive spanning sets. It does not use derivatives.

We first introduce the worst case complexity theory of direct search, and discuss how to choose the positive spanning set to minimize the complexity bound. The discussion leads us to a long-standing open
problem in Discrete Geometry. A recent result on this problem enables us to establish the optimal order for the worst case complexity of direct search.

We then show how to achieve even lower complexity bound by using random polling directions. It turns out that polling along two random directions at each iteration is sufficient to guarantee the convergence
of direct search for any dimension, and the resultant algorithm enjoys lower complexity both in theory and in practice.

The last part of the talk is devoted to direct search based on inaccurate function values. We address three questions:
i) what kind of solution can we obtain by direct search if the function values are inaccurate? 
ii) what is the worst case complexity to attain such a solution? iii) given
the inaccuracy in the function values, when to stop the algorithm in order
to guarantee the quality of the solution and also avoid “over-optimization”?

This talk is based on joint works with F. Delbos, M. Dodangeh, S. Gratton, B. Pauwels, C. W. Royer, and L. N. Vicente.

I will review some recent work on the problem of reliable automatic detection of blow-up behaviour for nonlinear parabolic PDEs. The adaptive algorithms developed are based on rigorous conditional a posteriori error bounds. The use of space-time adaptivity is crucial in making the problem computationally tractable. The results presented are applicable to quite general spatial operators, rendering the approach potentially useful in informing respective PDE theory. The new adaptive algorithm is shown to accurately estimate the blow-up time of a number of problems, including ones exhibiting regional blow-up. 

The talk will introduce the concepts of multilevel optimization and motivate them in the context of problems arising from the discretization of infinite dimensional applications. It will be shown how optimization methods can accomodate a number of useful (and classical) ideas from the multigrid
community, and thereby produce substantial efficiency improvements compared to existing large-scale minimization techniques.  Two different classes of multilevel methods will be discussed: trust-region and linesearch algorithms.
The first class will be described in the context of a multilevel generalization of the well-known trust-region-Newton method.  The second will focus on limited-memory quasi-Newton algorithms.  Preliminary numerical results will be presented which indicate that both types of multilevel algorithms may be practically very advantageous.

In numerical atmosphere models, values of relevant physical parameters are often uncertain by more than 100% and weather forecast skill is significantly reduced after a couple of days. Still, numerical operations are typically calculated in double precision with 15 significant decimal digits. If we reduce numerical precision, we can reduce power consumption and increase computational performance significantly. If savings are reinvested to build larger supercomputers, this would allow an increase in resolution in weather and climate models and might lead to better predictions of future weather and climate. 
I will discuss approaches to reduce numerical precision beyond single precision in high performance computing and in particular in weather and climate modelling. I will present results that show that precision can be reduced significantly in atmosphere models and that potential savings can be huge. I will also discuss how rounding errors will impact model dynamics and interact with model uncertainty and predictability.

We propose the use of a constraint preconditioner for the iterative solution of the linear system arising from the finite element discretization of the coupled Stokes-Darcy system. The Stokes-Darcy system is a set of coupled PDEs that can be used to model a freely flowing fluid over porous media flow. The fully coupled system matrix is large, sparse, non-symmetric, and of saddle point form. We provide for exact versions of the constraint preconditioner spectral and field-of-values bounds that are independent of the underlying mesh width. We present several numerical experiments, using the deal.II finite element library, that illustrate our results in both two and three dimensions. We compare exact and inexact versions of the constraint preconditioner against standard block diagonal and block lower triangular preconditioners to illustrate its favorable properties.

When formulated appropriately, the broad families of sequential quadratic programming, augmented Lagrangian and interior-point methods all require the solution of symmetric saddle-point linear systems. When regularization is employed, the systems become symmetric and quasi definite. The latter are
indefinite but their rich structure and strong relationships with definite systems enable specialized linear algebra, and make them prime candidates for matrix-free implementations of optimization methods. In this talk, I explore various formulations of the step equations in optimization and corresponding
iterative methods that exploit their structure.

Security Constrained Optimal Power Flow is an increasingly important problem for power systems operation both in its own right and as a subproblem for more complex problems such as transmission switching or
unit commitment.

The structure of the problem resembles stochastic programming problems in that one aims to find a cost optimal operation schedule that is feasible for all possible equipment outage scenarios
(contingencies). Due to the presence of power flow constraints (in their "DC" or "AC" version), the resulting problem is a large scale linear or nonlinear programming problem.

However it is known that only a small subset of the contingencies is active at the solution. We show how Interior Point methods can exploit this structure both by simplifying the linear algebra operations as
well as generating necessary contingencies on the fly and integrating them into the algorithm using IPM warmstarting techniques. The final problem solved by this scheme is significantly smaller than the full
contingency constrained problem, resulting in substantial speed gains.

Numerical and theoretical results of our algorithm will be presented.

Can we extend the FEM to general polytopic, i.e. polygonal and polyhedral, meshes while retaining 
the ease of implementation and computational cost comparable to that of standard FEM? Within this talk, I present two approaches that achieve just  that (and much more): the Virtual Element Method (VEM) and an hp-version discontinuous Galerkin (dG) method.

The Virtual Element spaces are like the usual (polynomial) finite element spaces with the addition of suitable non-polynomial functions. This is far from being a novel idea. The novelty of the VEM approach is that it avoids expensive evaluations of the non-polynomial "virtual" functions by basing all 
computations solely on the method's carefully chosen degrees of freedom. This way we can easily deal 
with complicated element geometries and/or higher continuity requirements (like C1, C2, etc.), while 
maintaining the computational complexity comparable to that of standard finite element computations.

As you might expect, the choice and number of the degrees of freedom depends on such continuity 
requirements. If mesh flexibility is the goal, while one is ready to  give up on global regularity, other approaches can be considered. For instance, dG approaches are naturally suited to deal with polytopic meshes. Here I present an extension of the classical Interior Penalty dG method which achieves optimal rates of convergence on polytopic meshes even under elemental edge/face degeneration. 

The next step is to exploit mesh flexibility in the efficient resolution of problems characterised by 
complicated geometries and solution features, for instance within the framework of automatic FEM 
adaptivity. I shall finally introduce ongoing work in this direction.

In this seminar I will present a semi-langrangian discretisation of the Monge-Ampère operator, which is of interest in optimal transport 
and differential geometry as well as in related fields of application.

I will discuss the proof of convergence to viscosity solutions. To address the challenge of uniqueness and convexity we draw upon the classical relationship with Hamilton-Jacobi-Bellman equations, which we extend to the viscosity setting. I will explain that the monotonicity of semi-langrangian schemes implies that they possess large stencils, which in turn requires careful treatment of the boundary conditions.

The contents of the seminar is based on current work with X Feng from the University of Tennessee.

The Singular Value Decomposition (SVD) of matrices and the related Principal Components Analysis (PCA) express a matrix in terms of singular vectors, which are linear combinations of all the input data and lack an intuitive physical interpretation. Motivated by the application of PCA and SVD in the analysis of populations genetics data, we will discuss the notion of leverage scores: a simple statistic that reveals columns/rows of a matrix that lie in the subspace spanned by the top principal components (left/right singular vectors). We will then use the leverage scores to present matrix decompositions that express the structure in a matrix in terms of actual columns (and/or rows) of the matrix. Such decompositions are easier to interpret in applications, since the selected columns and rows are subsets of the data. We will also discuss extensions of the leverage scores to reveal influential entries of a matrix.

We present a trust region algorithm for solving nonconvex optimization problems that, in the worst-case, is able to drive the norm of the gradient of the objective below a prescribed threshold $\epsilon > 0$ after at most ${\cal O}(\epsilon^{-3/2})$ function evaluations, gradient evaluations, or iterations.  Our work has been inspired by the recently proposed Adaptive Regularisation framework using Cubics (i.e., the ARC algorithm), which attains the same worst-case complexity bound.  Our algorithm is modeled after a traditional trust region algorithm, but employs modified step acceptance criteria and a novel trust region updating mechanism that allows it to achieve this desirable property.  Importantly, our method also maintains standard global and fast local convergence guarantees.

Although Toeplitz matrices are often dense, matrix-vector products with Toeplitz matrices can be quickly performed via circulant embedding and the fast Fourier transform. This makes their solution by preconditioned Krylov subspace methods attractive. 

For a wide class of symmetric Toeplitz matrices, symmetric positive definite circulant preconditioners that cluster eigenvalues have been proposed. MINRES or the conjugate gradient method can be applied to these problems and descriptive convergence theory based on eigenvalues guarantees fast convergence. 

In contrast, although circulant preconditioners have been proposed for nonsymmetric Toeplitz systems, guarantees of fast convergence are generally only available for CG for the normal equations (CGNE). This is somewhat unsatisfactory because CGNE has certain drawbacks, including slow convergence and a larger condition number. In this talk we discuss a simple alternative symmetrization of nonsymmetric Toeplitz matrices, that allows us to use MINRES to solve the resulting linear system. We show how existing circulant preconditioners for nonsymmetric Toeplitz matrices can be straightforwardly adapted to this situation and give convergence estimates similar to those in the symmetric case.

Numerical simulation tools for fluid and solid mechanics are often based on the discretisation of coupled systems of partial differential equations, which can easily be identified in terms of physical
conservation laws.  In contrast, much physical insight is often gained from the equivalent formulation of the relevant energy or free-energy functional, possibly subject to constraints.  Motivated by the
nonlinear static and dynamic behaviour of nematic liquid crystals and of magnetosensitive elastomers, we propose a finite-element framework for minimising these free-energy functionals, using Lagrange multipliers to enforce additional constraints.  This talk will highlight challenges, limitations, and successes, both in the formulation of these models and their use in numerical simulation.
This is joint work with PhD students Thomas Benson, David Emerson, and Dong Han, and with James Adler, Timothy Atherton, and Luis Dorfmann.

Preconditioning is of significant importance in the solution of large dimensional systems of linear equations such as those that arise from the numerical solution of partial differential equation problems. In this talk we will attempt a broad ranging review of preconditioning.

In this talk, we discuss the development of fast iterative solvers for matrix systems arising from various constrained optimization problems. In particular, we seek to exploit the saddle point structure of these problems to construct powerful preconditioners for the resulting systems, using appropriate approximations of the (1,1)-block and Schur complement.

The problems we consider arise from two well-studied subject areas within computational optimization. Specifically, we investigate the
numerical solution of PDE-constrained optimization problems, and the interior point method (IPM) solution of linear/quadratic programming
problems. Indeed a particular focus in this talk is the interior point method solution of PDE-constrained optimization problems with
additional inequality constraints on the state and control variables.

We present a range of optimization problems which we seek to solve using our methodology, and examine the theoretical and practical
convergence properties of our iterative methods for these problems.
 

Stability of the hp-Raviart-Thomas projection operator as a mapping H^1(K) -> H^1(K) on the unit cube K in R^3 has been addressed e.g. in [2], see also [1]. These results are suboptimal with respect to the polynomial degree. In this talk we present quasi-optimal stability estimates for the hp-Raviart-Thomas projection operator on the cube. The analysis involves elements of the polynomial approximation theory on an interval and the real method of Banach space interpolation.

(Joint work with Herbert Egger, TU Darmstadt)

[1] Mark Ainsworth and Katia Pinchedez. hp-approximation theory for BDFM and RT finite elements on quadrilaterals. SIAM J. Numer. Anal., 40(6):2047–2068 (electronic) (2003), 2002.

[2] Dominik Schötzau, Christoph Schwab, and Andrea Toselli. Mixed hp-DGFEM for incompressible flows. SIAM J. Numer. Anal., 40(6):2171–2194 (electronic) (2003), 2002.

Nonlinear systems of partial differential equations (PDEs) may permit several distinct solutions. The typical current approach to finding distinct solutions is to start Newton's method with many different initial guesses, hoping to find starting points that lie in different basins of attraction. In this talk, we present an infinite-dimensional deflation algorithm for systematically modifying the residual of a nonlinear PDE problem to eliminate known solutions from consideration. This enables the Newton--Kantorovitch iteration to converge to several different solutions, even starting from the same initial guess. The deflated Jacobian is dense, but an efficient preconditioning strategy is devised, and the number of Krylov iterations is observed not to grow as solutions are deflated. The technique is then applied to computing distinct solutions of nonlinear PDEs, tracing bifurcation diagrams, and to computing multiple local minima of PDE-constrained optimisation problems.

In this talk we will describe the steps towards self-consistent computer simulations of the evolution of the universe beginning soon after the Big Bang and ending with the formation of realistic stellar systems like the Milky Way. This is a multi-scale problem of vast proportions. The first step has been the Millennium Simulation, one of the largest and most successful numerical simulations of the Universe ever carried out. Now we are in the midst of the next step, where this is carried to a much higher level of physical fidelity on the latest supercomputing platforms. This talk will be illustrate how the role of mathematics is essential in this endeavor. Also computer simulations will be shown. This is joint work among others with Volker Springel.

Some problems in scientific computing, like the forward simulation of electromagnetic waves in geophysical prospecting, can be
solved via approximation of f(A)b, the action of a large matrix function f(A) onto a vector b. Iterative methods based on rational Krylov
spaces are powerful tools for these computations, and the choice of parameters in these methods is an active area of research.
We provide an overview of different approaches for obtaining optimal parameters, with an emphasis on the exponential and resolvent function, and the square root. We will discuss applications of the rational Arnoldi method for iteratively generating near-optimal absorbing boundary layers for indefinite Helmholtz problems, and for rational least squares vector fitting.

Gaussian quadrature rules are of theoretical and practical interest because of their role in numerical integration and interpolation. For general weighting functions, their computation can be performed with the Golub-Welsch algorithm or one of its refinements. However, for the specific case of Gauss-Legendre quadrature, computation methods based on asymptotic series representations of the Legendre polynomials have recently been proposed. 
For large quadrature rules, these methods provide superior accuracy and speed at the cost of generality. We provide an overview of the progress that was made with these asymptotic methods, focusing on the ideas and asymptotic formulas that led to them. 
Finally, the limited generality will be discussed with Gauss-Jacobi quadrature rules as a prominent possibility for extension.

The usual variational formulations of the Helmholtz equation are sign-indefinite (i.e. not coercive). In this talk, I will argue that this indefiniteness is not an inherent feature of the Helmholtz equation itself, only of its standard formulations. I will do this by presenting new sign-definite formulations of several Helmholtz boundary value problems.

This is joint work with Andrea Moiola (Reading).
 

Incomplete Cholesky factorizations are commonly used as black-box preconditioners for the iterative solution of large sparse symmetric positive definite linear systems. Traditionally, incomplete 
factorizations are obtained by dropping (i.e., replacing by zero) some entries of the factors during the factorization process. Here we consider a less common way to approximate the factors : through low-rank approximations of some off-diagonal blocks. We focus more specifically on approximation schemes that satisfy the orthogonality condition: the approximation should be orthogonal to the corresponding approximation error.

The resulting incomplete Cholesky factorizations have attractive theoretical properties. First, the underlying factorization process can be shown breakdown-free. Further, the condition number of the 
preconditioned system, that characterizes the convergence rate of standard iterative schemes, can be shown bounded as a function of the accuracy of individual approximations. Hence, such a bound can benefit from better approximations, but also from some algorithmic peculiarities. Eventually, the above results can be shown to hold for any symmetric positive definite system matrix.

On the practical side, we consider a particular variant of the preconditioner. It relies on a nested dissection ordering of unknowns to  insure an attractive memory usage and operations count. Further, it exploits in an algebraic way the low-rank structure present in system matrices that arise from PDE discretizations. A preliminary implementation of the method is compared with similar Cholesky and 
incomplete Cholesky factorizations based on dropping of individual entries.