Past Numerical Analysis Group Internal Seminar

We provide a simple convergence analysis unified for TR and a new ARC algorithms, which we name ARCq and which is very close in spirit to trust region methods, closer than the original ARC is. We prove global convergence to second order points. We also obtain as a corollary the convergence of the original ARC method. Since one of our aims is to achieve a simple presentation, we sacrifice some generality which we discuss at the end of our developments. In this simplified setting, we prove the optimal complexity property for the ARCq and identify the key elements which allow it. We then propose efficient implementations using a Cholesky like factorization as well as a scalable version based on conjugate gradients.

TBA

In this talk I will give an overview of the algorithms used by Chebfun to numerically compute polynomial and trigonometric minimax approximations of continuous functions. I'll also present Chebfun's capabilities to compute best approximations on compact subsets of an interval and how these methods can be used to design digital filters.

In this talk we will explore the convergence of Krylov methods when used to solve $Lu = f$ where $L$ is an unbounded linear operator.  We will show that for certain problems, methods like Conjugate Gradients and GMRES still converge even though the spectrum of $L$ is unbounded. A theoretical justification for this behavior is given in terms of polynomial approximation on unbounded domains.    

Sparse matrices occur in numerical simulations throughout science and engineering. In particular, it is often desirable to solve systems of the form Ax=b, where A is a sparse matrix with 100,000+ rows and columns. The order that the rows and columns occur in can have a dramatic effect on the viability of a direct solver e.g., the time taken to find x, the amount of memory needed, the quality of x,... We shall consider symmetric matrices and, with the help of playdough, explore how best to order the rows/columns using a nested dissection strategy. Starting with a straightforward strategy, we will discover the pitfalls and develop an adaptive strategy with the aim of coping with a large variety of sparse matrix structures.

Some of the talk will involve the audience playing with playdough, so bring your inner child along with you!

We review a class of systems of non-linear PDEs, derived from the Cahn--Hilliard and Ohta--Kawasaki functionals, that describe the energy evolution of diblock copolymers. These are long chain molecules that can self assemble into repeating patterns as they cool. We are particularly interested in finite element numerical methods that approximate these PDEs in the two-phase (in which we model the polymer only) and three-phase (in which we imagine the polymer surrounded by, and interacting with, a void) cases.

We present a brief derivation of the underlying models, review a class of numerical methods to approximate them, and showcase some early results from our codes.

The talk will present ongoing work on medical image reconstruction from x-ray scanners. A suitable method for reconstruction of these undersampled systems is compressed sensing. The presentation will show respective reconstruction methods and their analysis. Furthermore, work in progress about extensions of the standard approach will be shown.

When the Computing Laboratory discarded its hardcopy journals around 2008, I kept the first ten years or so of each of six classic numerical analysis journals, starting from volume 1, number 1.  This will not be a seminar in the usual sense but a mutual exploration.  Come prepared to look through a few of these old volumes yourself and perhaps to tell the group of something interesting you find.  Bring a pen and paper.  All are welcome.

Mathematics of Computation, from 1943
SIAM Journal, from 1953
Numerische Mathematik, from 1959
BIT, from 1961
SIAM Journal on Numerical Analysis, from 1964
Journal of the IMA, from 1965

In recent years, several preconditioning strategies have been proposed for control problems in fluid dynamics. These are a special case of the general saddle point problem in optimisation. Here, we will extend a preconditionier for distributed Stokes control problems, developed by Rees and Wathen, to the case of boundary control. We will show the usefulness of low-rank structures in constructing a good approximation for the Schur complement of the saddle point matrix. In the end, we will discuss the applicability of this strategy for Navier-Stokes control.

We develop a spectral method for solving univariate singular integral equations over unions of intervals and circles, by utilizing Chebyshev, ultraspherical and Laurent polynomials to reformulate the equations as banded infinite-dimensional systems. Low rank approximations are used to obtain compressed representations of the bivariate kernels. The resulting system can be solved in linear time using an adaptive QR factorization, determining an optimal number of unknowns needed to resolve the solution to any pre-determined accuracy. Applications considered include fracture mechanics, the Faraday cage, and acoustic scattering. The Julia software package https://github.com/ApproxFun/SIE.jl implements our method with a convenient, user-friendly interface.

I present a new method for optimizing quadratic functions subject to simple bound constraints.  If the problem happens to be strictly convex, the algorithm reduces to a highly efficient method by Dostal and Schoberl.  Our algorithm, however, is also able to efficiently solve nonconcex problems. During this talk I will present the algorithm, a sketch of the convergence theory, and numerical results for convex and nonconvex problems.

We will begin by discussing the risk parity portfolio selection problem, which aims to find  portfolios for which the contributions of risk from all assets are equally weighted. The risk parity may be satisfied over either individual assets or groups of assets. We show how convex optimization techniques can find a risk parity solution in the nonnegative  orthant, however, in general cases the number of such solutions can be anywhere between zero and  exponential in the dimension. We then propose a nonconvex least-squares formulation which allows us to consider and possibly solve the general case. 

Motivated by this problem we present several alternating direction schemes for specially structured nonlinear nonconvex problems. The problem structure allows convenient 2-block variable splitting.  Our methods rely on solving convex subproblems at each iteration and converge to a local stationary point. Specifically, discuss approach  alternating directions method of multipliers and the alternating linearization method and we provide convergence rate results for both classes of methods. Moreover, global optimization techniques from polynomial optimization literature are applied to complement our local methods and to provide lower bounds.

Many numerical solvers of ordinary differential equations require problems to be posed as a system of first order differential equations. This means that if one wishes to solve higher order problems, the system have to be rewritten, which is a cumbersome and error-prone process. This talk presents a technique for automatically doing such reformulations.

Choreographies are periodic solutions of the n-body problem in which all of the bodies have unit masses, share a common orbit and are uniformly spread along it. In this talk, I will present an algorithm for numerical computation and stability analysis of choreographies.  It is based on approximations by trigonometric polynomials, minimization of the action functional using a closed-form expression of the gradient, quasi-Newton methods, automatic differentiation and Floquet stability analysis.

The first FFT-based algorithm for numerical homogenization from high-resolution images was proposed by Moulinec and Suquet in 1994 as an alternative to finite elements and twenty years later, it is still widely used in computational micromechanics of materials. The method is based on an iterative solution to an integral equation of the Lippmann-Schwinger type, whose kernel can be explicitly expressed in the Fourier domain. Only recently, it has been recognized that the algorithm has a variational structure arising from a Fourier-Galerkin method. In this talk, I will show how this insight can be used to significantly improve the performance of the original Moulinec-Suquet solver. In particular, I will focus on (i) influence of an iterative solver used to solve the system of linear algebraic equations, (ii) effects of numerical integration of the Galerkin weak form, and (iii) convergence of an a-posteriori bound on the solution during iterations.

A year ago I gave a talk raising questions about Faraday shielding which stimulated discussion with John Ockendon and others and led to a collaboration with Jon Chapman and Dave Hewett.  The problem is one of harmonic functions subject to constant-potential boundary conditions.  A year later, we are happy with the solution we have found, and the paper will appear in SIAM Review.  Though many assume as we originally did that Faraday shielding must be exponentially effective, and Feynman even argues this explicitly in his Lectures, we have found that in fact, the shielding is only linear.  Along the way to explaining this we make use of Mikhlin's numerical method of series expansion, homogenization by multiple scales analysis, conformal mapping, a phase transition, Brownian motion, some ideas recollected from high school about electrostatic induction, and a constrained quadratic optimization problem solvable via a block 2x2 KKT matrix.

When high-frequency acoustic or electromagnetic waves are incident upon an obstacle, the resulting scattered field is composed of rapidly oscillating waves. Conventional numerical methods for such problems use piecewise-polynomial approximation spaces which are not well-suited to capture the oscillatory solution. Hence these methods are prohibitively expensive in the high-frequency regime. Much work has been done in developing “hybrid numerical-asymptotic” (HNA) boundary element methods which utilise approximation spaces containing oscillatory functions carefully chosen to capture the high-frequency asymptotic behaviour of the solution. The computational cost of this approach is significantly smaller than that of conventional methods, and for many problems it is independent of the frequency. In this talk, I will outline the HNA method and discuss its extension to scattering by penetrable obstacles.​

All-at-once schemes aim to solve all time-steps of parabolic PDE-constrained optimization problems in one coupled computation, leading to exceedingly large linear systems requiring efficient iterative methods. We present a new block diagonal preconditioner which is both optimal with respect to the mesh parameter and parallelizable over time, thus can provide significant speed- up. We will present numerical results to demonstrate the effectiveness of this preconditioner.

There is a beautiful problem resulting from arithmetic number theory where a continuous and compactly supported function's 3-fold autoconvolution is constant. In this talk, we optimize the coefficients of a Chebyshev series multiplied by an endpoint singularity to obtain a highly accurate approximation to this constant. Convolving functions with endpoint singularities turns out to be a challenge for standard quadrature routines. However, variable transformations inducing double exponential endpoint decay are used to effectively annihilate the singularities in a way that keeps accuracy high and complexity low.

We present an algorithm, Parallel-$\ell_0$, for {\em combinatorial compressed sensing} (CCS), where the sensing matrix $A \in \mathbb{R}^{m\times n}$ is the adjacency matrix of an expander graph. The information preserving nature of expander graphs allow the proposed algorithm to provably recover a $k$-sparse vector $x\in\mathbb{R}^n$ from $m = \mathcal{O}(k \log (n/k))$ measurements $y = Ax$ via $\mathcal{O}(\log k)$ simple and parallelizable iterations when the non-zeros in the support of the signal form a dissociated set, meaning that all of the $2^k$ subset sums of the support of $x$ are pairwise different. In addition to the low computational cost, Parallel-$\ell_0$ is observed to be able to recover vectors with $k$ substantially larger than previous CCS algorithms, and even higher than $\ell_1$-regularization when the number of measurements is significantly smaller than the vector length.

Flow thought a porous media is usually described by assuming the superficial velocity can be expressed in terms of a constant permeability and a pressure gradient. In poroelastic flows the underlying elastic matrix responds to changes in the fluid pressure. When the elastic deformation is allowed to influence the permeability through the elastic strain, it becomes possible for increased fluid pressure gradient not to result in increased flow, but to decrease the permeability and potentially this may close off or choke the flow. I will talk about a simple model problem for a number of different elastic constitutive models and a number of different permeability-strain models and examine whether there is a general criterion that can be derived to show when, or indeed if, choking can occur for different elasticity-permeability combinations.