Workshop on fast solvers for fractional diffusion problems

UPDATE: 16 March, 2020

 Please note that, in light of the covid-19 situation, this workshop will not go ahead as a face-to-face meeting. updates will be sent to all registered attendees.

2 April, 2020
Technology and Innovation Centre,
University of Strathclyde, Glasgow, G1 1RD
This workshop will bring together those actively working on solvers for fractional diffusion problems, and researchers whose tools and expertise may help to solve some of the remaining challenges in the field. 

Invited speakers:

The programme and abstracts are below. 


There will also be a conference dinner on the evening of April 2. More details will be available shortly. 

There is no registration fee, but for catering purposes please register using the link below.

For additional information, contact Jennifer Pestana

Spectral Fractional Diffusion via Decoupled Reaction Diffusion Equations
Lehel Banjai
Joint work with J. Markus Melenk and Ch. Schwab
We discuss the discretization of the spectral fractional Laplacian and the efficient solution of the arising linear systems. We consider two representations of the fractional Laplacian - one as the Dirichlet-to-Neumann operator of a degenerate but local elliptic problem posed in an infinite cylinder (sometimes called the Caffarelli-Silvestre extension) and the second the representation by the Dunford-Taylor integral. We show how a combination of an hp-FEM in the extended direction and a diagonalization procedure leads to a decoupled system of reaction diffusion systems for the first approach. A careful choice of parameters in the hp-FEM is needed in order to obtain stable results. A quadrature of the Dunford-Taylor integral leads directly to such a system as well. We describe how this decomposition into reaction diffusion problems aids both the analysis, e.g., of a fully hp-FEM discretization,  and the efficient solution. Finally we compare the two approaches.
Parallel solution of large sparse systems by direct and hybrid methods
Iain Duff

We discuss algorithms and codes for the solution of sparse systems that we have developed in an EU Horizon 2020 Project called NLAFET.

We used two approaches to get good single node performance. For symmetric sys- tems we used task-based algorithms based on an assembly tree representation of the factorization. We then used runtime systems for scheduling the computation on both multicore CPU nodes and GPU nodes [4]. The second approach was to design a new parallel threshold Markowitz algorithm [2]. This second approach was designed for highly unsymmetric systems, and we do not consider this further in this talk.

For the first approach, in addition to discussing the matrix factorization, we consider the efficient solution of systems using the factors from this factorization [1].

We then extended the scope of these two approaches to exploit distributed memory parallelism. In both cases, we partition the matrix and run our single-node code on each partition. We discuss the first case, where we base our work on the block Cimmino algorithm [3] using the ABCD software package coded by Zenadi in Toulouse [6]. We partition the matrix using a numerically aware partitioning algorithm [5]. The kernel for this algorithm is the direct factorization of a symmetric indefinite submatrix for which we use the symmetric tree-based code.

We show the performance of our codes on industrial strength large test problems on a heterogeneous platform. Our codes that are available on github are shown to outperform other state-of-the-art codes.

Rational Krylov for fractional diffusion problems: convergence and pole selection
Stefano Massei

Evaluating the action of a matrix function on a vector, that is x = f(M)v, is an ubiquitous task in applications. When the matrix M is large, subspace projection methods, such as the rational Krylov method, are usually employed. In this work, we provide a quasi-optimal pole choice for rational Krylov methods when f is either Cauchy-Stieltjes or Laplace-Stieltjes for a positive definite matrix M. This includes the case f(z) = z^(−α) with α > 0, which arises when solving fractional diffusion problems.

Then, we consider the case when the argument M has the Kronecker structure M = I ⊗ A + B ⊗ I, and v is obtained by vectorizing a low-rank matrix. This finds application, for instance, in solving two dimensional FDEs on rectangular domains. We introduce an error analysis for the numerical approximation of x. Pole choices and explicit convergence bounds are given also in this case.

Multigrids for both isotropic and anisotropic space-fractional diffusion equations
Mariarosa Mazza
Joint work with Marco Donatelli,  Rolf Krause and Ken Trotti

We are interested in a finite difference discretization of a two-dimensional time-dependent space- fractional diffusion equation (FDE) and in the solution of the related linear systems. Our focus is both on the case where the fractional orders are close to each other (isotropic case), and the one where they are not (anisotropic case). Driven by the fact that the discretization matrices have a block-Toeplitz-Toeplitz-blocks-like structure, and by the well-known negative results on multilevel circulant preconditioning, we opt for a multigrid approach. In this framework, the spectral properties of the matrices and the isotropic/anisotropic nature of the problem guide the selection of both projector and smoother.

In the isotropic case, based on the spectral analysis of the coefficient matrices, we define a multigrid method with classical linear interpolation as grid transfer operator and damped-Jacobi as smoother. In the anisotropic case, inspired by certain multigrid strategies for integer order anisotropic diffusion equations given in literature, we replace the classical linear interpolation with a semi-coarsening technique. Moreover, we estimate the Jacobi relaxation parameter by using an automatic spectral- based procedure. A further improvement in the robustness of the proposed method with respect to the anisotropy of the problem is attained employing the resulting V-cycle with semi-coarsening as smoother inside an outer full-coarsening.

Fast Solutions Methods for Convex Quadratic Fractional Differential Equation Optimization
Spyridon Pougkakiotis

Joint work with John W. Pearson, Santolo Leveque and Jacek Gondzio

We present an optimization method suitable for solving convex quadratic Partial Differential Equation (PDE) optimization problems with box constraints on the state and control variables. In particular, we assume that we are given an arbitrary PDE constrained inverse problem, an associated discretization method, and that the resulting sequences of discrete matrices belong to the class of Generalized Locally Toeplitz (GLT) sequences. Then, we propose the use of an Alternating Direction Method of Multipliers (ADMM) for solving the discretized optimization problems. We show that the linear systems required to be solved during the iterations of ADMM preserve the GLT structure of the initial problem matrices. Using this structure, we present and analyze some general methodologies for preconditioning efficiently such linear systems, and solving them using an appropriate Krylov subspace method.

Subsequently, we focus on a certain class of convex quadratic optimization problems with Fractional Differential Equation (FDE) constraints. Discretized versions of FDEs involve large dense linear systems. In order to overcome this difficulty, we design a recursive linear algebra, based on the fast Fourier trans- form (FFT). Focusing on time-dependent 2-dimensional FDEs, we manage to keep the storage requirements linear, with respect to the grid size N, while en- suring an order N log N computational complexity per iteration of the Krylov solver inside ADMM. We implement the proposed method, and demonstrate its robustness and efficiency, through a series of experiments over different setups of the FDE optimization problem.


Parallel preconditioning for time-dependent PDEs and PDE control
Andy Wathen

Joint work with Elle McDonald, Jennifer Pestana, Anthony Goddard and Federico Danieli

We present a novel approach to the solution of time-dependent PDEs via the so-called monolithic or all-at-once formulation. This approach will be explained for simple parabolic and hyperbolic problems and its utility in the context of PDE constrained optimization problems will be elucidated.

The underlying linear algebra includes circulant matrix approximations of Toeplitz-structured matrices and allows for effective parallel implementation. Supporting theoretical results and simple computational results will be shown for the heat equation and the wave equation which indicate the potential as a parallel-in-time method.


Jennifer Pestana

Department of Mathematics and Statistics

University of Strathclyde

Glasgow. UK

© 2020 by Jennifer Pestana. Supported by EPSRC Grant EP/R009821/1. Proudly created with