Past Computational Mathematics and Applications Seminar

Because of atmospheric turbulence, images of objects in outer space acquired via ground-based telescopes are usually blurry.  One way to estimate the blurring kernel or point spread function (PSF) is to make use of the aberration of wavefront received at the telescope, i.e., the phase. However only the low-resolution wavefront gradients can be collected by wavefront sensors. In this talk, I will discuss how to use regularization methods to reconstruct high-resolution phase gradients and then use them to recover the phase and the PSF in high accuracy. I will end by relating the problem to high-resolution image reconstruction and methods for solving it.
Joint work with Rui Zhao and research supported by HKRGC.

Several optimization problems combine nonlinear constraints with the integrality of a subset of variables. For an important class of problems  called Mixed Integer Second-Order Cone Optimization (MISOCO), with applications in facility location, robust optimization, and finance, among others, these nonlinear constraints are second-order (or Lorentz) cones.

For such problems, as for many discrete optimization problems, it is crucial to understand the properties of the union of two disjoint sets of feasible solutions. To this end, we apply the disjunctive programming paradigm to MISOCO and present conditions under which the convex hull of two disjoint sets can be obtained by intersecting the feasible set with a specially constructed second-order cone. Computational results show that such cone has a positive impact on the solution of MISOCO problems.

The mathematical analysis and numerical simulation of acoustic and electromagnetic wave scattering by planar screens is a classical topic. The standard technique involves reformulating the problem as a boundary integral equation on the screen, which can be solved numerically using a boundary element method. Theory and computation are both well-developed for the case where the screen is an open subset of the plane with smooth (e.g. Lipschitz or smoother) boundary. In this talk I will explore the case where the screen is an arbitrary subset of the plane; in particular, the screen could have fractal boundary, or itself be a fractal. Such problems are of interest in the study of fractal antennas in electrical engineering, light scattering by snowflakes/ice crystals in atmospheric physics, and in certain diffraction problems in laser optics. The roughness of the screen presents challenging questions concerning how boundary conditions should be enforced, and the appropriate function space setting. But progress is possible and there is interesting behaviour to be discovered: for example, a sound-soft screen with zero area (planar measure zero) can scatter waves provided the fractal dimension of the set is large enough. Accurate computations are also challenging because of the need to adapt the mesh to the fine structure of the fractal. As well as presenting numerical results, I will outline some of the outstanding open questions from the point of view of numerical analysis. This is joint work with Simon Chandler-Wilde (Reading) and Andrea Moiola (Pavia).
 

In this talk, I will present a novel solution strategy to efficiently and accurately compute approximate solutions to semilinear optimal control problems, focusing on the optimal control of phase field formulations of geometric evolution laws.
The optimal control of geometric evolution laws arises in a number of applications in fields including material science, image processing, tumour growth and cell motility.
Despite this, many open problems remain in the analysis and approximation of such problems.
In the current work we focus on a phase field formulation of the optimal control problem, hence exploiting the well developed mathematical theory for the optimal control of semilinear parabolic partial differential equations.
Approximation of the resulting optimal control problem is computationally challenging, requiring massive amounts of computational time and memory storage.
The main focus of this work is to propose, derive, implement and test an efficient solution method for such problems. The solver for the discretised partial differential equations is based upon a geometric multigrid method incorporating advanced techniques to deal with the nonlinearities in the problem and utilising adaptive mesh refinement.
An in-house two-grid solution strategy for the forward and adjoint problems, that significantly reduces memory requirements and CPU time, is proposed and investigated computationally.
Furthermore, parallelisation as well as an adaptive-step gradient update for the control are employed to further improve efficiency.
Along with a detailed description of our proposed solution method together with its implementation we present a number of computational results that demonstrate and evaluate our algorithms with respect to accuracy and efficiency.
A highlight of the present work is simulation results on the optimal control of phase field formulations of geometric evolution laws in 3-D which would be computationally infeasible without the solution strategies proposed in the present work.

Adjoint derivatives reveal the sensitivity of a computer program's output to changes in its inputs. These derivatives are useful as a building block for optimisation, uncertainty quantification, noise estimation, inverse design, etc., in many industrial and scientific applications that use PDE solvers or other codes.
Algorithmic differentiation (AD) is an established method to transform a given computation into its corresponding adjoint computation. One of the key challenges in this process is the efficiency of the resulting adjoint computation. This becomes especially pressing with the increasing use of shared-memory parallelism on multi- and many-core architectures, for which AD support is currently insufficient.
In this talk, I will present an overview of challenges and solutions for the differentiation of shared-memory-parallel code, using two examples: an unstructured-mesh CFD solver, and a structured-mesh stencil kernel, both parallelised with OpenMP. I will show how AD can be used to generate adjoint solvers that scale as well as their underlying original solvers on CPUs and a KNC XeonPhi. The talk will conclude with some recent efforts in using AD and formal verification tools to check the correctness of manually optimised adjoint solvers.
 

Gauss invented Gaussian quadrature following an approach entirely different from the one we now find in textbooks. I will describe leisurely the contents of Gauss's original memoir on quadrature, an impressive piece of mathematics, based on continued fractions, Padé approximation, generating functions, the hypergeometric series and more.

The talk will describe accelerated algorithms for computing full or partial matrix factorizations such as the eigenvalue decomposition, the QR factorization, etc. The key technical novelty is the use of  randomized projections to reduce the effective dimensionality of  intermediate steps in the computation. The resulting algorithms execute faster on modern hardware than traditional algorithms, and are particularly well suited for processing very large data sets.

The algorithms described are supported by a rigorous mathematical analysis that exploits recent work in random matrix theory. The talk will briefly review some representative theoretical results.

For many applications in science and engineering, the ability to efficiently and accurately approximate solutions to elliptic PDEs dictates what physical phenomena can be simulated numerically.  In this seminar, we present a high-order accurate discretization technique for variable coefficient PDEs with smooth coefficients.  The technique comes with a nested dissection inspired direct solver that scales linearly or nearly linearly with respect to the number of unknowns.  Unlike the application of nested dissection methods to classic discretization techniques, the constant prefactors do not grow with the order of the discretization.  The discretization is robust even for problems with highly oscillatory solutions.  For example, a problem 100 wavelengths in size can be solved to 9 digits of accuracy with 3.7 million unknowns on a desktop computer.  The precomputation of the direct solver takes 6 minutes on a desktop computer.  Then applying the computed solver takes 3 seconds.  The recent application of the algorithm to inverse media scattering also will be presented.

Structural optimization can be interpreted as the attempt to find the best mechanical structure to support specific load cases respecting some possible constraints. Within this context, topology optimization aims to obtain the connectivity, shape and location of voids inside a prescribed structural design domain. The methods for the design of stiff lightweight structures are well established and can already be used in a specific range of industries where such structures are important, e.g., in aerospace and automobile industries.

In this seminar, we will go through the basic engineering concepts used to quantify and analyze the computational models of mechanical structures. After presenting the motivation, the methods and mathematical tools used in structural topology optimization will be discussed. In our method, an implicit level set function is used to describe the structural boundaries. The optimization problem is approximated by linearization of the objective and constraint equations via Taylor’s expansion. Shape sensitivities are used to evaluate the change in the structural performance due to a shape movement and to feed the mathematical optimiser in an iterative procedure. Recent developments comprising multiscale and Multiphysics problems will be presented and a specific application proposal including acoustic-structure interaction will be discussed.

We describe a convergence acceleration technique for generic optimization problems. Our scheme computes estimates of the optimum from a nonlinear average of the iterates produced by any optimization method. The weights in this average are computed via a simple linear system, whose solution can be updated online. This acceleration scheme runs in parallel to the base algorithm, providing improved estimates of the solution on the fly, while the original optimization method is running. Numerical experiments are detailed on classical classification problems.
 

Abstract: We study nonuniform sampling in shift-invariant spaces whose generator is a totally positive function. For a subclass of such generators the sampling theorems can be formulated in analogy to the theorems of Beurling and Landau for bandlimited functions. These results are  optimal and validate  the  heuristic reasonings in the engineering literature. In contrast to the cardinal series, the reconstruction procedures for sampling in a shift-invariant space with a totally positive generator  are local and thus accessible to numerical linear algebra.

A subtle  connection between sampling in shift-invariant spaces and the theory of Gabor frames leads to new and optimal  results for Gabor frames.  We show that the set of phase-space shifts of  $g$ (totally positive with a Gaussian part) with respect to a rectangular lattice forms a frame, if and only if the density of the lattice  is strictly larger than 1. This solves an open problem going backto Daubechies in 1990 for the class of totally positive functions of Gaussian type.
 

Almost all real-world applications involve a degree of uncertainty. This may be the result of noisy measurements, restrictions on observability, or simply unforeseen events. Since many models in both engineering and the natural sciences make use of partial differential equations (PDEs), it is natural to consider PDEs with random inputs. In this context, passing from modelling and simulation to optimization or control results in stochastic PDE-constrained optimization problems. This leads to a number of theoretical, algorithmic, and numerical challenges.

 From a mathematical standpoint, the solution of the underlying PDE is a random field, which in turn makes the quantity of interest or the objective function an implicitly defined random variable. In order to minimize this distributed objective, one can use, e.g., stochastic order constraints, a distributionally robust approach, or risk measures. In this talk, we will make use of risk measures.

After motivating the approach via a model for the mitigation of an airborne pollutant, we build up an analytical framework and introduce some useful risk measures. This allows us to prove the existence of solutions and derive optimality conditions. We then present several approximation schemes for handling non-smooth risk measures in order to leverage existing numerical methods from PDE-constrained optimization. Finally, we discuss solutions techniques and illustrate our results with numerical examples.

During the last decade, the progress in the computational performance of commercial mixed-integer programming solvers have been significant. Part of this success is due to faster computers and better software engineering but a more significant part of it is due to the power of the cutting planes used in these solvers.
In the first part of this talk, we will discuss main components of a MIP solver and describe some classical families of valid inequalities (Gomory mixed integer cuts, mixed integer rounding cuts, split cuts, etc.) that are routinely used in these solvers. In the second part, we will discuss recent progress in cutting plane theory that has not yet made its way to commercial solvers. In particular, we will discuss cuts from lattice-free convex sets and answer a long standing question in the affirmative by deriving a finite cutting plane algorithm for mixed-integer programming.

In variational imaging and other inverse problem modeling, regularisation plays a major role. In recent years, high order regularizers such as the total generalised variation, the mean curvature and the Gaussian curvature are increasingly studied and applied, and many improved results over the widely-used total variation model are reported.
Here we first introduce the fractional order derivatives and the total fractional-order variation which provides an alternative  regularizer and is not yet formally analysed. We demonstrate that existence and uniqueness properties of the new model can be analysed in a fractional BV space, and, equally, the new model performs as well as the high order regularizers (which do not yet have much theory). 
In the usual framework, the algorithms of a fractional order model are not fast due to dense matrices involved. Moreover, written in a Bregman framework, the resulting Sylvester equation with Toeplitz coefficients can be solved efficiently by a preconditioned solver. Further ideas based on adaptive integration can also improve the computational efficiency in a dramatic way.
 Numerical experiments will be given to illustrate the advantages of the new regulariser for both restoration and registration problems.
 

We will present a very general framework for unconstrained stochastic optimization which is based on standard trust region framework using  random models. In particular this framework retains the desirable features such step acceptance criterion, trust region adjustment and ability to utilize of second order models. We make assumptions on the stochasticity that are different from the typical assumptions of stochastic and simulation-based optimization. In particular we assume that our models and function values satisfy some good quality conditions with some probability fixed, but can be arbitrarily bad otherwise. We will analyze the convergence and convergence rates of this general framework and discuss the requirement on the models and function values. We will will contrast our results with existing results from stochastic approximation literature. We will finish with examples of applications arising the area of machine learning. 
 

Work with Jemima Tabeart, Sarah Dance, Nancy Nichols, Joanne Waller (University of Reading) and Stefano Migliorini, Fiona Smith (Met Office). 
In environmental prediction variational data assimilation (DA) is a method for using observational data to estimate the current state of the system. The DA problem is usually solved as a very large nonlinear least squares problem, in which the fit to the measurements is balanced against the fit to a previous model forecast. These two terms are weighted by matrices describing the correlations of the errors in the forecast and in the observations. Until recently most operational weather and ocean forecasting systems assumed that the errors in the observations are uncorrelated. However, as we move to higher resolution observations then it is becoming more important to specify observation error correlations. In this work we look at the effect this has on the conditioning of the optimization problem. In the context of a linear system we develop bounds on the condition number of the problem in the presence of correlated observation errors. We show that the condition number is very dependent on the minimum eigenvalue of the observation error correlation matrix. We then present results using the Met Office data assimilation system, in which different methods for reconditioning the correlation matrix are tested. We investigate the effect of these different methods on the conditioning and the final solution of the problem.
 

Inverse problems are ubiquitous in many areas of Science and Engineering and, once discretised, they lead to ill-conditioned linear systems, often of huge dimensions: regularisation consists in replacing the original system by a nearby problem with better numerical properties, in order to find meaningful approximations of its solution. In this talk we will explore the regularisation properties of many iterative methods based on Krylov subspaces. After surveying some basic methods such as CGLS and GMRES, innovative approaches based on flexible variants of CGLS and GMRES will be presented, in order to efficiently enforce nonnegativity and sparsity into the solution.

We consider the Cauchy (or steepest descent) method with exact line search applied to a strongly convex function with Lipschitz continuous gradient. We establish the exact worst-case rate of convergence of this scheme, and show that this worst-case behavior is exhibited by a certain convex quadratic function. We also give worst-case complexity bound for a noisy variant of gradient descent method. Finally, we show that these results may be applied to study the worst-case performance of Newton's method for the minimization of self-concordant functions.

The proofs are computer-assisted, and rely on the resolution of semidefinite programming performance estimation problems as introduced in the paper [Y. Drori and M. Teboulle.  Performance of first-order methods for smooth convex minimization: a novel approach. Mathematical Programming, 145(1-2):451-482, 2014].

Joint work with F. Glineur and A.B. Taylor.

 Functions defined by evaluation programs involving smooth  elementals and absolute values as well as max and min are piecewise smooth. For this class we present first and second order, necessary and sufficient conditions for the functions to be locally optimal, or convex, or at least possess a supporting hyperplane. The conditions generalize the classical KKT and SSC theory and are constructive; though in the case of convexity they may be combinatorial to verify. As a side product we find that, under the Mangasarin-Fromowitz-Kink-Qualification, the well established nonsmooth concept of subdifferential regularity is equivalent to first order convexity. All results are based on piecewise linearization and suggest corresponding optimization algorithms.

We propose a multilevel paradigm for the global optimisation of polynomials with sparse support. Such polynomials arise through the discretisation of PDEs, optimal control problems and in global optimization applications in general. We construct projection operators to relate the primal and dual variables of the SDP relaxation between lower and higher levels in the hierarchy, and theoretical results are proven to confirm their usefulness. Numerical results are presented for polynomial problems that show how these operators can be used in a hierarchical fashion to solve large scale problems with high accuracy.

I will present a broad family of stochastic algorithms for inverting a matrix, including specialized variants which maintain symmetry or positive definiteness of the iterates. All methods in the family converge globally and linearly, with explicit rates. In special cases, the methods obtained are stochastic block variants of several quasi-Newton updates, including bad Broyden (BB), good Broyden (GB), Powell-symmetric-Broyden (PSB), Davidon-Fletcher-Powell (DFP) and Broyden-Fletcher-Goldfarb-Shanno (BFGS). After a pause for questions, I will then present a block stochastic BFGS method based on the stochastic method for inverting positive definite matrices. In this method, the estimate of the inverse Hessian matrix that is maintained by it, is updated at each iteration using a sketch of the Hessian, i.e., a randomly generated compressed form of the Hessian. I will propose several sketching strategies, present a new quasi-Newton method that uses stochastic block BFGS updates combined with the variance reduction approach SVRG to compute batch stochastic gradients, and prove linear convergence of the resulting method. Numerical tests on large-scale logistic regression problems reveal that our method is more robust and substantially outperforms current state-of-the-art methods.

Most current methods of Magnetic Resonance Imaging (MRI) reconstruction interpret raw signal values as samples of the Fourier transform of the object. Although this is computationally convenient, it neglects relaxation and off–resonance evolution in phase, both of which can occur to significant extent during a typical MRI signal. A more accurate model, known as Parameter Assessment by Recovery from Signal Encoding (PARSE), takes the time evolution of the signal into consideration. This model uses three parameters that depend on tissue properties: transverse magnetization, signal decay rate, and frequency offset from resonance. Two difficulties in recovering an image using this model are the low SNR for long acquisition times in single-shot MRI, and the nonlinear dependence of the signal on the decay rate and frequency offset. In this talk, we address the latter issue by using a second order approximation of the original PARSE model. The linearized model can be solved using convex optimization augmented with well-stablished regularization techniques such as total variation. The sensitivity of the parameters to noise and computational challenges associated with this approximation will be discussed.

We consider the problem of computing a nonnegative low rank factorization to a given nonnegative input matrix under the so-called "separabilty condition".  This assumption makes this otherwise NP hard problem polynomial time solvable, and we will use first order optimization techniques to compute such a factorization. The optimization model use is based on sparse regression with a self-dictionary, in which the low rank constraint is relaxed to the minimization of an l1-norm objective function.  We apply these techniques to endmember detection and classification in hyperspecral imaging data.

 The matrix logarithm, when applied to symmetric positive definite matrices, is known to satisfy a notable concavity property in the positive semidefinite (Loewner) order. This concavity property is a cornerstone result in the study of operator convex functions and has important applications in matrix concentration inequalities and quantum information theory.
In this talk I will show that certain rational approximations of the matrix logarithm remarkably preserve this concavity property and moreover, are amenable to semidefinite programming. Such approximations allow us to use off-the-shelf semidefinite programming solvers for convex optimization problems involving the matrix logarithm. These approximations are also useful in the scalar case and provide a much faster alternative to existing methods based on successive approximation for problems involving the exponential/relative entropy cone. I will conclude by showing some applications to problems arising in quantum information theory.

This is joint work with James Saunderson (Monash University) and Pablo Parrilo (MIT)

Rational Krylov methods are applicable to a wide range of scientific computing problems, and ​the rational Arnoldi algorithm is a commonly used procedure for computing an ​orthonormal basis of a rational Krylov space. Typically, the computationally most expensive component of this​ ​algorithm is the solution of a large linear system of equations at each iteration. We explore the​ ​option of solving several linear systems simultaneously, thus constructing the rational Krylov​ ​basis in parallel. If this is not done carefully, the basis being orthogonalized may become badly​ ​conditioned, leading to numerical instabilities in the orthogonalization process. We introduce the​ ​new concept of continuation pairs which gives rise to a near-optimal parallelization strategy that ​allows to control the growth of the condition number of this nonorthogonal basis. As a consequence we obtain a significantly more accurate and reliable parallel rational Arnoldi algorithm.
​ ​
The computational benefits are illustrated using several numerical examples from different application areas.
​ ​
This ​talk is based on joint work with Mario Berljafa  available as an Eprint at http://eprints.ma.man.ac.uk/2503/
 

We present global rates of convergence for a general class of methods for nonconvex smooth optimization that include linesearch, trust-region and regularisation strategies, but that allow inaccurate problem information. Namely, we assume the local (first- or second-order) models of our function are only sufficiently accurate with a certain probability, and they can be arbitrarily poor otherwise. This framework subsumes certain stochastic gradient analyses and derivative-free techniques based on random sampling of function values. It can also be viewed as a robustness
assessment of deterministic methods and their resilience to inaccurate derivative computation such as due to processor failure in a distribute framework. We show that in terms of the order of the accuracy, the evaluation complexity of such methods is the same as their counterparts that use deterministic accurate models; the use of probabilistic models only increases the complexity by a constant, which depends on the probability of the models being good. Time permitting, we also discuss the case of inaccurate, probabilistic function value information, that arises in stochastic optimization. This work is joint with Katya Scheinberg (Lehigh University, USA).
 

For linear dynamical systems, model reduction has achieved great success. In the case of linear dynamics,  we know how to construct, at a modest cost, (locally) optimal, input-independent reduced models; that is, reduced models that are uniformly good over all inputs having bounded energy. In addition, in some cases we can achieve this goal using only input/output data without a priori knowledge of internal  dynamics.  Even though model reduction has been successfully and effectively applied to nonlinear dynamical systems as well, in this setting,  bot the reduction process and the reduced models are input dependent and the high fidelity of the resulting approximation is generically restricted to the training input/data. In this talk, we will offer remedies to this situation.

To predict the behaviour of a dynamical system using a mathematical model, an accurate estimate of the current state of the system is needed in order to initialize the model. Complete information on the current state is, however, seldom available. The aim of optimal state estimation, known in the geophysical sciences as ‘data assimilation’, is to determine a best estimate of the current state using measured observations of the real system over time, together with the model equations. The problem is commonly formulated in variational terms as a very large nonlinear least-squares optimization problem. The lack of complete data, coupled with errors in the observations and in the model, leads to a highly ill-conditioned inverse problem that is difficult to solve.

To understand the nature of the inverse problem, we examine how different components of the assimilation system influence the conditioning of the optimization problem. First we consider the case where the dynamical equations are assumed to model the real system exactly. We show, against intuition, that with increasingly dense and precise observations, the problem becomes harder to solve accurately. We then extend these results to a 'weak-constraint' form of the problem, where the model equations are assumed not to be exact, but to contain random errors. Two different, but mathematically equivalent, forms of the problem are derived. We investigate the conditioning of these two forms and find, surprisingly, that these have quite different behaviour.

The requirement to compute Jordan blocks for multiple eigenvalues arises in a number of physical problems, for example panel flutter problems in aerodynamical stability, the stability of electrical power systems, and in quantum mechanics. We introduce a general method for computing a 2-dimensional Jordan block in a parameter-dependent matrix eigenvalue problem based on the so called Implicit Determinant Method. This is joint work with Alastair Spence (Bath).

In many applications we need to find or estimate the $p \ge 1$ largest elements of a matrix, along with their locations. This is required for recommender systems used by Amazon and Netflix, link prediction in graphs, and in finding the most important links in a complex network, for example. 

Our algorithm uses only matrix vector products and is based upon a power method for mixed subordinate norms. We have obtained theoretical results on the convergence of this algorithm via a comparison with rook pivoting for the LU  decomposition. We have also improved the practicality of the algorithm by producing a blocked version iterating on $n \times t$ matrices, as opposed to vectors, where $t$ is a tunable parameter. For $p > 1$ we show how deflation can be used to improve the convergence of the algorithm. 

Finally, numerical experiments on both randomly generated matrices and real-life datasets (the latter for $A^TA$ and $e^A$) show how our algorithms can reliably estimate the largest elements of a matrix whilst obtaining considerable speedups when compared to forming the matrix explicitly: over 1000x in some cases.

For many years, sum-factored algorithms for finite elements in rectangular reference geometry have combined low complexity with the mathematical power of high-order approximation.  However, such algorithms rely heavily on the tensor product structure inherent in the geometry and basis functions, and similar algorithms for simplicial geometry have proven elusive.

Bernstein polynomials are totally nonnegative, rotationally symmetric, and geometrically decomposed bases with many other remarkable properties that lead to optimal-complexity algorithms for element wise finite element computations.  The also form natural building blocks for the finite element exterior calculus bases for the de Rham complex so that H(div) and H(curl) bases have efficient representations as well.  We will also their relevance for explicit discontinuous Galerkin methods, where the element mass matrix requires special attention.

We will review the basic building blocks of iterative solvers, i.e. sparse matrix-vector multiplication, in the context of GPU devices such 
as the cards by NVIDIA; we will then discuss some techniques in preconditioning by approximate inverses, and we will conclude with an 
application to an image processing problem from the biomedical field.

The aim of this talk is to design an efficient multigrid method for constrained convex optimization problems arising from discretization  of  some  underlying  infinite  dimensional  problems. Due  to problem  dependency  of this approach, we only consider bound constraints with (possibly) a linear equality constraint. As our aim is to target large-scale problems, we want to avoid computation of second 
derivatives of the objective function, thus excluding Newton like methods. We propose a smoothing operator that only uses first-order information and study the computational efficiency of the resulting method. In the second part, we consider application of multigrid techniques to more general optimization problems, in particular, the topology design problem.

Quadrature is the term for the numerical evaluation of integrals.  It's a beautiful subject because it's so accessible, yet full of conceptual surprises and challenges.  This talk will review ten of these, with plenty of history and numerical demonstrations.  Some are old if not well known, some are new, and two are subjects of my current research.

Partial differential equations with more than three coordinates arise naturally if the model features certain kinds of stochasticity. Typical examples are the Schroedinger, Fokker-Planck and Master equations in quantum mechanics or cell biology, as well as quantification of uncertainty.
The principal difficulty of a straightforward numerical solution of such equations is the `curse of dimensionality': the storage cost of the discrete solution grows exponentially with the number of coordinates (dimensions).

One way to reduce the complexity is the low-rank separation of variables. One can see all discrete data (such as the solution) as multi-index arrays, or tensors. These large tensors are never stored directly.
We approximate them by a sum of products of smaller factors, each carrying only one of the original variables. I will present one of the simplest but powerful of such representations, the Tensor Train (TT) decomposition. The TT decomposition generalizes the approximation of a given matrix by a low-rank matrix to the tensor case. It was found that many interesting models allow such approximations with a significant reduction of storage demands.

A workhorse approach to computations with the TT and other tensor product decompositions is the alternating optimization of factors. The simple realization is however prone to convergence issues.
I will show some of the recent improvements that are indispensable for really many dimensions, or solution of linear systems with non-symmetric or indefinite matrices.

To face the advent of multicore processors and the ever increasing complexity of hardware architectures, programming
models based on DAG parallelism regained popularity in the high performance, scientific computing community. Modern runtime systems offer a programming interface that complies with this paradigm and powerful engines for scheduling the tasks into which the application is decomposed. These tools have already proved their effectiveness on a number of dense linear algebra applications. 

In this talk we present the design of task-based sparse direct solvers on top of runtime systems. In the context of the
qr_mumps solver, we prove the usability and effectiveness of our approach with the implementation of a sparse matrix multifrontal factorization based on a Sequential Task flow parallel programming model. Using this programming model, we developed features such as the integration of dense 2D Communication Avoiding algorithms in the multifrontal method allowing for better scalability compared to the original approach used in qr_mumps.

Following this approach, we move to heterogeneous architectures where task granularity and scheduling strategies are critical to achieve performance. We present, for the multifrontal method, a hierarchical strategy for data partitioning and a scheduling algorithm capable of handling the heterogeneity of resources.   Finally we introduce a memory-aware algorithm to control the memory behavior of our solver and show, in the context of multicore architectures, an important reduction of the memory footprint for the multifrontal QR factorization with a small impact on performance.

Functions are usually approximated numerically in a basis, a non-redundant and complete set of functions that span a certain space. In this talk we highlight a number of benefits of using overcomplete sets, in particular using the more general notion of a "frame". The main 

benefit is that frames are easily constructed even for functions of several variables on domains with irregular shapes. On the other hand, allowing for possible linear depencies naturally leads to ill-conditioning of approximation algorithms. The ill-conditioning is 

potentially severe. We give some useful examples of frames and we first address the numerical stability of best approximations in a frame. Next, we briefly describe special point sets in which interpolation turns out to be stable. Finally, we review so-called Fourier extensions and an efficient algorithm to approximate functions with spectral accuracy on domains without structure.

When assigned with the task of extracting information from given image data the first challenge one faces is the derivation of a truthful model for both the information and the data. Such a model can be determined by the a-priori knowledge about the image (information), the data and their relation to each other. The source of this knowledge is either our understanding of the type of images we want to reconstruct and of the physics behind the acquisition of the data or we can thrive to learn parametric models from the data itself. The common question arises: how can we customise our model choice to a particular application? Or better how can we make our model adaptive to the given data?

Starting from the first modelling strategy this talk will lead us from nonlinear diffusion equations and subdifferential inclusions of total variation type functionals as the most successful image modeltoday to non-smooth second- and third-order variational models, with data models for Gaussian and Poisson distributed data as well as impulse noise. These models exhibit solution-dependent adaptivities in form of nonlinearities or non-smooth terms in the PDE or the variational problem, respectively. Applications for image denoising, inpainting and surface reconstruction are given. After a critical discussion of these different image and data models we will turn towards the second modelling strategy and propose to combine it with the first one using a PDE constrained optimisation method that customises a parametrised form of the model by learning from examples. In particular, we will consider optimal parameter derivation for total variation denoising with multiple noise distributions and optimising total generalised variation regularisation for its application in photography.

Equations of quantum mechanics in the semiclassical regime present an enduring challenge for numerical analysts, because their solution is highly oscillatory and evolves on two scales. Standard computational approaches to the semiclassical Schrödinger equation do not allow for long time integration as required, for example, in quantum control of atoms by short laser bursts. This has motivated our approach of asymptotic splittings. Combining techniques from Lie-algebra theory and numerical algebra, we present a new computational paradigm of symmetric Zassenhaus splittings, which lends itself to a very precise discretisation in long time intervals, at very little cost. We will illustrate our talk by examples of quantum phenomena – quantum tunnelling and quantum scattering – and their computation and, time allowing, discuss an extension of this methodology to time-dependent semiclassical systems using Magnus expansions

This talk focuses on the direct search method, arguably one of the simplest optimization algorithms. The algorithm minimizes an objective function by iteratively evaluating it along a number of (polling) directions, which are typically taken from so-called positive spanning sets. It does not use derivatives.

We first introduce the worst case complexity theory of direct search, and discuss how to choose the positive spanning set to minimize the complexity bound. The discussion leads us to a long-standing open
problem in Discrete Geometry. A recent result on this problem enables us to establish the optimal order for the worst case complexity of direct search.

We then show how to achieve even lower complexity bound by using random polling directions. It turns out that polling along two random directions at each iteration is sufficient to guarantee the convergence
of direct search for any dimension, and the resultant algorithm enjoys lower complexity both in theory and in practice.

The last part of the talk is devoted to direct search based on inaccurate function values. We address three questions:
i) what kind of solution can we obtain by direct search if the function values are inaccurate? 
ii) what is the worst case complexity to attain such a solution? iii) given
the inaccuracy in the function values, when to stop the algorithm in order
to guarantee the quality of the solution and also avoid “over-optimization”?

This talk is based on joint works with F. Delbos, M. Dodangeh, S. Gratton, B. Pauwels, C. W. Royer, and L. N. Vicente.