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.

# Past Computational Mathematics and Applications Seminar

The Dynamic Dictionary of Mathematical Functions (or DDMF, http://ddmf.msr-inria.inria.fr/) is an interactive website on special functions inspired by reference books such as the NIST Handbook of Special Functions. The originality of the DDMF is that each of its “chapters” is automatically generated from a short mathematical description of the corresponding function.

To make this possible, the DDMF focuses on so-called D-finite (or holonomic) functions, i.e., complex analytic solutions of linear ODEs with polynomial coefficients. D-finite functions include in particular most standard elementary functions (exp, log, sin, sinh, arctan...) as well as many of the classical special functions of mathematical physics (Airy functions, Bessel functions, hypergeometric functions...). A function of this class can be represented by a finite amount of data (a differential equation along with sufficiently many initial values),

and this representation makes it possible to develop a computer algebra framework that deals with the whole class in a unified way, instead of ad hoc algorithms and code for each particular function. The DDMF attempts to put this idea into practice.

In this talk, I will present the DDMF, some of the algorithms and software libraries behind it, and ongoing projects based on similar ideas, with an emphasis on symbolic-numeric algorithms.

The coefficients in mathematical models of physical processes are often impossible to determine fully or accurately, and are hence subject to uncertainty. It is of great importance to quantify the uncertainty in the model outputs based on the (uncertain) information that is available on the model inputs. This invariably leads to very high dimensional quadrature problems associated with the computation of statistics of quantities of interest, such as the time it takes a pollutant plume in an uncertain subsurface flow problem to reach the boundary of a safety region or the buckling load of an airplane wing. Higher order methods, such as stochastic Galerkin or polynomial chaos methods, suffer from the curse of dimensionality and when the physical models themselves are complex and computationally costly, they become prohibitively expensive in higher dimensions. Instead, some of the most promising approaches to quantify uncertainties in continuum models are based on Monte Carlo sampling and the “multigrid philosophy”. Multilevel Monte Carlo (MLMC) Methods have been introduced recently and successfully applied to many model problems, producing significant gains. In this talk I want to recall the classical MLMC method and then show how the gains can be improved further (significantly) by using quasi-Monte Carlo (QMC) sampling rules. More importantly the dimension independence and the improved gains can be justified rigorously for an important model problem in subsurface flow. To achieve uniform bounds, independent of the dimension, it is necessary to work in infinite dimensions and to study quadrature in sequence spaces. I will present the elements of this new theory for the case of lognormal random coefficients in a diffusion problem and support the theory with numerical experiments.

For many tomographic imaging problems there are explicit inversion formulas, and depending on the completeness of the data these are unstable to differing degrees. Increasingly we are solving tomographic problems as though they were any other linear inverse problem using numerical linear algebra. I will illustrate the use of numerical singular value decomposition to explore the (in)stability for various problems. I will also show how standard techniques from numerical linear algebra, such as conjugate gradient least squares, can be employed with systematic regularization compared with the ad hoc use of slowly convergent iterative methods more traditionally used in computed tomography. I will mainly illustrate the talk with examples from three dimensional x-ray tomography but I will also touch on tensor tomography problems.

We outline a path from polynomial numerical hulls to multicentric calculus for evaluating f(A). Consider

$$Vp(A) = {z ∈ C : |p(z)| ≤ kp(A)k}$$

where p is a polynomial and A a bounded linear operator (or matrix). Intersecting these sets over polynomials of degree 1 gives the closure of the numerical range, while intersecting over all polynomials gives the spectrum of A, with possible holes filled in.

Outside any set Vp(A) one can write the resolvent down explicitly and this leads to multicentric holomorphic functional calculus.

The spectrum, pseudospectrum or the polynomial numerical hulls can move rapidly in low rank perturbations. However, this happens in a very controlled way and when measured correctly one gets an identity which shows e.g. the following: if you have a low-rank homotopy between self-adjoint and quasinilpotent, then the identity forces the nonnormality to increase in exact compensation with the spectrum shrinking.

In this talk we shall mention how the multicentric calculus leads to a nontrivial extension of von Neumann theorem

$$kf(A)k ≤ sup |z|≤1

kf(z)k$$

where A is a contraction in a Hilbert space, and conclude with some new results on (nonholomorphic) functional calculus for operators for which p(A) is normal at a nontrivial polynomial p. Notice that this is always true for matrices.

In numerical analysis the design and analysis of computational methods is often based on, and closely linked to, a well-posedness result for the underlying continuous problem. In particular the continuous dependence of the continuous model is inherited by the computational method when such an approach is used. In this talk our aim is to design a stabilised finite element method that can exploit continuous dependence of the underlying physical problem without making use of a standard well-posedness result such as Lax-Milgram's Lemma or The Babuska-Brezzi theorem. This is of particular interest for inverse problems or data assimilation problems which may not enter the framework of the above mentioned well-posedness results, but can nevertheless satisfy some weak continuous dependence properties. First we will discuss non-coercive elliptic and hyperbolic equations where the discrete problem can be ill-posed even for well posed continuous problems and then we will discuss the linear elliptic Cauchy problem as an example of an ill-posed problem where there are continuous dependence results available that are suitable for the framework that we propose.

Gradient-based optimisation techniques have become a common tool in the analysis of fluid systems. They have been applied to replace and extend large-scale matrix decompositions to compute optimal amplification and optimal frequency responses in unstable and stable flows. We will show how to efficiently extract linearised and adjoint information directly from nonlinear simulation codes and how to use this information for determining common flow characteristics. We also extend this framework to deal with the optimisation of less common norms. Examples from aero-acoustics and mixing will be presented.

Mathematical imaging is not only a multidisciplinary research area but also a major cross-discipline subject within mathematical sciences as image analysis techniques involve differential geometry, optimization, nonlinear partial differential equations (PDEs), mathematical analysis, computational algorithms and numerical analysis. Segmentation refers to the essential problem in imaging and vision of automatically detecting objects in an image.

In this talk I first review some various models and techniques in the variational framework that are used for segmentation of images, with the purpose of discussing the state of arts rather than giving a comprehensive survey. Then I introduce the practically demanding task of detecting local features in a large image and our recent segmentation methods using energy minimization and PDEs. To ensure uniqueness and fast solution, we reformulate one non-convex functional as a convex one and further consider how to adapt the additive operator splitting method for subsequent solution. Finally I show our preliminary work to attempt segmentation of blurred images in the framework of joint deblurring and segmentation.

This talk covers joint work with Jianping Zhang, Lavdie Rada, Bryan Williams, Jack Spencer (Liverpool, UK), N. Badshah and H. Ali (Pakistan). Other collaborators in imaging in general include T. F. Chan, R. H. Chan, B. Yu, L. Sun, F. L. Yang (China), C. Brito (Mexico), N. Chumchob (Thailand), M. Hintermuller (Germany), Y. Q. Dong (Denmark), X. C. Tai (Norway) etc. [Related publications from http://www.liv.ac.uk/~cmchenke ]

Mathematical imaging is not only a multidisciplinary research area but also a major cross-discipline subject within mathematical sciences as image analysis techniques involve differential geometry, optimization, nonlinear partial differential equations (PDEs), mathematical analysis, computational algorithms and numerical analysis. Segmentation refers to the essential problem in imaging and vision of automatically detecting objects in an image.

In this talk I first review some various models and techniques in the variational framework that are used for segmentation of images, with the purpose of discussing the state of arts rather than giving a comprehensive survey. Then I introduce the practically demanding task of detecting local features in a large image and our recent segmentation methods using energy minimization and PDEs. To ensure uniqueness and fast solution, we reformulate one non-convex functional as a convex one and further consider how to adapt the additive operator splitting method for subsequent solution. Finally I show our preliminary work to attempt segmentation of blurred images in the framework of joint deblurring and segmentation.

This talk covers joint work with Jianping Zhang, Lavdie Rada, Bryan Williams, Jack Spencer (Liverpool, UK), N. Badshah and H. Ali (Pakistan). Other collaborators in imaging in general include T. F. Chan, R. H. Chan, B. Yu, L. Sun, F. L. Yang (China), C. Brito (Mexico), N. Chumchob (Thailand), M. Hintermuller (Germany), Y. Q. Dong (Denmark), X. C. Tai (Norway) etc.

[Related publications from http://www.liv.ac.uk/~cmchenke ]