University of Manchester

Although Toeplitz matrices are often dense, matrix-vector products with Toeplitz matrices can be quickly performed via circulant embedding and the fast Fourier transform. This makes their solution by preconditioned Krylov subspace methods attractive. 

For a wide class of symmetric Toeplitz matrices, symmetric positive definite circulant preconditioners that cluster eigenvalues have been proposed. MINRES or the conjugate gradient method can be applied to these problems and descriptive convergence theory based on eigenvalues guarantees fast convergence. 

In contrast, although circulant preconditioners have been proposed for nonsymmetric Toeplitz systems, guarantees of fast convergence are generally only available for CG for the normal equations (CGNE). This is somewhat unsatisfactory because CGNE has certain drawbacks, including slow convergence and a larger condition number. In this talk we discuss a simple alternative symmetrization of nonsymmetric Toeplitz matrices, that allows us to use MINRES to solve the resulting linear system. We show how existing circulant preconditioners for nonsymmetric Toeplitz matrices can be straightforwardly adapted to this situation and give convergence estimates similar to those in the symmetric case.

Hazardous geophysical mass flows, such as snow avalanches, debris-flows and pyroclastic flows, often spontaneously develop large particle rich levees that channelize the flow and enhance their run-out. Measurements of the surface velocity near an advancing flow front have been made at the United States Geological Survey (USGS) debris-flow flume, where 10m^3 of water saturated sand and gravel are allowed to flow down an 80m chute onto a run-out pad. In the run-out phase the flow front is approximately invariant in shape and advances at almost constant speed. By tracking the motion of surface tracers and using a simple kinematic model, it was possible to infer bulk motion as incoming material is sheared towards the front, over-run and shouldered to the side. At the heart of the levee formation process is a subtle segregation-mobility feedback effect. Simple models for particle segregation and the depth-averaged motion of granular avalanches are described and one of the first attempts is made to couple these two types of models together. This process proves to be non-trivial, yielding considerable complexity as well as pathologies that require additional physics to be included.

We present a block LU factorization with panel rank revealing

pivoting (block LU_PRRP), an algorithm based on strong

rank revealing QR for the panel factorization.

Block LU_PRRP is more stable than Gaussian elimination with partial

pivoting (GEPP), with a theoretical upper bound of the growth factor

of $(1+ \tau b)^{(n/ b)-1}$, where $b$ is the size of the panel used

during the block factorization, $\tau$ is a parameter of the strong

rank revealing QR factorization, and $n$ is the number of columns of

the matrix. For example, if the size of the panel is $b = 64$, and

$\tau = 2$, then $(1+2b)^{(n/b)-1} = (1.079)^{n-64} \ll 2^{n-1}$, where

$2^{n-1}$ is the upper bound of the growth factor of GEPP. Our

extensive numerical experiments show that the new factorization scheme

is as numerically stable as GEPP in practice, but it is more resistant

to some pathological cases where GEPP fails. We note that the block LU_PRRP

factorization does only $O(n^2 b)$ additional floating point operations

compared to GEPP.

Convex regularization has become a popular approach to solve large scale inverse or data separation problems. A prominent example is the problem of identifying a sparse signal from linear samples my minimizing the l_1 norm under linear constraints. Recent empirical research indicates that many convex regularization problems on random data exhibit a phase transition phenomenon: the probability of successfully recovering a signal changes abruptly from zero to one as the number of constraints increases past a certain threshold. We present a rigorous analysis that explains why phase transitions are ubiquitous in convex optimization. It also describes tools for making reliable predictions about the quantitative aspects of the transition, including the location and the width of the transition region. These techniques apply to regularized linear inverse problems, to demixing problems, and to cone programs with random affine constraints. These applications depend on a new summary parameter, the statistical dimension of cones, that canonically extends the dimension of a linear subspace to the class of convex cones.

Joint work with Dennis Amelunxen, Mike McCoy and Joel Tropp.

The use of formal mathematical models in sociology started in the 1940s and attracted mathematicians such as Frank Harary in the 1950s. The idea is to take the rather intuitive ideas described in social theory and express these in formal mathematical terms. Social network analysis is probably the best known of these and it is the area which has caught the imagination of a wider audience and has been the subject of a number of popular books. We shall give a brief over view of the field of social networks and will then look at three examples which have thrown up problems of interest to the mathematical community. We first look at positional analysis techniques and give a formulation that tries to capture the notion of social role by using graph coloration. We look at algebraic structures, properties, characterizations, algorithms and applications including food webs. Our second and related example looks at core-periphery structures in social networks. Our final example relates to what the network community refer to as two-mode data and a general approach to analyzing networks of this form. In all cases we shall look at the mathematics involved and discuss some open problems and areas of research that could benefit from new approaches and insights.

Computing the eigenvalue decomposition of a symmetric matrix and the singular value decomposition of a general matrix are two of the central tasks in numerical linear algebra. There has been much recent work in the development of linear algebra algorithms that minimize communication cost. However, the reduction in communication cost sometimes comes at the expense of significantly more arithmetic and potential instability.

\\

\\

In this talk I will describe algorithms for the two decompositions that have optimal communication cost and arithmetic cost within a small factor of those for the best known algorithms. The key idea is to use the best rational approximation of the sign function, which lets the algorithm converge in just two steps. The algorithms are backward stable and easily parallelizable. Preliminary numerical experiments demonstrate their efficiency.

The Cholesky factorization approach to solving the symmetric definite generalized eigenvalue problem

$Ax = \lambda Bx$, where $A$ is symmetric and $B$ is symmetric positive definite, computes a Cholesky factorization $B = LL^T$ and solves the equivalent standard symmetric eigenvalue problem $C y = \l y$ where $C = L^{-1} A L^{-T}$. Provided that a stable eigensolver is used, standard error analysis says that the computed eigenvalues are exact for $A+\dA$ and $B+\dB$ with $\max( \normt{\dA}/\normt{A}, \normt{\dB}/\normt{B} )$

bounded by a multiple of $\kappa_2(B)u$, where $u$ is the unit roundoff. We take the Jacobi method as the eigensolver and explain how backward error bounds potentially much smaller than $\kappa_2(B)u$ can be obtained.

To show the practical utility of our bounds we describe a vibration problem from structural engineering in which $B$ is ill conditioned yet the error bounds are small. We show how, in cases of instability, iterative refinement based on Newton's method can be used to produce eigenpairs with small backward errors.

Our analysis and experiments also give insight into the popular Cholesky--QR method used in LAPACK, in which the QR method is used as the eigensolver.

In this presentation we examine the convergence characteristics of a

Krylov subspace solver preconditioned by a new indefinite

constraint-type preconditioner, when applied to discrete systems

arising from low-order mixed finite element approximation of the

classical biharmonic problem. The preconditioning operator leads to

preconditioned systems having an eigenvalue distribution consisting of

a tightly clustered set together with a small number of outliers. We

compare the convergence characteristics of a new approach with the

convergence characteristics of a standard block-diagonal Schur

complement preconditioner that has proved to be extremely effective in

the context of mixed approximation methods.

\\

\\

In the second part of the presentation we are concerned with the

efficient parallel implementation of proposed algorithm on modern

shared memory architectures. We consider use of the efficient parallel

"black-box'' solvers for the Dirichlet Laplacian problems based on

sparse Cholesky factorisation and multigrid, and for this purpose we

use publicly available codes from the HSL library and MGNet collection.

We compare the performance of our algorithm with sparse direct solvers

from the HSL library and discuss some implementation related issues.