Past Numerical Analysis Group Internal Seminar

We investigate theoretical and numerical properties of sparse sketching for both dense and sparse Linear Least Squares (LLS) problems. We show that, sketching with hashing matrices --- with one nonzero entry per column and of size proportional to the rank of the data matrix --- generates a subspace embedding with high probability, provided the given data matrix has low coherence; thus optimal residual values are approximately preserved when the LLS matrix has similarly important rows. We then show that using $s-$hashing matrices, with $s>1$ nonzero entries per column, satisfy similarly good sketching properties for a larger class of low coherence data matrices. Numerically, we introduce our solver Ski-LLS for solving generic dense or sparse LLS problems. Ski-LLS builds upon the successful strategies employed in the Blendenpik and LSRN solvers, that use sketching to calculate a preconditioner before applying the iterative LLS solver LSQR. Ski-LLS significantly improves upon these sketching solvers by judiciously using sparse hashing sketching while also allowing rank-deficiency of input; furthermore, when the data matrix is sparse, Ski-LLS also applies a sparse factorization to the sketched input. Extensive numerical experiments show Ski-LLS is also competitive with other state-of-the-art direct and preconditioned iterative solvers for sparse LLS, and outperforms them in the significantly over-determined regime.

A link for this talk will be sent to our mailing list a day or two in advance.  If you are not on the list and wish to be sent a link, please contact @email.

Riemannian optimization is a powerful and active area of research that studies the optimization of functions defined on manifolds with structure. A class of functions of interest is the set of geodesically convex functions, which are functions that are convex when restricted to every geodesic. In this talk, we will present an accelerated first-order method, nearly achieving the same rates as accelerated gradient descent in the Euclidean space, for the optimization of smooth and g-convex or strongly g-convex functions defined on the hyperbolic space or a subset of the sphere. We will talk about accelerated optimization of another non-convex problem, defined in the Euclidean space, that we solve as a proxy. Additionally, for any Riemannian manifold of bounded sectional curvature, we will present reductions from optimization methods for smooth and g-convex functions to methods for smooth and strongly g-convex functions and vice versa.

This talk is based on the paper https://arxiv.org/abs/2012.03618.

-

A link for this talk will be sent to our mailing list a day or two in advance.  If you are not on the list and wish to be sent a link, please contact @email.

We propose a randomized algorithm for solving a linear system $Ax = b$ with a highly numerically rank-deficient coefficient matrix $A$ with nearly consistent right-hand side possessing a small-norm solution. Our algorithm finds a small-norm solution with small residual in $O(N_r + nrlogn + r^3 )$ operations, where $r$ is the numerical rank of $A$ and $N_r$ is the cost of multiplying an $n\times r$ matrix to $A$. 

Joint work with Marcus Webb (Manchester). 

A link for this talk will be sent to our mailing list a day or two in advance.  If you are not on the list and wish to be sent a link, please contact @email.

Often in scattering applications it is advantageous to reformulate the problem as an integral equation, discretize, and then solve the resulting linear system using a fast direct solver. The computational cost of this approach is typically dominated by the work needed to compress the coefficient matrix into a rank-structured format. In this talk, we present a novel technique which exploits the bandlimited-nature of solutions to the Helmholtz equation in order to accelerate this procedure in environments where multiple frequencies are of interest.

This talk is based on joint work with Gunnar Martinsson (UT Austin).

-

A link for this talk will be sent to our mailing list a day or two in advance.  If you are not on the list and wish to be sent a link, please contact @email.

Motivated by the work of K.R. Rajagopal, the objective of the talk is to discuss the construction and analysis of numerical approximations to a class of models that fall outside the realm of classical Cauchy elasticity. The models under consideration are implicit and nonlinear, and are referred to as strain-limiting, because the linearised strain remains bounded even when the stress is very large, a property that cannot be guaranteed within the framework of classical elastic or nonlinear elastic models. Strain-limiting models can be used to describe, for example, the behavior of brittle materials in the vicinity of fracture tips, or elastic materials in the neighborhood of concentrated loads where there is concentration of stress even though the magnitude of the strain tensor is limited.

We construct a finite element approximation of a strain-limiting elastic model and discuss the theoretical difficulties that arise in proving the convergence of the numerical method. The analytical results are illustrated by numerical experiments.

The talk is based on joint work with Andrea Bonito (Texas A&M University) and Vivette Girault (Sorbonne Université, Paris).

--

A link for this talk will be sent to our mailing list a day or two in advance.  If you are not on the list and wish to be sent a link, please contact @email.

For high order finite elements (continuous piecewise polynomials) the conditioning of the basis is important. However, so far there seems no generally accepted concept of "a well-conditioned basis”,  or a general strategy for how to obtain such representations. In this presentation, we use the $L^2$ condition number as a measure of the conditioning, and construct representations by frames such that the associated $L^2$ condition number is bounded independently of the polynomial degree. The main tools include the bubble transform, which is a stable decomposition of functions into local modes, and orthogonal polynomials on simplexes.  We also include a brief discussion on potential applications in preconditioning. This is a joint work with Ragnar Winther. 

--

A link for this talk will be sent to our mailing list a day or two in advance.  If you are not on the list and wish to be sent a link, please contact @email.

That neural networks may be pruned to high sparsities and retain high accuracy is well established. Recent research efforts focus on pruning immediately after initialization so as to allow the computational savings afforded by sparsity to extend to the training process. In this work, we introduce a new `DCT plus Sparse' layer architecture, which maintains information propagation and trainability even with as little as 0.01% trainable kernel parameters remaining. We show that standard training of networks built with these layers, and pruned at initialization, achieves state-of-the-art accuracy for extreme sparsities on a variety of benchmark network architectures and datasets. Moreover, these results are achieved using only simple heuristics to determine the locations of the trainable parameters in the network, and thus without having to initially store or compute with the full, unpruned network, as is required by competing prune-at-initialization algorithms. Switching from standard sparse layers to DCT plus Sparse layers does not increase the storage footprint of a network and incurs only a small additional computational overhead.

--

A link for this talk will be sent to our mailing list a day or two in advance.  If you are not on the list and wish to be sent a link, please contact @email.

A wide variety of fixed-point iterative methods for the solution of nonlinear operator equations in Hilbert spaces exists. In many cases, such schemes can be interpreted as iterative local linearisation methods, which can be obtained by applying a suitable preconditioning operator to the original (nonlinear) equation. Based on this observation, we will derive a unified abstract framework which recovers some prominent iterative methods. It will be shown that for strongly monotone operators this unified iteration scheme satisfies an energy contraction property. Consequently, the generated sequence converges to a solution of the original problem.

--

A link for this talk will be sent to our mailing list a day or two in advance.  If you are not on the list and wish to be sent a link, please contact @email.

Point cloud registration is the task of finding the transformation that aligns two data sets. We make the assumption that the data lies on a low-dimensional algebraic variety.  The task is phrased as an optimization problem over the special orthogonal group of rotations. We solve this problem using Riemannian optimization algorithms and show numerical examples that illustrate the efficiency of this approach for point cloud registration. 

--

A link for this talk will be sent to our mailing list a day or two in advance.  If you are not on the list and wish to be sent a link, please contact @email.

A link for this talk will be sent to our mailing list a day or two in advance.  If you are not on the list and wish to be sent a link, please contact @email.

Topology optimisation finds the optimal material distribution of a fluid or solid in a domain, subject to PDE, volume, and box constraints. The optimisation problem is normally nonconvex and can support multiple local minima. In recent work [1], the authors developed an algorithm for systematically discovering multiple minima of two-dimensional problems through a combination of barrier methods, active-set strategies, and deflation. The bottleneck of the algorithm is solving the Newton systems that arise. In this talk, we will present preconditioning methods for these linear systems as they occur in the topology optimization of Stokes flow. The strategies involve a mix of block preconditioning and specialized multigrid relaxation schemes that reduce the computational work required and allow the application of the algorithm to three-dimensional problems.

[1] “Computing multiple solutions of topology optimization problems”, I. P. A. Papadopoulos, P. E. Farrell, T. M. Surowiec, 2020, https://arxiv.org/abs/2004.11797

A link for this talk will be sent to our mailing list a day or two in advance.  If you are not on the list and wish to be sent a link, please contact @email.

Identifying discrete patterns in binary data is an important dimensionality reduction tool in machine learning and data mining. In this paper, we consider the problem of low-rank binary matrix factorisation (BMF) under Boolean arithmetic. Due to the NP-hardness of this problem, most previous attempts rely on heuristic techniques. We formulate the problem as a mixed integer linear program and use a large scale optimisation technique of column generation to solve it without the need of heuristic pattern mining. Our approach focuses on accuracy and on the provision of optimality guarantees. Experimental results on real world datasets demonstrate that our proposed method is effective at producing highly accurate factorisations and improves on the previously available best known results for 16 out of 24 problem instances.

--

A link for this talk will be sent to our mailing list a day or two in advance.  If you are not on the list and wish to be sent a link, please contact @email.

In this talk, we will present an approach to constrained optimisation when the set of constraints is a smooth manifold. This setting is of particular interest in data science applications, as many interesting sets of matrices have a manifold structure. We will show how we may couple classic ideas from differential geometry with modern methods such as autodifferentiation to simplify optimisation problems from spaces with a difficult topology (e.g. problems with orthogonal or fixed-rank constraints) to problems on ℝⁿ where we can use any classical optimisation methods to solve them. We will also show how to use these methods to automatically compute quantities such as the Riemannian gradient and Hessian. We will present the library GeoTorch that allows for putting these kind of constraints within models written in PyTorch by adding just one line to the model. We will also comment on some convergence results if time allows.

A link for this talk will be sent to our mailing list a day or two in advance.  If you are not on the list and wish to be sent a link, please contact @email.

A link for this talk will be sent to our mailing list a day or two in advance.  If you are not on the list and wish to be sent a link, please send email to @email.

Preconditioning techniques are widely used for speeding up the iterative solution of systems of linear equations, often by transforming the system into one with lower condition number. Even though the condition number also serves as the determining constant in simple bounds for the numerical error of the solution, simple experiments and bounds show that such preconditioning on the matrix level is not guaranteed to reduce this error. Transformations on the operator level, on the other hand, improve both accuracy and speed of iterative methods as predicted by the change of the condition number. We propose to investigate such methods under a common framework, which we call full operator preconditioning, and show practical examples.

A link for this talk will be sent to our mailing list a day or two in advance.  If you are not on the list and wish to be sent a link, please send an email to @email.

A link for this talk will be sent to our mailing list a day or two in advance.  If you are not on the list and wish to be sent a link, please send email to @email.

A link for this talk will be sent to our mailing list a day or two in advance.  If you are not on the list and wish to be sent a link, please send email to @email.

Smectic A liquid crystals are of great interest in physics for their striking defect structures, including curvature walls and focal conics. However, the mathematical modeling of smectic liquid crystals has not been extensively studied. This work takes a step forward in understanding these fascinating topological defects from both mathematical and numerical viewpoints. In this talk, we propose a new (two- and three-dimensional) mathematical continuum model for the transition between the smectic A and nematic phases, based on a real-valued smectic order parameter for the density perturbation and a tensor-valued nematic order parameter for the orientation. Our work expands on an idea mentioned by Ball & Bedford (2015). By doing so, physical head-to-tail symmetry in half charge defects is respected, which is not possible with vector-valued nematic order parameter.

A link for this talk will be sent to our mailing list a day or two in advance.  If you are not on the list and wish to be sent a link, please send email to @email.

Motivated by the advent of machine learning, the last few years saw the return of hardware-supported low-precision computing. Computations with fewer digits are faster and more memory and energy efficient, but can be extremely susceptible to rounding errors. An application that can largely benefit from the advantages of low-precision computing is the numerical solution of partial differential equations (PDEs), but a careful implementation and rounding error analysis are required to ensure that sensible results can still be obtained. In this talk we study the accumulation of rounding errors in the solution of the heat equation, a proxy for parabolic PDEs, via Runge-Kutta finite difference methods using round-to-nearest (RtN) and stochastic rounding (SR). We demonstrate how to implement the numerical scheme to reduce rounding errors and we present \emph{a priori} estimates for local and global rounding errors. Let $u$ be the roundoff unit. While the worst-case local errors are $O(u)$ with respect to the discretization parameters, the RtN and SR error behaviour is substantially different. We show that the RtN solution is discretization, initial condition and precision dependent, and always stagnates for small enough $\Delta t$. Until stagnation, the global error grows like $O(u\Delta t^{-1})$. In contrast, the leading order errors introduced by SR are zero-mean, independent in space and mean-independent in time, making SR resilient to stagnation and rounding error accumulation. In fact, we prove that for SR the global rounding errors are only $O(u\Delta t^{-1/4})$ in 1D and are essentially bounded (up to logarithmic factors) in higher dimensions.

A link for this talk will be sent to our mailing list a day or two in advance.  If you are not on the list and wish to be sent a link, please send email to @email.

Since Weierstrass it has been known that there are functions that are continuous but nowhere differentiable.  A beautiful example (with probability 1) is any Brownian path.  Brownian paths can be constructed either in space, via Brownian bridge, or in Fourier space, via random Fourier series.

What about functions, which we call "smoothies", that are $C^\infty$ but nowhere analytic?  This case is less familiar but analogous, and again one can do the construction either in space or Fourier space.  We present the ideas and illustrate them with the new Chebfun $\tt{smoothie}$ command.  In the complex plane, the same idea gives functions analytic in the open unit disk and $C^\infty$ on the unit circle, which is a natural boundary.

Reconstruction of 3D images from a set of 2D X-ray projections is a standard inverse problem, particularly in medical imaging. Improvements in imaging technologies have enabled the development of a flat-panel X-ray source, comprised of an array of low-power emitters that are fired in quick succession. During a complete firing sequence, there may be shifts in the patient’s resting position which ultimately create artifacts in the final reconstruction. We present a method for correcting images with respect to unknown body motion, focusing on the case of simple rigid body motion. Image reconstructions are obtained by solving a sparse linear inverse problem, with respect to not only the underlying body but also the unknown velocity. Results find that reconstructions of a moving body can be much better than those obtained by measuring a stationary body, as long as the underlying motion is well approximated.

The standard rounding procedure in floating-point computations is round to nearest (RN). In this talk we consider an alternative rounding strategy called stochastic rounding (SR) which has the appealing property of being exact (actually exact!) in expectation. In the first part of the talk we discuss recent developments in probabilistic rounding error analysis and we show how rounding errors grow at an O(\sqrt{n}) rate rather than O(n) when SR is employed. This shows that Wilkinson's rule of thumb provably holds for this type of rounding. In the second part of the talk we consider the application of SR to parabolic PDEs admitting a steady state solution. We show that when the heat equation is solved in half precision RN fails to compute an accurate solution, while SR successfully solves the problem to decent accuracy.
 

We consider the problem of dynamically pricing multiple products on a network of resources, such as that faced by an airline selling tickets on its route network. For computational reasons this inherently stochastic problem is often approximated deterministically, however even the deterministic dynamic pricing problem can be impractical to solve. For this reason we have derived a novel iterative Quadratic Programming approximation to the deterministic dynamic pricing problem that is not only very fast to solve in practice but also exhibits a provably linear rate of convergence. This is joint work with Saksham Jain and Raphael Hauser.
 

Robust principal component analysis and low-rank matrix completion are extensions of PCA that allow for outliers and missing entries, respectively. Solving these problems requires a low coherence between the low-rank matrix and the canonical basis. However, in both problems the well-posedness issue is even more fundamental; in some cases, both Robust PCA and matrix completion can fail to have any solutions due to the fact that the set of low-rank plus sparse matrices is not closed. Another consequence of this fact is that the lower restricted isometry property (RIP) bound cannot be satisfied for some low-rank plus sparse matrices unless further restrictions are imposed on the constituents. By restricting the energy of one of the components, we close the set and are able to derive the RIP over the set of low rank plus sparse matrices and operators satisfying concentration of measure inequalities. We show that the RIP of an operator implies exact recovery of a low-rank plus sparse matrix is possible with computationally tractable algorithms such as convex relaxations or line-search methods. We propose two efficient iterative methods called Normalized Iterative Hard Thresholding (NIHT) and Normalized Alternative Hard Thresholding (NAHT) that provably recover a low-rank plus sparse matrix from subsampled measurements taken by an operator satisfying the RIP.
 

Randomized SVD has become an extremely successful approach for efficiently computing a low-rank approximation of matrices. In particular the paper by Halko, Martinsson (who is speaking twice this week), and Tropp (SIREV 2011) contains extensive analysis, and made it a very popular method. 
The complexity for $m\times n$ matrices is $O(Nr+(m+n)r^2)$ where $N$ is the cost of a (fast) matrix-vector multiplication; which becomes $O(mn\log n+(m+n)r^2)$ for dense matrices. This work uses classical results in numerical linear algebra to reduce the computational cost to $O(Nr)$ without sacrificing numerical stability. The cost is essentially optimal for many classes of matrices, including $O(mn\log n)$ for dense matrices. The method can also be adapted for updating, downdating and perturbing the matrix, and is especially efficient relative to previous algorithms for such purposes.  

We introduce a new and generic approximation to Schur complements arising from inf-sup stable mixed finite element discretizations of self-adjoint multi-physics problems. The approximation exploits the discretization mesh by forming local, or element, Schur complements and projecting them back to the global degrees of freedom. The resulting Schur complement approximation is sparse, has low construction cost (with the same order of operations as a general finite element matrix), and can be solved using off-the-shelf techniques, such as multigrid. Using the Ladyshenskaja-Babu\v{s}ka-Brezzi condition, we show that this approximation has favorable eigenvalue distributions with respect to the global Schur complement. We present several numerical results to demonstrate the viability of this approach on a range of applications. Interestingly, numerical results show that the method gives an effective approximation to the non-symmetric Schur complement from the steady state Navier-Stokes equations.
 

In STFC's Computational Mathematics Group, we work alongside scientists at large-scale National Facilities, such as ISIS Neutron and Muon source, Diamond Light Source, and the Central Laser Facility. For each of these groups, non-linear least squares fitting is a vital part of their scientific workflow. In this talk I will describe FitBenchmarking, a software tool we have developed to benchmark the performance of different non-linear least squares solvers on real-world data. It is designed to be easily extended, so that new data formats and new minimizers can be added. FitBenchmarking will allow (i) scientists to determine objectively which fitting engine is optimal for solving their problem on their computing architecture, (ii) scientific software developers to quickly test state-of-the-art algorithms in their data analysis software, and (iii) mathematicians and numerical software developers to test novel algorithms against realistic datasets, and to highlight characteristics of problems where the current best algorithms are not sufficient.
 

ARC methods are Newton-type solvers for unconstrained, possibly nonconvex, optimisation problems. In this context, a novel variant based on dynamic inexact Hessian information is discussed. The approach preserves the optimal complexity of the basic framework and the main probabilistic results on the complexity and convergence analysis in the finite-sum minimisation setting will be shown. At the end, some numerical tests on machine learning problems, ill-conditioned databases and real-life applications will be given, in order to confirm the theoretical achievements. Joint work with Stefania Bellavia and Benedetta Morini (Università di Firenze). 

We extend the lightning Laplace solver (Gopal and Trefethen, SINUM 2019) to unbounded domains and to the biharmonic equation. Illustrating the high accuracy of such methods, we get beautiful contour plots of Moffatt eddies.

Many computational techniques in data science involve the factorization of a data matrix into the product of two or more structured matrices. Examples include PCA, which relies on computing an SVD, recommendation systems, which leverage non-negative matrix factorization, infilling missing entries with low rank matrix completion, and finding sparse representations via dictionary learning. In our work we study a new matrix factorization problem, involving the recovery of $\textbf{A}$ and $\textbf{X}$ from $\textbf{Y} := \textbf{A}\textbf{X}$ under the following assumptions; $\textbf{A}$ is an $m \times n$ sparse binary matrix with a fixed number $d$ of nonzeros per column and $\textbf{X}$ is an $n \times N$ sparse real matrix whose columns have $k$ nonzeros and are dissociated. This setup is inspired and motivated by similar models studied in the dictionary learning literature as well as potential connections both with stochastic block models and combinatorial compressed sensing. In this talk we present a new algorithm, EBR, for solving this problem, as well as recovery guarantees in the context of a particular probabilistic data model. Using the properties of expander graphs we are able to show, under certain assumptions, that with just $N = \textit{O}( \log^2(n))$ samples then EBR recovers the factorization up to permutation with high probability. 

We show that the scalability challenges of Global Optimisation algorithms can be overcome for functions with low effective dimensionality, which are constant along certain linear subspaces. Such functions can often be found in applications, for example, in hyper-parameter optimization for neural networks, heuristic algorithms for combinatorial optimization problems and complex engineering simulations. We propose the use of random subspace embeddings within a(ny) global minimisation algorithm, extending the approach in Wang et al. (2016). Using tools from random matrix theory and conic integral geometry, we investigate the efficacy and convergence of our random subspace embeddings approach, in a static and/or adaptive formulation. We illustrate our algorithmic proposals and theoretical findings numerically, using state of the art global solvers. This work is joint with Coralia Cartis.
 

The QDWH algorithm can compute the polar decomposition of a matrix in a stable and efficient way. We generalize this method in order to compute generalized polar decompositions with respect to signature matrices. Here, the role of the QR decomposition is played by the hyperbolic QR decomposition. However, it doesn't show the same favorable properties concerning stability as its orthogonal counterpart. Remedies are found by exploiting connections to the LDL^T factorization and by employing well-conditioned permuted graph bases. The computed polar decomposition is used to formulate a structure-preserving spectral divide-and-conquer method for pseudosymmetric matrices. Applications of this method are found in computational quantum physics, where eigenvalues and eigenvectors describe optical properties of condensed matter or molecules. Additional properties guarantee fast convergence and a reduction to symmetric definite eigenvalue problems after just one step of spectral divide-and-conquer.

Subspace correction (SSC) methods are based on a "divide and conquer" strategy, where a global problem is divided into a sequence of local ones. This framework includes iterative methods such as Jacobi iteration, domain decomposition, and multigrid methods. We are interested in nonlinear PDEs exhibiting highly localized nonlinearities, for which Newton's method can take many iterations. Instead of doing all this global work, nonlinear SSC methods tackle the localized nonlinearities within subproblems. In this presentation, we describe the SSC framework, and investigate combining Newton's method with nonlinear SSC methods to solve a regularized Stefan problem.
 

Vandermonde matrices are exponentially ill-conditioned, rendering the familiar “polyval(polyfit)” algorithm for polynomial interpolation and least-squares fitting ineffective at higher degrees. We show that Arnoldi orthogonalization fixes the problem.

We consider estimation of ordinary differential equation (ODE) models from noisy observations. For this problem, one conventional approach is to fit numerical solutions (e.g., Euler, Runge–Kutta) of ODEs to data. However, such a method does not account for the discretization error in numerical solutions and has limited estimation accuracy. In this study, we develop an estimation method that quantifies the discretization error based on data. The key idea is to model the discretization error as random variables and estimate their variance simultaneously with the ODE parameter. The proposed method has the form of iteratively reweighted least squares, where the discretization error variance is updated with the isotonic regression algorithm and the ODE parameter is updated by solving a weighted least squares problem using the adjoint system. Experimental results demonstrate that the proposed method improves estimation accuracy by accounting for the discretization error in a data-driven manner. This is a joint work with Yuto Miyatake (Osaka University).

A mechanism is described to symmetrize the ultraspherical spectral method for self-adjoint problems. The resulting discretizations are symmetric and banded. An algorithm is presented for an adaptive spectral decomposition of self-adjoint operators. Several applications are explored to demonstrate the properties of the symmetrizer and the adaptive spectral decomposition.

The Alternating Directions Method of Multipliers (ADMM) is a widely popular first-order method for solving convex optimization problems. Its simplicity is arguably one of the main reasons for its popularity. For a broad class of problems, ADMM iterates by repeatedly solving perhaps the two most standard Linear Algebra problems: linear systems and symmetric eigenproblems. In this talk, we discuss how employing standard Krylov-subspace methods for ADMM can lead to tenfold speedups while retaining convergence guarantees.

Subspace methods have the potential to outperform conventional methods, as the derivatives only need to be computed in a smaller dimensional subspace. The sub-problem that needs to be solved at each iteration is also smaller in size, and thus the Linear Algebra cost is also lower. However, if the subspace is not selected "properly", the progress per iteration can be significantly much lower than the progress of the equivalent full-space method. Thus, "improper" selection of the subspace results in subspace methods which are actually more computationally expensive per unit of progress than their full-space alternatives. The popular subspace selection methods (such as randomized) fall into this category when the objective function does not have a known (exploitable) structure. We provide a simple and effective rule to choose the subspace in the "right way" when the objective function does not have a structure. We focus on Gauss-Newton and Least-Squares, but the idea can be generalised to any other solvers and/or objective functions. We show theoretically that the cost of this strategy per unit progress lies in between (approximately) 50% and 100% of the cost of Gauss-Newton, and give an intuition why in practice, it should be closer to the favorable end of the spectrum. We confirm these expectations by running numerical experiments on the CUTEst32 test set. We also compare the proposed selection method with randomized subspace selection. We briefly show that the method is globally convergent and has a 2-step quadratic asymptotic rate of convergence for zero-residual problems.
 

The Approximate Message Passing (AMP) algorithm is a powerful iterative method for reconstructing undersampled sparse signals. Unfortunately, AMP is sensitive to the type of sensing matrix employed and frequently encounters convergence problems. One case where AMP tends to fail is compressed sensing MRI, where Fourier coefficients of a natural image are sampled with variable density. An AMP-inspired algorithm constructed specifically for MRI is presented that exhibits a 'state evolution', where at every iteration the image estimate before thresholding behaves as the ground truth corrupted by Gaussian noise with known covariance. Numerical experiments explore the practical benefits of such effective noise behaviour.
 

Boundary integral equations (BIEs) are well established for solving scattering at bounded infinitely thin objects, so-called screens, which are modelled as “open surfaces” in 3D and as “open curves” in 2D. Moreover, the unknowns of these BIEs are the jumps of traces across $\Gamma$. Things change considerably when considering scattering at multi-screens, which are arbitrary arrangements of thin panels that may not be even locally orientable because of junction points (2D) or junction lines (3D). Indeed, the notion of jumps of traces is no longer meaningful at these junctions. This issue can be solved by switching to a quotient space perspective of traces, as done in recent work by Claeys and Hiptmair. In this talk, we present the extension of the quotient-space approach to the Galerkin boundary element (BE) discretization of first-kind BIEs. Unlike previous approaches, the new quotient-space BEM relies on minimal geometry information and does not require any special treatment at junctions. Moreover, it allows for a rigorous numerical analysis.
 

We describe the main geometric tools required to work on the manifold of fixed-rank symmetric positive-semidefinite matrices: we present expressions for the Riemannian logarithm and the injectivity radius, to complement the already known Riemannian exponential. This manifold is particularly relevant when dealing with low-rank approximations of large positive-(semi)definite matrices. The manifold is represented as a quotient of the set of full-rank rectangular matrices (endowed with the Euclidean metric) by the orthogonal group. Our results allow understanding the failure of some curve fitting algorithms, when the rank of the data is overestimated. We illustrate these observations on a dataset made of covariance matrices characterizing a wind field.

Topology optimisation finds the optimal material distribution of a fluid or solid in a domain, subject to PDE and volume constraints. There are many formulations and we opt for the density approach which results in a PDE, volume and inequality constrained, non-convex, infinite-dimensional optimisation problem without a priori knowledge of a good initial guess. Such problems can exhibit many local minima or even no minima. In practice, heuristics are used to obtain the global minimum, but these can fail even in the simplest of cases. In this talk, we will present an algorithm that solves such problems and systematically discovers as many of these local minima as possible along the way.  

Ocean biogeochemical models used in climate change predictions are very computationally expensive and heavily parameterised. With derivatives too costly to compute, we optimise the parameters within one such model using derivative-free algorithms with the aim of finding a good optimum in the fewest possible function evaluations. We compare the performance of the evolutionary algorithm CMA-ES which is a stochastic global optimization method requiring more function evaluations, to the Py-BOBYQA and DFO-LS algorithms which are local derivative-free solvers requiring fewer evaluations. We also use initial Latin Hypercube sampling to then provide DFO-LS with a good starting point, in an attempt to find the global optimum with a local solver. This is joint work with Coralia Cartis and Samar Khatiwala.
 

The variational data assimilation (VarDA) problem is usually solved using a method equivalent to Gauss-Newton (GN) to obtain the initial conditions for a numerical weather forecast. However, GN is not globally convergent and if poorly initialised, may diverge such as when a long time window is used in VarDA; a desirable feature that allows the use of more satellite data. To overcome this, we apply two globally convergent GN variants (line search & regularisation) to the long window VarDA problem and show when they locate a more accurate solution versus GN within the time and cost available.
Joint work with Coralia Cartis, Amos S. Lawless, Nancy K. Nichols.

Complex spatial patterns such as superlattice patterns and quasipatterns occur in a variety of physical systems ranging from vibrated fluid layers to crystallising soft matter. Reduced order models that describe such systems are usually PDEs. Close to a phase transition, modal expansion along with perturbation methods can be applied to convert the PDEs to normal form equations in the form of coupled ODEs. I use equivariant bifurcation theory along with homotopy methods (developed in computational algebraic geometry) to obtain all solutions of the normal form equations in a non-iterative method. I want to talk about how this approach allows us to ask new questions about the physical systems of interest and what extensions to this method might be possible. This forms a step in my long-term interest to explore how to better ‘complete’ a bifurcation diagram!

We discuss sketching techniques for sparse Linear Least Squares (LLS) problems, that perform a randomised dimensionality reduction for more efficient and scalable solutions. We give theoretical bounds for the accuracy of the sketched solution/residual when hashing matrices are used for sketching, quantifying carefully the trade-off between the coherence of the original, un-sketched matrix and the sparsity of the hashing matrix. We then use these bounds to quantify the success of our algorithm that employs a sparse factorisation of the sketched matrix as a preconditioner for the original LLS, before applying LSQR. We extensively compare our algorithm to state-of-the-art direct and iterative solvers for large-scale and sparse LLS, with encouraging results.

In this talk, I will present a strong (or pathwise) approximation of standard Brownian motion by a class of orthogonal polynomials. Most notably, the coefficients obtained from this expansion are independent Gaussian random variables. This will enable us to generate approximate Brownian paths by matching certain polynomial moments. To conclude the talk, I will discuss related works and applications to numerical methods for SDEs.
 

In this work, we propose a class of numerical schemes for solving semilinear Hamilton-Jacobi-Bellman-Isaacs (HJBI) boundary value problems which arise naturally from exit time problems of diffusion processes with controlled drift. We exploit policy iteration to reduce the semilinear problem into a sequence of linear Dirichlet problems, which are subsequently approximated by a multilayer feedforward neural network ansatz. We establish that the numerical solutions converge globally in the H^2 -norm, and further demonstrate that this convergence is superlinear, by interpreting the algorithm as an inexact Newton iteration for the HJBI equation. Moreover, we construct the optimal feedback controls from the numerical value functions and deduce convergence. The numerical schemes and convergence results are then extended to oblique derivative boundary conditions. Numerical experiments on the stochastic Zermelo navigation problem and the perpetual American option pricing problems are presented to illustrate the theoretical results and to demonstrate the effectiveness of the method.
 

Partial differential equations governing unknown or evolving geometries are ubiquitous in applications and challenging to discretize. A great deal of numerical work in this area has focused on extrinsic geometric flows, where the evolving geometry is a curve or surface embedded in Euclidean space. Much less attention has been paid to the discretization of intrinsic geometric flows, where the evolving geometry is described by a Riemannian metric. This talk will present finite element discretizations of such flows.