Past Numerical Analysis Group Internal Seminar

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.

D-finite functions are solutions of linear differential equations with polynomial coefficients.  They have drawn a lot of attention, both in Computer Algebra--because of their numerous (algorithmic) closure properties--but also in Numerical Analysis, because their defining ODEs can be numerically solved very efficiently.  In this talk, I will show how a mix of symbolic and numerical methods yields fast and well-conditioned spectral methods on various domains and using different bases of functions.

It is well-known that a matrix $A$ is Hurwitz stable if and only if there exists a positive definite solution to the Lyapunov matrix equation $A X + X A^* = B$, where $B$ is Hermitian negative definite. We present a verified numerical algorithm to rigorously prove the stability of a given matrix $A$ in the presence of rounding errors.  The computational cost of the algorithm is cubic and it is fast since we can cast almost all operations in level 3 BLAS for which interval arithmetic can be implemented very efficiently.  This is a joint work with Andreas Frommer and the results are already published in ETNA in 2013.

This talk concerns the numerical solution of elliptic partial differential equations posed on general smooth surfaces by the Closest Point Method. Based on the closest point representation of the surface, we formulate an embedding equation in a narrow band surrounding the surface, then discretize it using standard finite differences and interpolation schemes. Numerical convergence of the method will be discussed. In order to solve the resulting large sparse linear systems, we propose a specific geometric multigrid method which makes use of the closest point representation of the surface.
 

We discuss the iterative solution of a finite element discretisation of the magma dynamics equations.  These equations share features of the Stokes equations, however, Elman-Silvester-Wathen (ESW) preconditioners for the magma dynamics equations are not optimal. By introducing a new field, the compaction pressure, into the magma dynamics equations, we have developed a new three-field preconditioner which is optimal in terms of problem size and less sensitive to physical parameters compared to the ESW preconditioners.

There is a well established body of research on quadratic optimization problems based on reformulations of the original problem as a conic program over the cone of completely positive matrices, or its conic dual, the cone of copositive matrices. As a result of this reformulation approach, novel solution schemes for quadratic polynomial optimization problems have been designed by drawing on conic programming tools, and the extensively studied cones of completely positive and of copositive matrices. In particular, this approach has been applied to address key combinatorial optimization problems. Along this line of research, we consider quadratically constrained quadratic programs and provide sufficient and necessary conditions for
this type of problems to be reformulated as a conic program over the cone of completely positive matrices. Thus, recent related results for quadratic problems can be further strengthened. Moreover, these results can be generalized to optimization problems involving higher order polynomias.

A key question to develop our understanding of turbulence in shear flows is: what is the smallest perturbation to the laminar flow that causes a transition to turbulence, and how does this change with the Reynolds number, R?  Finding this so-called ``minimal seed'' is as yet unachievable in direct numerical simulations of the Navier-Stokes equations. We search for the minimal seed in a low-dimensional model analogue to the full Navier-Stokes in plane sinusoidal flow, developed by Waleffe (1997). A previous such calculation found the minimal seed as the least distance (energy norm) from the origin (laminar flow) to the basin of attraction of another fixed point (turbulent attractor).  However, using a non-linear optimization technique, we found an internal boundary of the basin of attraction of the origin that separates flows which directly relaminarize from flows which undergo transient turbulence. It is this boundary which contains the minimal seed, and we find it to be smaller than the previously calculated minimal seed. We present results over a range of Reynolds numbers up to 2000 and find an R^{-1} scaling law fits reasonably well. We propose a new scaling law which asymptotes to R^{-1} for large R but, using some additional information, matches the minimal seed scaling better at low R.

Given a subset E of R^n with empty interior, what is the maximum regularity exponent s for which there exist non-zero distributions in the Bessel potential Sobolev space H^s_p(R^n) that are supported entirely inside E? This question has arisen many times in my recent investigations into boundary integral equation formulations of linear wave scattering by fractal screens, and it is closely related to other fundamental questions concerning Sobolev spaces defined on ``rough'' (i.e. non-Lipschitz) domains. Roughly speaking, one expects that the ``fatter'' the set, the higher the maximum regularity that can be supported. For sets of zero Lebesgue measure one can show, using results on certain set capacities from classical potential theory, that the maximum regularity (if it exists) is negative, and is (almost) characterised by the fractal (Hausdorff) dimension of E. For sets with positive measure the maximum regularity (if it exists) is non-negative,but appears more difficult to characterise in terms of geometrical properties of E.  I will present some partial results in this direction, which have recently been obtained by studying the asymptotic behaviour of the Fourier transform of the characteristic functions of certain fat Cantor sets.

It is well known that piecewise smooth signals are approximately sparse in a wavelet basis. However, other sparse representations are possible, such as the discrete gradient basis. It turns out that signals drawn from a random piecewise constant model have sparser representations in the discrete gradient basis than in Haar wavelets (with high probability). I will talk about this result and its implications, and also show some numerical experiments in which the use of the gradient basis improves compressive signal reconstruction.

One general approach to random number generation is to take a uniformly distributed (0,1) random variable and then invert the cumulative distribution function (CDF) to generate samples from another distribution.  This talk follows this approach, approximating the inverse CDF for the Poisson distribution in a way which is particularly efficient for vector execution on NVIDIA GPUs.

Rational Krylov subspaces have been proven to be useful for many applications, like the approximation of matrix functions or the solution of matrix equations. It will be shown that extended and rational Krylov subspaces —under some assumptions— can be retrieved without any explicit inversion or system solves involved. Instead we do the necessary computations of $A^{-1} v$ in an implicit way using the information from an enlarged standard Krylov subspace.

It is well-known that both for classical and extended Krylov spaces, direct unitary similarity transformations exist providing us the matrix of recurrences. In practice, however, for large dimensions computing time is saved by making use of iterative procedures to gradually gather the recurrences in a matrix. Unfortunately, for extended Krylov spaces one is required to frequently solve, in some way or another a system of equations. In this talk both techniques will be integrated. We start with an orthogonal basis of a standard Krylov subspace of dimension $m+m+p$. Then we will apply a unitary similarity built by rotations compressing thereby significantly the initial subspace and resulting in an orthogonal basis approximately spanning an extended or rational Krylov subspace.

Numerical experiments support our claims that this approximation is very good if the large Krylov subspace contains $A^{-(m+1)} v$, …, $A^{-1} v$ and thus can culminate in nonneglectable dimensionality reduction and as such also can lead to time savings when approximating, e.g., matrix functions.

A stable algorithm to compute the roots of polynomials is presented. The roots are found by computing the eigenvalues of the associated companion matrix by Francis's implicitly-shifted $QR$ algorithm.  A companion matrix is an upper Hessenberg matrix that is unitary-plus-rank-one, that is, it is the sum of a unitary matrix and a rank-one matrix.  These properties are preserved by iterations of Francis's algorithm, and it is these properties that are exploited here. The matrix is represented as a product of $3n-1$ Givens rotators plus the rank-one part, so only $O(n)$ storage space is required.  In fact, the information about the rank-one part is also encoded in the rotators, so it is not necessary to store the rank-one part explicitly.  Francis's algorithm implemented on this representation requires only $O(n)$ flops per iteration and thus $O(n^{2})$ flops overall.  The algorithm is described, backward stability is proved under certain conditions on the polynomial coefficients, and an extensive set of numerical experiments is presented.  The algorithm is shown to be about as accurate as the (slow) Francis $QR$ algorithm applied to the companion matrix without exploiting the structure.  It is faster than other fast methods that have been proposed, and its accuracy is comparable or better.

Traditionally threshold partial pivoting is used to ensure stability of sparse LDL^T factorizations of symmetric matrices. This involves comparing a candidate pivot with all entries in its row/column to ensure that growth in the size of the factors is limited by a threshold at each stage of the factorization. It is capabale of delivering a scaled backwards error on the level of machine precision for practically all real world matrices. However it has significant flaws when used in a massively parallel setting, such as on a GPU or modern supercomputer. It requires all entries of the column to be up-to-date and requires an all-to-all communication for every column. The latter requirement can be performance limiting as the factorization cannot proceed faster than k*(communication latency), where k is the length of the longest path in the sparse elimination tree.

We introduce a new family of communication-avoiding pivoting techniques that reduce the number of messages required by a constant factor allowing the communication cost to be more effectively hidden by computation. We exhibit two members of this family. The first deliver equivalent stability to threshold partial pivoting, but is more pessimistic, leading to additional fill in the factors. The second provides similar fill levels as traditional techniques and, whilst demonstrably unstable for pathological cases, is able to deliver machine accuracy on even the worst real world examples.

The relationship of diagonal dominance ideas to the convergence of stationary iterations is well known. There are a multitude of situations in which such considerations can be used to guarantee convergence when the splitting matrix (the preconditioner) is positive definite. In this talk we will describe and prove sufficient conditions for convergence of a stationary iteration based on a splitting with an indefinite preconditioner. Simple examples covered by this theory coming from Optimization and Economics will be described.

This is joint work with Michael Ferris and Tom Rutherford

This brief lecture will highlight several best-practice observations and
research for writing software for mathematical research, drawn from a
number of sources, including; Best Practices for Scientific Computing
[BestPractices], Code Complete [CodeComplete], and personal observation
from the presenter.  Specific focus will be given to providing the
definition of important concepts, then describing how to apply them
successfully in day-to-day research settings.  Following the outline from
Best Practices, we will cover the following topics:

* Write Programs for People, Not Computers
* Let the Computer Do the Work
* Make Incremental Changes
* Don't Repeat Yourself (or Others)
* Plan for Mistakes
* Optimize Software Only after It Works Correctly
* Document Design and Purpose, Not Mechanics
* Collaborate

[BestPractices]
http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001745
[CodeComplete] http://www.cc2e.com/Default.aspx

We present advances in several fundamental fields of computer vision: image segmentation, object tracking, stereo reconstruction for depth map estimation and full 3D multi-view reconstruction. The basic method applied to these fields is convex relaxation. Convex relaxation methods allow for global optimization of numerous energy functionals and provide a step towards less user input and more automation. We will show how the respective computer vision problems can be formulated in this convex optimization framework. Efficient parallel implementations of the arising numerical schemes using graphics processing units allow for interactive applications.

We investigate an X-ray imaging system that fires multiple point sources of radiation simultaneously from close proximity to a probe. Radiation traverses the probe in a non-parallel fashion, which makes it necessary to use tomosynthesis as a preliminary step to calculating a 2D shadowgraph. The system geometry requires imaging techniques that differ substantially from planar X-rays or CT tomography. We present a proof of concept of such an imaging system, along with relevant artefact removal techniques.  This work is joint with Kishan Patel.

Economic dispatch is a critical part of electricity planning and 
operation. Enhancing the dispatch problem to improve its robustness 
in the face of equipment failures or other contingencies is standard 
practice, but extremely time intensive, leading to restrictions on 
the richness of scenarios considered. We model post-contingency 
corrective actions in the security-constrained economic dispatch 
and consider multiple stages of rescheduling to meet different 
security constraints. The resulting linear program is not solvable
by traditional LP methods due to its large size. We devise and 
implement a series of algorithmic enhancements based on the Benders'
decomposition method to ameliorate the computational difficulty.
In addition, we propose a set of online measures to diagnose
and correct infeasibility issues encountered in the solution process.

The overall solution approach is able to process the ``N-1'' 
contingency list in ten minutes for all large network cases
available for experiments. Extensions to the nonlinear setting will 
be discussed via a semidefinite relaxation.

We present our recent results  on the fluctuation of Optimal Alignments of random sequences and Longest Common Subsequences (LCS). We show how OA and LCS are special cases of certain Last Passage Percolation models which can also be viewed as growth models. this is joint work with Saba Amsalu, Raphael Hauser and Ionel Popescu.

Incomplete Cholesky (IC) factorizations have long been an important tool in the armoury of methods for the numerical solution of large sparse symmetric linear systems Ax = b. In this talk, I will explain the use of intermediate memory (memory used in the construction of the incomplete factorization but is subsequently discarded)  and show how it can significantly improve the performance of the resulting IC preconditioner. I will then focus on extending the approach to sparse symmetric indefinite systems in saddle-point form. A limited-memory signed IC factorization of the form LDLT is proposed, where the diagonal matrix D has entries +/-1. The main advantage of this approach is its simplicity as it avoids the use of numerical pivoting.  Instead, a global shift strategy is used to prevent breakdown and to improve performance. Numerical results illustrate the effectiveness of the signed incomplete Cholesky factorization as a preconditioner.

Ever wondered how the log function in your code is computed? This talk, which was prepared for the 400th anniversary of Napier's development of logarithms, discusses the computation of reciprocals, exponentials and logs, and also my own work on some special functions which are important in Monte Carlo simulation.

Abstract: We propose a term structure power price model that, in contrast to widely accepted no-arbitrage based approaches, accounts for the non-storable nature of power. It belongs to a class of equilibrium game theoretic models with players divided into producers and consumers. Consumers' goal is to maximize a mean-variance utility function subject to satisfying inelastic demand of their own clients (e.g households, businesses etc.) to whom they sell the power on. Producers, who own a portfolio of power plants each defined by a running fuel (e.g. gas, coal, oil...) and physical characteristics (e.g. efficiency, capacity, ramp up/down times, startup costs...), would, similarly, like to maximize a mean-variance utility function consisting of power, fuel, and emission prices subject to production constraints. Our goal is to determine the term structure of the power price at which production matches consumption. In this talk we outline that such a price exists and develop conditions under which it is also unique. Under condition of existence, we propose a tractable quadratic programming formulation for finding the equilibrium term structure of the power price. Numerical results show performance of the algorithm when modeling the whole system of UK power plants.

We study the behaviour of orthogonal polynomials on triangles and their coefficients in the context of spectral approximations of partial differential equations.  For spectral approximation we consider series expansions $u=\sum_{k=0}^{\infty} \hat{u}_k \phi_k$ in terms of orthogonal polynomials $\phi_k$. We show that for any function $u \in C^{\infty}$ the series expansion converges faster than with any polynomial order.  With these result we are able to employ the polynomials $\phi_k$ in the spectral difference method in order to solve hyperbolic conservation laws.

It is a well known fact that discontinuities can arise leading to oscillatory numerical solutions. We compare standard filtering and the super spectral vanishing viscosity methods, which uses exponential filters build from the differential operator of the respective orthogonal polynomials.  We will extend the spectral difference method for unstructured grids by using 
 classical orthogonal polynomials and exponential filters. Finally we consider some numerical test cases.

A fast method for computing eigenpairs of positive definite matrices using GPUs is presented. The method uses Chebyshev polynomial spectral transformations to map the desired eigenvalues of the original matrix $A$ to exterior eigenvalues of the transformed matrix $p(A)$, making them easily computable using existing Krylov methods. The construction of the transforming polynomial $p(z)$ can be done efficiently and only requires knowledge of the spectral radius of $A$. Computing $p(A)v$ can be done using only the action of $Av$. This requires no extra memory and is typically easy to parallelize. The method is implemented using the highly parallel GPU architecture and for specific problems, has a factor of 10 speedup over current GPU methods and a factor of 100 speedup over traditional shift and invert strategies on a CPU.

Partial differential equations defined on surfaces appear in various applications, for example image processing and reconstruction of non-planar images. In this talk, I will present a penalty method for evolution equations, based on an implicit representation of the surface. I will derive a simple equation in the surrounding space, formulated with an extension operator, and then show some analysis and applications of the method.

Several problems lead to the question of how well can a fine grid function be approximated by a coarse grid function, such as preconditioning in finite element methods or data compression and image processing. Particular challenges in answering this question arise when the functions may be only piecewise-continuous, or when the coarse space is not nested in the fine space. In this talk, we solve the problem by using a stable approximation from a space of globally smooth functions as an intermediate step, thereby allowing the use of known approximation results to establish the approximability by a coarse space. We outline the proof, which relies on techniques from the theory of discontinuous Galerkin methods and on the theory of Helmholtz decompositions. Finally, we present an application of our to nonoverlapping domain decomposition preconditioners for hp-version DGFEM.

Given a model dynamical system, a model of any measuring apparatus relating states to observations, and a prior assessment of uncertainty, the probability density of subsequent system states, conditioned upon the history of the observations, is of some practical interest.

When observations are made at discrete times, it is known that the evolving probability density is a solution of the Bayesian filtering equations. This talk will describe the difficulties in approximating the evolving probability density using a Gaussian mixture (i.e. a sum of Gaussian densities). In general this leads to a sequence of optimisation problems and related high-dimensional integrals. There are other problems too, related to the necessity of using a small number of densities in the mixture, the requirement to maintain sparsity of any matrices and the need to compute first and, somewhat disturbingly, second derivatives of the misfit between predictions and observations. Adjoint methods, Taylor expansions, Gaussian random fields and Newton’s method can be combined to, possibly, provide a solution. The approach is essentially a combination of filtering methods and '4-D Var’ methods and some recent progress will be described.

Everybody has heard of the Faraday cage effect, in which a wire mesh does a good job of blocking electric fields and electromagnetic waves. For example, the screen on the front of your microwave oven keeps the microwaves from getting out, while light with its smaller wavelength escapes so you can see your burrito.  Surely the mathematics of such a famous and useful phenomenon has been long ago worked out and written up in the physics books, right?

Well, maybe.   Dave Hewett and I have communicated with dozens of mathematicians, physicists, and engineers on this subject so far, and we've turned up amazingly little.   Everybody has a view of why the Faraday cage mathematics is obvious, and most of their views are different.  Feynman discusses the matter in his Lectures on Physics, but so far as we can tell, he gets it wrong. 

For the static case at least (the Laplace equation), Hewett and I have made good progress with numerical explorations based on  Mikhlin's method backed up by a theorem.   The effect seems to much weaker than we had imagined -- are we missing something?  For time-harmonic waves (the Helmholtz equation), our simulations lead to further puzzles.  We need advice!  Where in the world is the literature on this problem?