## SIAM National Student Chapter Conference 2012 (SNSCC12)

### Plenary Talks

Prof. Nick Trefethen (Oxford) - How Chebfun Solves ODEs and Eigenvalue Problems

Chebfun has become a remarkably powerful tool for solving differential equations and eigenvalue problems in one dimension. This talk will describe some of the algorithms that make this possible, which are mainly due to Asgeir Birkisson, Toby Driscoll, and Nick Hale. In particular we will focus on

• the "happiness iteration" underlying Chebfun adaptivity
• rectangular spectral collocation discretizations
• how EIGS selects its eigenfunctions
• Frechet derivative operators for nonlinear problems via AD
• the graphical user interface CHEBGUI

Edvin Deadman (NAG) - A Recursive Blocked Schur Algorithm for Computing the Matrix Square Root

The Schur algorithm for computing a matrix square root reduces the matrix to the Schur triangular form and then computes a square root of the triangular matrix. In this talk I will describe a recursive blocking technique in which the computation of the square root of the triangular matrix can be made rich in matrix multiplication. The excellent numerical stability of the point algorithm is shown to be preserved by recursive blocking. Numerical experiments using serial code show significant speedups over the point algorithm. However, in multithreaded implementations for the NAG Library, standard blocking is found to offer the greatest speedups.

Prof. Ernesto Estrada (Strathclyde) - An Invitation to Complex Networks

The field of complex network is introduced informally through several examples. The classical concepts of âsmall-worldâ and âscale-freenessâ are briefly discussed. Then, the problem of communicability in complex networks is motivated and analyzed by considering the use of matrix functions. The concept of network communicability is then applied to a few real-world situations. It motivates a new Euclidean metric for networks, which allows to embed every network into a hypersphere of certain radius. Finally we discuss dynamical processes on networks. In particular we show how to extend these concepts to the consideration of long-range interactions among the agents in complex networks. The consequences of these extensions for real-world situations are briefly discussed.

Dr Will Parnell (Manchester) - Elastodynamic Cloaking - How to Make an Invisible Building

Interest in cloaking theory (i.e. rendering objects near-invisible to incident waves) and its practical realization has grown significantly since the early theoretical work in 2006 of Leonhardt and the Pendry group in optics and electromagnetism respectively. Methods have largely been based on the idea of coordinate transformations, which motivate the design of cloaking metamaterials. These materials are able to guide waves around a specific region of space. Research has subsequently focused on the possibility of cloaking in the contexts of acoustics, surface waves in fluids, heat transfer, fluid flow and linear elastodynamics.

Elastic waves are important in a number of application areas, e.g. vibration control, seismic surveying and earthquakes. In particular as with lots of waves, it is frequently of interest to be able to guide and control the propagation of elastic waves.

It was shown by Milton and co-workers that cloaking from elastic waves (elastodynamic cloaking) is made difficult due to the lack of invariance of Navier's equations under general coordinate transformations which retain the symmetries of the elastic modulus tensor. Invariance of the governing equations can be achieved if assumptions are relaxed on the minor symmetries of the elastic modulus tensor but commonly occurring elastic materials do not possess this property.

In this talk we shall describe the notion of an elastic wave cloak', how the transformation method works and describe some of the existing difficulties with the theory and its practical realization. We shall also describe some recent work that has shown that it is theoretically possible to construct elastodynamic cloaks by pre-stressing hyperelastic (e.g. rubbery) elastic solids. We shall discuss an initial simple case (antiplane waves) as was studied in [1] and describe various generalizations including finite cloaks [2] and more general elastodynamic wave cloaking problems.

[1] Parnell, W.J. 2012 "Nonlinear pre-stress for cloaking from antiplane elastic waves", Proc. Roy. Soc. A 468, 563-580
[2] Parnell, W.J., Norris, A.N. and Shearer, T. 2012 "Employing pre-stress to generate finite cloaks for antiplane elastic waves". Appl. Phys. Letters 100, 171907.

### PhD / Postdoc Talks

#### 11.20 - 12.40 Room 1

Martin Takac (Edinburgh) - How to Climb a Billion Dimensional Hill using a Coin and a Compass and Count the Steps Before Departure

With growing digitization of the world it is increasingly easier to collect monstrous amounts of data. Often, this data is analyzed using an optimization algorithm, and leads to difficult huge-scale problems with millions or billions of variables.

Existing optimization algorithms, which are perfectly suited for solving problems of medium size, such as polynomial-time interior-point methods, are often not useful in this new setting due to the bad dependence of their complexity on the problem dimension. Hence, there is a pressing need to devise and analyze new methods, or adapt classical methods, that would be capable of working in the emerging huge-scale setting. An entirely different approach to the scale problem is via acceleration of existing methods on parallel computing architectures such as many-core computers and clusters thereof, or systems based on graphical processing units (GPUs).

In this talk we describe a new method that combines the two approaches outlined above. Our method has both i) a good dependence on the problem dimension and b) is parallel in nature, and hence is well-suited for solving certain structured optimization problems of huge dimension. In particular, we show that randomized block coordinate descent methods, such as those developed in [1, 2], can be accelerated by parallelization when applied to the problem of minimizing the sum of a partially block separable smooth convex function and a simple block separable convex function.

Many problems of current interest in diverse communities (statistics, optimal experimental design, machine learning, mechanical engineering), can be cast in this form, including least-squares, L1 regularized least-squares (LASSO), group and sparse group LASSO , computing c and A-optimal designs of statistical experiments, training (sparse) linear support vector machines (SVM) and truss topology design (TTD).

We describe a generic parallel randomized block coordinate descent algorithm (PR-BCD) and several variants thereof based on the way parallelization is performed. In all cases we prove iteration complexity results, i.e., we give bounds on the number of iterations sufficient to approximately solve the problem with high probability. Our results generalize the intuitive observation that in the separable case the theoretical speedup caused by parallelization should be equal to the number of processors. We show that the speedup increases with the number of processors and with the degree of partial separability of the smooth component of the objective function. Our analysis also works in the mode when the number of blocks being updated at each iteration is random, which allows for modelling situations with variable (busy or unreliable) number of processors.

We conclude with some encouraging computational results applied to huge-scale LASSO , SVM and TTD instances.

All results are based on joint work with Peter Richtarik (Edinburgh)

[1] Yurii Nesterov. Efficiency of coordinate descent methods on huge-scale optimization problems. CORE Discussion Paper #2010/2, February 2010.
[2] Peter Richtarik and Martin Takac. Iteration complexity of randomized block-coordinate descent methods for minimizing a composite function. After 1st revision in Mathematical Programming A, 2011.

Waqquas Bukhsh (Edinburgh) - Local Solutions of Optimal Power Flow

Optimal power flow (OPF) is a well studied optimization problem in electricity industry. Over the last two decades it has become a standard tool for planning, real time operations and market auctions.

OPF is nonlinear optimization problem and the existence of locally optimal solutions has been a question of interest for decades. Often it is conjectured that OPF feasible region is convex. In this talk, we present examples of local solutions of OPF on a range of power systems networks. We also show that a recent reformulation of OPF as SDP problem sometimes fails to recover feasible solutions of OPF.

Andreas Papayiannis (Manchester) - Continuous-Time Revenue Management in Carparks

We study optimal revenue management applied to carparks, with primary objective to maximize revenues under a continuous-time framework. We develop a stochastic discrete-time model and propose a rejection algorithm that makes optimal decisions (accept/reject) according to the future expected revenues generated and on the opportunity cost that arises before each sale. For this aspect of the problem, a Monte Carlo approach is used to derive optimal rejection policies. We then extend this approach to show that there exists an equivalent continuous-time methodology that yields to a partial differential equation (PDE). The nature of the PDE, as opposed to the Monte Carlo approach, generates the rejection policies quicker and causes the optimal surfaces to be significantly smoother. However, because the solution to the PDE is considered not to solve the 'full' problem, we propose an approach to generate optimal revenues using the discrete-time model by exploiting the information coming from the PDE. We give a worked example of how to generate near-optimal revenues with an order of magnitude decrease in computation speed.

Mingliang Cheng (Manchester) - Real Options with Cash Holding - The Forgotten Dimension

This presentation provides an overview of how cash holdings can be used within Real Options Analysis, to solve long-standing corporate finance problems. We argue that the inclusion of (path dependent) cash reserves should have been used since the beginning of Real Options theory, which began back in the late seventies (because the risk structure of a firm is inherently linked to the operational flexibility a cash reserve provides). However, most of the literature has (often intentionally) ignored the dimension of cash holdings. This is not because it was viewed as unimportant, but rather because of the additional complexity it brings to the underlying mathematics. We show how this complexity may be surmounted, via suitable modelling and numerical techniques, and how this allows many unsolved problems to finally render themselves tractable.

#### 11.20 - 12.40 Room 2

Nicholas Bird (Reading) - A Velocity-Based Moving-Mesh Finite Element Method

We present a velocity-based moving-mesh finite element method for use in approximating solutions to the fourth order nonlinear PDE $u_t = -(uu_{xxx})_x$. The method uses a local conservation of mass principle and maintains scale invariance.

By using a combination of piecewise quartic, quadratic, and linear basis functions we are able to show that when the initial conditions are taken from a self-similar solution the approximation matches the self-similar solution to within rounding error for all time.

Damon McDougall (Warwick) - Uncertainty Quantification in Lagrangian Data Assimilation

It is of great interest within the scientific community to understand uncertainty. Moreover, understanding uncertainty regarding prediction in weather, climate and ocean models is of utmost importance for the general public. In this talk we will set up the Lagrangian data assimilation problem in a Bayesian framework for a testbed model for ocean gliders in a kinematic travelling wave. We will then expose the effect on the variance of the posterior distribution as we force ocean gliders to cross oceanic transport boundaries and into different flow regimes.

Tommaso Biancalani (Manchester) - Stochastic Patterns in Population Models

The emergence of spatial-temporal patterns in populations is (often) mathematically understood via a linear instability of an homogeneous state. This mechanism is sometimes called "Turing instability" and it has been developed in the framework of partial differential equations (PDEs). However, using PDEs requires that every population is considered as a continuous fluid, a non-appropriate description when the number of degrees of freedom is not macroscopic. To overcome this, stochastic models have been introduced, appropriately given the noise created as a consequence of the discreteness of the system. We show that in stochastic models the parameter region for which patterns arise is highly enlarged. These âstochastic patternsâ appear to be robust and they do not need the requirement of the large separation of diffusivities, like their deterministic counterpart. To study them, we present a simple 2-species model which displays both steady patterns and travelling waves. We show numerically the existence of the stochastic patterns and we analysed them analytically.

[1] T. Biancalani, F. Di Patti and D. Fanelli âStochastic Turing patterns in a Brusselator modelâ, Phys. Rev. E 81, 046215 (2010)
[2] T. Biancalani, T. Galla and A. J. McKane âStochastic waves in a Brusselator model with nonlocal interactionâ, Phys. Rev. E 84, 026201 (2011) Abstract to appear soon.

Robin Thompson (Cambridge) - Modelling the Biological Control of Crop Disease

Global food security concerns are likely to remain prominent as population growth and climate change increase pressures on food production. With commonly-used techniques for controlling crop diseases, such as fungicides and crop rotation, subject to ecological, economic and environmental disadvantages, experimentalists have become increasingly interested in biological control, where a natural enemy of a pathogen is deployed to reduce disease.

We initially consider the probability of invasion in the stochastic SIS epidemiological model, focusing on the case where the population size is relatively small. We demonstrate analytically how the probability of an epidemic occurring depends sensitively on the exact definition used to define an epidemic, and how the response is affected by both the population size and the basic reproductive number of the pathogen. Our analysis is then extended to deterministic and stochastic adaptations of the SIS model where an agent of biological control is incorporated. In particular, we determine conditions for biological control to be able to eradicate a particular pathogen, and how eradication is conditioned on the epidemiology of the underlying host-pathogen interaction as well as the level and timing of application. Our results have clear significance in terms of optimising the deployment of biological control and explaining the variability that has all too often beset performance in the field.

#### 14.00 - 15.20 Room 1

Dmitry Savostyanov (Moscow) - Superfast Inversion of Triangular Toeplitz Matrices in QTT Format with Applications to Fractal Calculus

The Riemann-Lioville integral operator, which plays the central role in fractional calculus, is defined as follows
$D^{-\alpha}f(t)=\frac{1}{\Gamma(\alpha)}\int\limits_0^t(t-\tau)^{\alpha-1}f(\tau)d\tau.$

Due to the properties of the core function, the discretization of the operator gives triangular $n\timesn$ matrix, which can be inverted in $\mathcal{O}(n^2)$ operations. If uniforms grid is used, the matrix $\mathcal{A}=[a(j,k)]$ is also Toeplitz, i.e., $a(j,k)=a(j-k)$, which allows $\mathcal{O}(n{\log}n)$ inversion algorithms based on Fast Fourier transform.

We approach the inversion algorithms for triangular Toeplitz matrices using the recently proposed QTT format [1,2], i.e.,
$a(k)=a(\overbar{k_1\dotsk_d})=\mathcal{A}_{k_1}^{(1)}\dots\mathcal{A}_{k_d}^{(d)},$
where $k=0,\dots,n-1,\quad{n=2^d},\quad{k=\overbar{k_1\dotsk_d}}=\sum_{b=1}^dk_b2^{b-1}$, and each $\mathcal{A}_{k_b}^{(b)}$ is $r_{b-1}{\times}r_{b-1}$ matrix, with $r_0=r_d=1$. This reduces the storage to $\mathcal{O}(dr^2)$ and complexity to $\mathcal{O}(dr^4)$, where $d=\log_2n$ and $r=\max{r_i}$. Therefore, for a class of vectors with $r\le{\log}n$, we result in algorithms of sublinear complexity, which allow the use of very fine grids.

By numerical examples we demonstrate the advances of QTT-based algorithms for fractional calculus problems and compare them with algorithms based on FTT as well as 'finite memory' approach.

[1] I. V. Oseledets. Tensor-train decomposition. SIAM J. Sci. Comput., 33(5):2295-2317, 2011
[2] B. N. Khoromskij. $\mathcal{O}(d{\log}n)$-Quantics approximation of N-d tensors in high-dimensional numerical modeling. Constr. Appr., 34(2):257-280, 2011.

Shengxin Zhu (Oxford) - Which is the Most Important Kernel for Sparse Linear Solves on Heterogeneous Supercomputers?

Recent years have witnessed that traditional iterative solvers are not suitable for distributed super-computers due to intensively global communications. Unlike on share memory computers in which the sparse matrix vector multiplication is usually the most time-consuming part, inner product computation on most of current heterogeneous supercomputer could dominate the entire simulation time and become the bottle-neck due to global communications. Theoretical analysis and experiments will let us rethink which is the most important kernel for next generation sparse linear solvers: matrix-vector multiplication or inner product computation? Numerical experiments on a Dawning 5000A, one of the most fast heterogeneous supercomputers, will be reported.

Nick Dingle (Manchester) - Reducing the Influence of Tiny Normwise Relative Errors on Performance Profiles

It is a widespread but little-noticed phenomenon that the normwise relative error of vectors x and y of floating point numbers, where y is an approximation to x, can be many orders of magnitude smaller than the unit roundoff. In this talk I will explain why such a seemingly counter-intuitive result can be observed and highlight why it causes problems when we use performance profiles to compare the accuracy of our algorithms. I will then present a remedy that reduces the influence of these extreme errors in a controlled manner while preserving the monotonicity of the underlying data.

Bubacarr Bah (Edinburgh) - Compressed Sensing with Sparse Random Matrices from Expander Graphs

Using a probabilistic construction and a novel technique of set-splitting we achieved a great improvement in the parameters of a lossless expander graph and consequently an improvement on the Compressed Sensing recovery guarantees of sparse non-mean zero binary matrices which are adjacency matrices of lossless expander graphs.

#### 14.00 - 15.20 Room 2

Simon Gaulter (Reading) - Acoustic Trapped Modes of a Slowly Varying 3D Waveguide

A trapped mode is an eigensolution to the time-harmonic wave equation, characterised by a localised region of finite oscillatory energy and exponential decay elsewhere. We present an asymptotic scheme to approximate the trapped mode solutions in a longitudinally symmetric three-dimensional waveguide with a smooth but otherwise arbitrarily shaped cross-section and a single bulge varying over a large length scale, $1/\epsilon$ for $0<\epsilon{\ll}1$.

A set of so-called quasi-mode' solutions are constructed via a modified WKBJ-type expansion around the small parameter, $\epsilon$, and used to find the cut-off and cut-on common of the waveguide. A uniform solution, valid in both these `outer' regions and at the turning point, is proposed and utilising the symmetry of the waveguide about $x=0$, we use this to reduce the calculation of the leading order trapped mode wavenumbers to the solution of a scalar transcendental equation. These approximate solutions are compared with those obtained via a three-dimensional pseudospectral numerical scheme and seen to be extremely accurate.

Christopher Rowlatt (Cardiff) - On the Immersed Boundary (IB) Method and its Applications to a Viscoelastic Fluid

In a classical fluid-structure interaction problem, the fluid and the structure are considered separately then coupled together via some suitable conditions. In the immersed boundary (IB) method [C.S. Peskin, J. Comp. Phys. 10:252-271, 1972] the structure is considered as part of the fluid and additional mass/forces may be imposed in the region where the structure resides. In this talk we apply a spectral element approach to the finite element immersed boundary (FE-IB) method of Boffi et al. [D. Boffi, L. Gastaldi, L. Heltai and C.S. Peskin, Comput. Meth. Appl. Mech. Engrg. 197:2210-2231, 2008]. We present results for a closed elastic membrane immersed in a Newtonian fluid before applying the method to a viscoelastic fluid. A proposed extension to allow for discontinuous viscosity will be briefly discussed at the end of the talk, if time allows.

Maria Veretennikova (Warwick) - Control Fractional Dynamics

Firstly, you will be introduced to fractional calculus which has recently gained popularity in a wide range of fields, in particular establishing itself useful in modeling anomalous diffusion by suitable continuous time random walks (CTRWs). Secondly, you will see how to write a dynamic programming equation for the payoff function for a process in our consideration which is derived from a controlled CTRW, and how scaling affects it. You will see the new equations derived in my research for the different versions on the process. We will then discuss resulting fractional Hamilton Jacobi Bellman type equations for the payoff functions.

Andrew Lam (Warwick) - Diffuse Interface Modelling of Surfactants in Two-Phase Flow

We present the derivation of a diffuse interface model of surfactants in a mixture of two Newtonain fluids with different densities. It can be shown that the model satisifies an energy inequality and it approximates a free boundary model in the sharp interface limit. We then present 1D numerics to compare the diffuse interface model and the Ward-Tordai problem.

### Posters

Samuel Groth - Frequency Independent Models for Scattering Atmospheric Particles

Martin Takac - A Randomized Coordinate Descent Method For Large-Scale Truss Topology Design

Waqquas Bukhsh - An MILP Model for Preventative Islanding of Power Systems

Ndifreke Udosen - Automated Optimisation of Measurement Locations to solve Meta Inverse Problems

Pablo Gonzalez-Brevis - A Natural Stabilized Column Generation Method

Yiming Yan - Optimal Active-Set Prediction for Interior Point Methods

Nick Dingle - GPU-Accelerated Solution of Systems of Sparse Linear Equations

Haojie Lin - Sensitivity Analysis and Shape Optimisation with an Isogeometric Boundary Element Method

Ramaseshan Kannan - Sparse Matrix Vector Multiplication with a Novel Block Format