Past Numerical Analysis Group Internal Seminar

Composite dilation wavelets are affine systems which extend the notion of wavelets by incorporating a second set of dilations.  The addition of a second set of dilations allows the composite system to capture directional information in addition to time and frequency information.  We classify admissible dilation groups at two extremes: frequency localization through minimally supported frequency composite dilation wavelets and time localization through crystallographic Haar-type composite dilation wavelets. 

Convolution is widely-used and fundamental mathematical operation
in signal processing, statistics, and PDE theory.

Unfortunately the CONV() method in Chebfun for convolving two chebfun 
objects has long been one of the most disappointingly slow features of 
the project. In this talk we will present a new algorithm, which shows 
performance gains on the order of a factor 100.

The key components of the new algorithm are:
* a convolution theorem for Legendre polynomials 
* recurrence relations satisfied by spherical Bessel functions
* recent developments in fast Chebyshev-Legendre transforms [1]

Time-permitting, we shall end with an application from statistics,
using the fact that the probability distribution of the sum of two 
independent random variables is the convolution of their individual 
distributions.

[1] N. Hale and A. Townsend, "A fast, simple, and stable Chebyshev-
Legendre transform using an asymptotic formula”, SISC (to appear).

The accurate and efficient solution of linear systems $Ax=b$ is very important in many engineering and technological applications, and systems of this form also arise as subproblems within other algorithms. In particular, this is true for interior point methods (IPM), where the Newton system must be solved to find the search direction at each iteration. Solving this system is a computational bottleneck of an IPM, and in this talk I will explain how preconditioning and deflation techniques can be used, to lessen this computational burden.  This work is joint with Jacek Gondzio.

A new benchmark, High Performance Conjugate Gradient (HPCG), finally was introduced recently for the Top500 list and the Green500 list. This will draw more attention to performance of sparse iterative solvers on distributed supercomputers and energy efficiency of hardware and software. At the same time, this will more widely promote the concept that communications are the bottleneck of performance of iterative solvers on distributed supercomputers, here we will go a little deeper, discussing components of communications and discuss which part takes a dominate share. Also discussed are mathematics tricks to detect some metrics of an underlying supercomputer.

We discuss the development of finite element techniques and solvers for magma dynamics computations. These are implemented within the FEniCS framework. This approach allows for user-friendly, expressive, high-level code development, but also provides access to powerful, scalable numerical solvers and a large family of finite element discretizations. The ability to easily scale codes to three dimensions with large meshes means that efficiency of the numerical algorithms is vital. We therefore describe our development and analysis of preconditioners designed specifically for finite element discretizations of equations governing magma dynamics. The preconditioners are based on Elman-Silvester-Wathen methods for the Stokes equation, and we extend these to flows with compaction.  This work is joint with Andrew Wathen and Richard Katz from the University of Oxford and Laura Alisic, John Rudge and Garth Wells from the University of Cambridge.

The antitriangular factorisation of real symmetric indefinite matrices recently proposed by Mastronardi and van Dooren has several pleasing properties. It is backward stable, preserves eigenvalues and reveals the inertia, that is, the number of positive, zero and negative eigenvalues. 

In this talk we show that the antitriangular factorization simplifies for saddle point matrices, and that solving a saddle point system in antitriangular form is equivalent to applying the well-known nullspace method. We obtain eigenvalue bounds for the saddle point matrix and discuss the role of the factorisation in preconditioning. 

From cryptography to the proof of Fermat's Last Theorem, elliptic curves (those curves of the form y^2 = x^3 + ax+b) are ubiquitous in modern number theory.  In particular, much activity is focused on developing techniques to discover rational points on these curves. It turns out that finding a rational point on an elliptic curve is very much like finding the proverbial needle in the haystack -- in fact, there is currently no algorithm known to completely determine the group of rational points on an arbitrary elliptic curve.

 I'll introduce the ''real'' picture of elliptic curves and discuss why the ambient real points of these curves seem to tell us little about finding rational points. I'll summarize some of the story of elliptic curves over finite and p-adic fields and tell you about how I study integral points on (hyper)elliptic curves via p-adic integration, which relies on doing a bit of p-adic linear algebra.  Time permitting, I'll also give a short demo of some code we have to carry out these algorithms in the Sage Math Cloud.

In this talk we explore continuous analogues of matrix factorizations.  The analogues we develop involve bivariate functions, quasimatrices (a matrix whose columns are 1D functions), and a definition of triangular in the continuous setting.  Also, we describe why direct matrix algorithms must become iterative algorithms with pivoting for functions. New applications arise for function factorizations because of the underlying assumption of continuity. One application is central to Chebfun2. 

Crouzeix's conjecture is an exasperating problem of linear algebra that has been open since 2004: the norm of p(A) is bounded by twice the maximum value of p on the field of values of A, where A is a square matrix and p is a polynomial (or more generally an analytic function).  I'll say a few words about the conjecture and
show the beautiful proof of Pearcy in 1966 of a special case, based on a vector-valued barycentric interpolation formula.

Hessians of functionals of PDE solutions have important applications in PDE-constrained optimisation (Newton methods) and uncertainty quantification (for accelerating high-dimensional Bayesian inference).  With current techniques, a typical cost for one Hessian-vector product is 4-11 times the cost of the forward PDE solve: such high costs generally make their use in large-scale computations infeasible, as a Hessian solve or eigendecomposition would have costs of hundreds of PDE solves.

In this talk, we demonstrate that it is possible to exploit the common structure of the adjoint, tangent linear and second-order adjoint equations to greatly accelerate the computation of Hessian-vector products, by trading a large amount of computation for a large amount of storage. In some cases of practical interest, the cost of a Hessian-
vector product is reduced to a small fraction of the forward solve, making it feasible to employ sophisticated algorithms which depend on them.

Finding a sparse signal solution of an underdetermined linear system of measurements is commonly solved in compressed sensing by convexly relaxing the sparsity requirement with the help of the l1 norm. Here, we tackle instead the original nonsmooth nonconvex l0-problem formulation using projected gradient methods. Our interest is motivated by a recent surprising numerical find that despite the perceived global optimization challenge of the l0-formulation, these simple local methods when applied to it can be as effective as first-order methods for the convex l1-problem in terms of the degree of sparsity they can recover from similar levels of undersampled measurements. We attempt here to give an analytical justification in the language of asymptotic phase transitions for this observed behaviour when Gaussian measurement matrices are employed. Our approach involves novel convergence techniques that analyse the fixed points of the algorithm and an asymptotic probabilistic analysis of the convergence conditions that derives asymptotic bounds on the extreme singular values of combinatorially many submatrices of the Gaussian measurement matrix under matrix-signal independence assumptions.

This work is joint with Andrew Thompson (Duke University, USA).

We propose a new algorithm for the approximate solution of large-scale high-dimensional tensor-structured linear systems. It can be applied to high-dimensional differential equations, which allow a low-parametric approximation of the multilevel matrix, right-hand side and solution in a tensor product format. We apply standard one-site tensor optimisation algorithm (ALS), but expand the tensor manifolds using the classical iterative schemes (e.g. steepest descent).  We obtain the rank--adaptive algorithm with the theoretical convergence estimate not worse than the one of the steepest descent, and fast practical convergence, comparable or even better than the convergence of more expensive two-site optimisation algorithm (DMRG).
The method is successfully applied for a high--dimensional problem of quantum chemistry, namely the NMR simulation of a large peptide.

This is a joint work with S.Dolgov (Max-Planck Institute, Leipzig, Germany), supported by RFBR and EPSRC grants.

Keywords: high--dimensional problems, tensor train format, ALS, DMRG, steepest descent, convergence rate, superfast algorithms, NMR.

We study a system of partial differential equations describing a steady flow of an incompressible generalized Newtonian fluid, wherein the Cauchy stress depends on concentration. Namely, we consider a coupled system of the generalized Navier-Stokes equations (viscosity of power-law type with concentration dependent power index) and convection-diffusion equation with non-linear diffusivity. We focus on the existence analysis of a weak solution for certain class of models by using a generalization of the monotone operator theory which fits into the framework of generalized Sobolev spaces with variable exponent (class of Sobolev-Orlicz spaces). Such results is then adapted for a suitable FEM approximation, for which the main tool of proof is a generalization of the Lipschitz approximation method.

We present a multilevel preconditioner for the mixed finite element discretization of the biharmonic equation of first kind. While for the interior degrees of freedom a standard multigrid methods can be applied, a different approach is required on the boundary. The construction of the preconditioner is based on a BPX type multilevel representation in fractional Sobolev spaces. Numerical examples illustrate the obtained theoretical results.

Linear wave scattering problems (e.g. for acoustic, electromagnetic and elastic waves) are ubiquitous in science and engineering applications. However, conventional numerical methods for such problems (e.g. FEM or BEM with piecewise polynomial basis functions) are prohibitively expensive when the wavelength of scattered wave is small compared to typical lengthscales of the scatterer (the so-called "high frequency" regime). This is because the solution possesses rapid oscillations which are expensive to capture using conventional approximation spaces. In this talk I will outline some of my recent work in the development of "hybrid numerical-asymptotic" methods, which incur significantly reduced computational cost. These methods use approximation spaces containing oscillatory basis functions, carefully chosen to capture the high frequency asymptotic behaviour. In particular I will discuss some of the interesting challenges arising from non convex, penetrable and three-dimensional scatterers.