Conference on Fast Direct Solvers (Online)

Department of Mathematics, Purdue University, October 23-24, 2021

Program

To receive the online meeting information, please register.

 

(Time zone: US Eastern Time. Please click on the talk titles for the abstracts)

Saturday, October 23
Session 1 (Chair: Jianlin Xia)
8:20-8:30am Welcome Remark
8:30-9:20

George Biros, University of Texas at Austin

An Overlapping Domain Decomposition Preconditioner for Integral Equations

9:20-9:45

Robert Krasny, University of Michigan

Fast Summation Methods based on Barycentric Lagrange Interpolation

9:45-10:10

Yuanzhe Xi, Emory University

Data-Driven Linear Complexity Low-Rank Approximation of General Kernel Matrices

10:10-10:35

Alice Cortinovis, Ecole Polytechnique Federale de Lausanne

Divide and Conquer Methods for Functions of Matrices with Banded or Hierarchical Low-Rank Structure

10:35-10:50 Break
Session 2 (Chair: Bryan Quaife)
10:50-11:40

Adrianna Gillman, University of Colorado, Boulder

Solvers for Applications of Integral Equations

11:40-12:05

Stefano Massei, Eindhoven University of Technology

Hierarchical Adaptive Low-Rank Format with Applications to Discretized PDEs

12:05-12:30

Daria Sushnikova, NYU

Direct Solution of Systems with Rank-Structured Matrices

12:30-1:40 Lunch Break
Session 3 (Chair: Robert Krasny)
1:40-2:30pm

Michael W. Mahoney, UC Berkeley

Recent Developments in the Theory and Practice of Second Order Optimization for Machine Learning

2:30-2:55

Alex Townsend, Cornell University

Learning Elliptic Partial Differential Equations with Randomized Linear Algebra

2:55-3:20

Vivak Patel, University of Wisconsin -- Madison

Incrementalizing Random Sketching for Solving Consistent Linear Systems

3:20-3:45

Nick Polydorides, University of Edinburgh, UK

Low Variance Sketched Finite Elements for Elliptic Equations

3:45-4:00 Break
Session 4 (Chair: Yuanzhe Xi)
4:00-4:50

Xiaoye S. Li, Lawrence Berkeley Lab

Scalable and Robust Hierarchical Matrix Factorizations via Randomization

4:50-5:15

Chao Chen, University of Texas at Austin

Distributed-Memory Parallel Algorithms Based on Rank-Structured Matrices

5:15-5:40

Heather Wilber, University of Texas at Austin

Designing Low Rank Methods for Matrices with Displacement Structure

5:40-6:05

Xin Xing, University of California, Berkeley

Efficient Construction of An HSS Preconditioner for Symmetric Positive Definite H2 Matrices

6:05-6:30pm

Mikhail Lepilov, Purdue University

Proxy Point Method for Low-Rank Approximations of Analytic Kernels with Applications to Some Toeplitz Problems

 
Sunday, October 24
Session 5 (Chair: Jie Shen)
8:30-9:20am

Michael Ng, University of Hong Kong

Numerical Methods for Fractional Diffusion Equations

9:20-9:45

Shan Zhao, University of Alabama

A Fourth Order Finite Difference Method for Solving Elliptic Interface Problems with the FFT Acceleration

9:45-10:10

Daniel Fortunato, Flatiron Institute

A Fully Adaptive Poisson Solver for Smooth Two-Dimensional Domains

10:10-10:35

Yong Zhang, Tianjin University

A Spectrally Accurate Numerical Method for Computing the Bogoliubov-de Gennes Excitations of Dipolar Bose-Einstein Condensates

10:35-10:50 Break
Session 6 (Chair: Shan Zhao)
10:50-11:40

Jiang Yang, Southern University of Science and Technology

A Second Order Globally Convergent Method for Computing Top Eigenpair of Irreducible Non-negative Matrices and its Applications

11:40-12:05

Manki Cho, University of Houston - Clear Lake

Applications of the Sloshing Problem to Elliptic Problems via Steklov Eigenpairs

12:05-12:30

Carmen Scalone, University of L'Aquila

Approximating Low-Rank Rightmost Eigenpairs of a Class of Matrix-Valued Linear Operators

12:30-1:40 Lunch Break
Session 7 (Chair: Alex Townsend)
1:40-2:30pm

Qiang Ye, University of Kentucky

Structured Preconditioning for Neural Network Training

2:30-3:20

Lexing Ying, Stanford University

Solving Inverse Problems with Deep Learning

3:20-3:35 Break
Session 8 (Chair: Qiang Ye)
3:35-4:25 Olaf Schenk, Universita della Svizzera italiana, Switzerland
(Talk will be delivered by Aryan Eftekhari, Universita della Svizzera italiana)
Large-Scale Sparse Precision Matrix Estimation
4:25-4:50

Thijs Steel, KU Leuven

A Fast Algorithm for Hessenberg-Triangular Pencil Factorization

(Cancelled) Pablo Brubeck, University of Oxford

Sparse Vertex-Star Relaxation for High-Order FEM

4:50-5:15

Jun Liu, Southern Illinois University Edwardsville

A Well-Conditioned Direct Parallel-in-Time (PinT) Algorithm for Evolutionary Differential Equations

5:15-5:40

Tianyi Shi, Cornell University

Low Tensor-Train Rank Methods to Solve Sylvester Tensor Equations

5:40-6:05

Xiaofeng Ou, Purdue University

SuperDC: Stable Superfast Divide-and-Conquer Eigenvalue Decomposition 

 

Abstracts


 

An Overlapping Domain Decomposition Preconditioner for Integral Equations

George Biros, University of Texas at Austin

The discretization of integral equations leads to systems that dense and often ill-conditioned. We introduce a new preconditioner based on a novel overlapping domain decomposition that can be combined efficiently with fast direct solvers. Empirically, we observe that the condition number of the preconditioned system is O(1), independent of the problem size. Our domain decomposition is designed so that we can construct approximate factorizations of the subproblems efficiently. In particular, we apply the recursive skeletonization algorithm to subproblems associated with every subdomain. We present numerical results on problem sizes with 288 million unknowns in 2D and 17 million unknowns in 3D, which were factorized in less than 4 hours and one hour, respectively, on a single Intel Xeon Platinum 8280M socket.


 

Sparse Vertex-Star Relaxation for High-Order FEM

Pablo Brubeck, University of Oxford

Pavarino proved that the additive Schwarz method with vertex patches and a low-order coarse space gives a $p$-robust solver for symmetric and coercive problems. However, for very high polynomial degree it is not feasible to assemble or factorize the matrices for each patch. In this work we introduce a direct solver for separable patch problems that scales to very high polynomial degree on tensor product cells. The solver constructs a tensor product basis that diagonalizes the blocks in the stiffness matrix for the internal degrees of freedom of each individual cell. As a result, the non-zero structure of the cell matrices is that of the graph connecting internal degrees of freedom to their projection onto the facets. In the new basis, the patch problem is as sparse as a low-order finite difference discretization, while having a sparser Cholesky factorization. We can thus afford to assemble and factorize the matrices for the vertex-patch problems, even for very high polynomial degree. In the non-separable case, the method can be applied as a preconditioner by approximating the problem with a separable surrogate.


 

Distributed-Memory Parallel Algorithms Based on Rank-Structured Matrices

Chao Chen, University of Texas at Austin

Consider solving large-scale linear systems that do not fit into the memory of a single compute node, and thus distributed-memory parallel algorithms are required. However, traditional direct solvers and iterative solvers may suffer from significant communication costs. A lot of recent research has shown that sequential fast direct solvers employing rank-structured matrices outperform traditional methods with much fewer flops. In this talk, we focus on the parallel performance of such fast direct solvers when applied to problems stored on distributed-memory machines. In particular, we show scalability results on up to a thousand processors for solving both sparse linear systems arising from the finite element discretization of elliptic partial differential equations and dense linear systems arising from the discretization of integral equations.


 

Applications of the Sloshing Problem to Elliptic Problems via Steklov Eigenpairs

Manki Cho, University of Houston - Clear Lake

There have been studies of the sloshing problem, a linear eigenvalue problem which could describe the oscillations of the free surface of fluid on a container. Such eigenpairs can be obtained either numerically or explicitly depending on the shape of domain. This talk will feature the use of Steklov eigenfunctions to construct bases of the space of solutions of linear elliptic PDEs. The corresponding spectral series converges rapidly from the boundary to interior of the region. Also this method can be a meshless scheme once explicit formulas of Steklov spectrum are found. Numerical experiments performed in solving Laplace equation subject to mixed boundary conditions arising from electrostatic problems will be presented.


 

Divide and Conquer Methods for Functions of Matrices with Banded or Hierarchical Low-Rank Structure

Alice Cortinovis, Ecole Polytechnique Federale de Lausanne

This talk is concerned with approximating matrix functions for banded matrices, hierarchically semiseparable matrices, and related structures. We propose new divide-and-conquer methods which exploit the fact that these matrices can be (recursively) decomposed as a sum A = D + R of a block diagonal matrix D and a low-rank correction R. While the update f(A) - f(D) often has low numerical rank and can be approximated via (rational) Krylov subspace projections, the block diagonal part f(D) is computed recursively for each diagonal block. We present a convergence analysis that relates the accuracy attained by the algorithm with best polynomial or rational approximations of the function. For the special case of a banded matrix, we show that the divide-and-conquer method reduces to a much simpler algorithm, which proceeds by computing matrix functions of small submatrices of A. When only the trace or the diagonal of the matrix function is of interest, we demonstrate - in practice and in theory - that convergence can be faster. Finally, we test the algorithms on a variety of matrices and functions; the numerical results demonstrate that, most of the time, the proposed methods outperform state-of-art techniques with respect to time consumption and offer a comparable accuracy. This is joint work with Daniel Kressner and Stefano Massei.


 

A Fully Adaptive Poisson Solver for Smooth Two-Dimensional Domains

Daniel Fortunato, Flatiron Institute

We present a new framework for the fast solution of inhomogeneous elliptic boundary value problems in domains with smooth boundaries in 2D. High-order solvers based on adaptive box codes or the FFT can efficiently treat the volumetric inhomogeneity but require care to be taken near the boundary to ensure that the volume data is globally smooth. We avoid function extension and cut-cell quadratures near the boundary by employing a smooth partition of unity defined by a variable-width annular region inside the domain that conforms to the boundary. We use a spectral element method in this region with an accompanying fast direct solver, and patch together the regions using layer potentials. The resulting solver allows for the adaptive resolution of volumetric data, boundary data, and geometric features across a range of length scales. This is joint work with Alex Barnett and David Stein.


 

Solvers for Applications of Integral Equations

Adrianna Gillman, University of Colorado, Boulder

Fast direct solvers are natural to apply to the linear systems that result from the discretization of integral equations. Due to the large upfront cost for constructing the fast solver, the application of these solution techniques has been limited to problems involving many sets of right hand sides and having to pay the up front cost of creating a new solver for each new problem even if it involves a minor change. In this talk, we will present recently developed techniques that allow for a single precomputed fast direct solver to be used for a range of problems related to an initial boundary value problem. The efficiency of these techniques make integral equations a viable solution technique for more applications. Numerical experiments will illustrate the potential impact of the presented fast solvers.


 

Fast Summation Methods based on Barycentric Lagrange Interpolation

Robert Krasny, University of Michigan

We present two fast summation methods based on barycentric Lagrange interpolation for particle interactions. The first is a treecode (BLTC) and the second is a fast multipole method in which the interaction lists are formed by dual tree traversal (BLDTT). The methods are kernel-independent and a distributed memory implementation running on multiple GPUs has been developed. The performance of the BLTC and BLDTT is demonstrated for several particle systems. This is joint work with Lei Wang, Svetlana Tlupova, Leighton Wilson, and Nathan Vaughn.


 

Hyperfast Rank-Structured Approximation of Some Toeplitz Matrices

Mikhail Lepilov, Purdue University

Toeplitz matrices with values taken from a kernel function arise frequently in applications such as image processing. Often, such matrices have low off-diagonal numerical ranks, allowing us to take advantage of their rank-structured approximations to obtain a substantial speedup in various tasks such as factoring, performing matrix-vector products, etc. The construction of such rank-structured approximations is often the bottleneck in this process. In this talk, we propose a sublinear-cost ("hyperfast") construction scheme for a hierarchically semiseparable (HSS) approximation of certain Toeplitz matrices which arise from kernel maps satisfying a certain conformal property. We also provide an analysis of the performance and accuracy of this scheme and propose some extensions of this idea.


 

Scalable and Robust Hierarchical Matrix Factorizations via Randomization

Xiaoye S. Li, Lawrence Berkeley Lab

Hierarchical matrix algebra, stemming from the H-matrix research, is the foundation of the low complexity data-sparse factorization algorithms. In addition, the multilevel algorithm paradigm provides a good vehicle for a scalable implementation on large-scale HPC systems. In this talk, we will discuss how randomized sketching can be used effectively in finding the hierarchical bases in the HSS and Butterfly representations. A key element is a robust adaptive sampling scheme with good stopping criteria. We will illustrate many practical issues of these algorithms using our parallel libraries STRUMPACK and ButterflyPACK, and demonstrate their effectiveness and scalability while solving the very challenging problems, such as high frequency wave equations.


 

A Well-Conditioned Direct Parallel-in-Time (PinT) Algorithm for Evolutionary Differential Equations

Jun Liu, Southern Illinois University Edwardsville

In this talk, we introduce a direct parallel-in-time (PinT) algorithm for evolutionary differential equations with first- or second-order derivative in time. We use a second-order boundary value method as the time integrator and explicitly diagonalize the time discretization matrix, which yields a direct parallel implementation across all $n$ time levels. A crucial issue on this methodology is how the condition number of the eigenvector matrix $V$ of the time discretization matrix behaves as $n$ grows. Based on a novel connection between the characteristic equation and the Chebyshev polynomials, we present explicit formulas for $V$ and $V^{-1}$, by which we proved that $Cond2(V) = O(n^2)$. This implies that the diagonalization process is well-conditioned and the round-off error only increases moderately as $n$ grows. Numerical results on parallel machine are given to support our findings, where over 60 times speedup is achieved with 256 cores.


Recent Developments in the Theory and Practice of Second Order Optimization for Machine Learning

Michael W. Mahoney, ICSI and Department of Statistics, UC Berkeley

Second order optimization methods have been ubiquitous in scientific computing, but they are less used in machine learning, compared to their first order counterparts.  Motivated by well-known problems with first order methods, recent work has begun to experiment with second order methods for machine learning problems.  Here, we review recent developments.  On the theory side, we'll describe how methods from randomized numerical linear algebra (RandNLA) can be used to overcome the so-called inversion bias (basically, if we want to compute a quantity that depends on the inverse of a sum of distributed matrices, then the sum of the inverses does not equal the inverse of the sum), which is a problem in distributed optimization, distributed statistical ensembling, etc.  On the practical side, we'll describe how to use ideas from RandNLA to engineer a ADAHESSIAN, a second order stochastic optimization algorithm which dynamically incorporates the curvature of the loss function via ADAptive estimates of the Hessian, in order to obtain state-of-the-art performance in computer vision, natural language processing, and recommendation systems, with minimal w fiddling.  We'll also discuss future directions.

 

Hierarchical Adaptive Low-Rank Format with Applications to Discretized PDEs

Stefano Massei, Eindhoven University of Technology

A novel compressed matrix format is proposed that combines an adaptive hierarchical partitioning of the matrix with low-rank approximation. One typical application is the approximation of discretized functions on rectangular domains; the flexibility of the format makes it possible to deal with functions that feature singularities in small, localized regions. To deal with time evolution and relocation of singularities, the partitioning can be dynamically adjusted based on features of the underlying data. Our format can be leveraged to efficiently solve linear systems with Kronecker product structure, as they arise from discretized partial differential equations (PDEs). For this purpose, these linear systems are rephrased as linear matrix equations and a recursive solver is derived from low-rank updates of such equations. We demonstrate the effectiveness of our framework for stationary and time-dependent, linear and nonlinear PDEs, including the Burgers’ and Allen-Cahn equations. This is joint work with Daniel Kressner (EPFL) and Leonardo Robol (Università di Pisa).


 

Numerical Methods for Fractional Diffusion Equations

Michael Ng, University of Hong Kong

In this talk, we discuss some preconditioning methods for fractional diffusion equations. Also preconditioning for time fractional diffusion inverse source problems is studied. Numerical examples are reported to demonstrate these preconditioning techniques. Also some deep learning methods are discussed in the talk.


 

SuperDC: Stable Superfast Divide-and-Conquer Eigenvalue Decomposition

Xiaofeng Ou, Purdue University

For dense Hermitian matrices with small off-diagonal (numerical) ranks and in a hierarchically semiseparable form, we give a stable divide-and-conquer eigendecomposition method with nearly linear complexity (called SuperDC) that significantly improves an earlier basic algorithm in [Vogel, Xia, et al., SIAM J. Sci. Comput., 38 (2016)]. We incorporate a sequence of key stability techniques and provide many improvements in the algorithm design. Various stability risks in the original algorithm are analyzed, including potential exponential norm growth, cancellations, loss of accuracy with clustered eigenvalues or intermediate eigenvalues, etc. In the dividing stage, we give a new structured low-rank update strategy with balancing that eliminates the exponential norm growth and also minimizes the ranks of low-rank updates. In the conquering stage with low-rank updated eigenvalue solution, the original algorithm directly uses the regular fast multipole method (FMM) to accelerate function evaluations, which has the risks of cancellation, division by zero, and slow convergence. Here, we design a triangular FMM to avoid cancellation. Furthermore, when there are clustered intermediate eigenvalues or when updates to existing eigenvalues are very small, we design a novel local shifting strategy to integrate FMM accelerations into the solution of shifted secular equations so as to achieve both the efficiency and the reliability. We also provide several improvements or clarifications on some structures and techniques that are missing or unclear in the previous work. The resulting SuperDC eigensolver has significantly better stability while keeping the nearly linear complexity for finding the entire eigenvalue decomposition. In a set of comprehensive tests, SuperDC shows dramatically lower runtime and storage than the Matlab eig function. The stability benefits are also confirmed with both analysis and numerical comparisons.


 

Incrementalizing Random Sketching for Solving Consistent Linear Systems

Vivak Patel, University of Wisconsin -- Madison

Random sketching solvers of linear systems often require allocating a new linear system before applying an efficient direct solver. Unfortunately, the allocation of this intermediate system is a critical bottleneck in the efficiency of such methods. In this talk, we discuss how to incrementally and implicitly construct the randomly sketched system, while simultaneously applying a direct solver to calculate the solution. As a result, we can avoid the costly intermediate step of constructing the sketched system, yet still benefit from the computational advantages offered by sketching.


 

Low Variance Sketched Finite Elements for Elliptic Equations

Nick Polydorides, University of Edinburgh, UK

We present a direct solver for parameter-dependent linear systems of large dimension arising from the application of the finite element method on elliptic boundary value problems. The solver is particularly suited to the many-parameter-query context, typically encountered in uncertainty quantification and inverse problems, when computations must be done on-the-fly and/or without substantial computational resources. Our approach involves low-dimensional subspace projection and randomized sketching of the induced equations, that exploits the positive definite structure of the coefficients matrix and the prior knowledge of this matrix for a uniformly distributed parameter. To suppress the random component of the error in the projected solution we introduce an estimator based on control variates that reduces the variance of the sketched system whilst preserving its well-posedness. Joint work with Robert Lung.


 

Approximating Low-Rank Rightmost Eigenpairs of a Class of Matrix-Valued Linear Operators

Carmen Scalone, University of L'Aquila

We are interested in computing the rightmost eigenpairs of a linear matrix valued operator. We make a priori the hypothesis (as this is a common property of several applications) for which the corresponding eigenmatrix has quickly decaying singular values. This suggests to constrain the search of approximate eigensolutions to a low-rank manifold. We start introducing an appropriate ordinary differential equation, whose solution allows us to approximate the rightmost eigenpair of the linear operator. From the analysis of the behaviour of such ODE on the whole space, we conclude that, under generic assumptions, the solution converges globally to its leading eigenmatrix, when the rightmost eigenvalue is simple and real. After that, we leave the whole space for projecting the differential equation on a low-rank manifold of prescribed rank. The projected operator is nonlinear and this makes the analysis more subtle. Finally, we propose two explicit numerical methods. The numerical experiments show that the approach is effective and competitive. This is a joint work with N. Guglielmi (GSSI) and D. Kressner (EPFL).


 

Large-Scale Sparse Precision Matrix Estimation

Olaf Schenk, Universita della Svizzera italiana

The L1-regularized Gaussian maximum likelihood method is a common approach for sparse precision matrix estimation but one that poses a computational challenge for high-dimensional datasets. I will present a novel method for performant large-scale sparse precision matrix estimation utilizing the block structures and parallelism in the underlying computations for added performance and scalability. Our numerical examples and comparative results with various modern open-source packages reveal that the proposed algorithm provides orders of magnitude faster runtimes and scales to significantly higher-dimensional datasets.


 

Low Tensor-Train Rank Methods to Solve Sylvester Tensor Equations

Tianyi Shi, Cornell University

In this talk, we focus on designing algorithms to solve Sylvester tensor equations in tensor-train (TT) format. The cornerstone of our solvers is the factored alternating direction implicit (fADI) method for Sylvester matrix equations. For different rank structures of the displacement tensor, we combine fADI with TTSVD or parallel-TTSVD to obtain the TT cores of the solution tensor directly. The solvers have optimal running complexity, and are practical in applications such as solving 3D Poisson equations on the cube. This talk is based on joint work with Alex Townsend and Maximilian Ruth.


 

A Fast Algorithm for Hessenberg-Triangular Pencil Factorization

Thijs Steel, KU Leuven

Reducing a pencil to Hessenberg-triangular (HT) form is an important step in the QZ algorithm (which solves dense generalized eigenvalue problems). Due to improvements to the QZ algorithm, the reduction step is now the dominant cost in dense generalized eigenvalue computations. In this talk, we introduce a new reduction algorithm that relies the connection between generalized and standard eigenvalue problems. Numerical experiments show that the algorithm is fast, accurate and easily portable to different architectures.


 

Direct Solution of Systems with Rank-Structured Matrices

Daria Sushnikova, NYU

My primary interest is FMM-based matrices (HSS, H2, etc.), and direct solvers for systems with such matrices. Methods that I am going to present in my talk are based on the idea of the low-rank approximation of the fill-in that appears in the course of the direct solution. The proposed approach allows building direct solvers for the general case of the FMM-based matrix. It has linear scaling and good accuracy for both in 2d and 3d problems on the complex geometries. Solvers successfully applied to numerical solutions of some differential and integral equations and applied statistics problems.

 


 

Learning Elliptic Partial Differential Equations with Randomized Linear Algebra

Alex Townsend, Cornell University

Can one learn a differential operator from pairs of solutions and righthand sides? If so, how many pairs are required? These two questions have received significant research attention in differential equation learning. Given input-output pairs from an unknown elliptic partial differential equation in three dimensions, we will derive a theoretically rigorous scheme for learning the associated Green's function G. By exploiting the hierarchical low-rank structure of Green’s functions and randomized linear algebra, we will describe a rigorous scheme for PDE learning with a provable learning rate. This talk is joint work with Nicolas Boullé.


 

Designing Low Rank Methods for Matrices with Displacement Structure

Heather Wilber, University of Texas at Austin

In this talk, we use rational functions to design low rank methods for computing with matrices that have special displacement structures, including Toeplitz, Cauchy, Vandermonde, and more. The main workhorse of our approach is the alternating direction implicit (ADI) method. Expanding the capabilities of this classic method, we use it to explain the low rank properties of various families of matrices and develop efficient ADI-based algorithms.


 

Data-Driven Linear Complexity Low-Rank Approximation of General Kernel Matrices

Yuanzhe Xi, Emory University

A Geometric Approach Dense kernel matrices arise frequently in many scientific computing and machine learning applications. For large-scale problems, it is common to exploit the block or hierarchically low-rank structures in kernel matrices to accelerate the computation. In this paper, we introduce an efficient data-driven framework for compressing numerically low-rank rectangular kernel matrices. To derive a rank-r approximation efficiently, we propose to first compress the given dataset to O(r) in size and then apply algebraic techniques on the submatrix corresponding to the compressed data. In this framework, matrix compression is driven by data compression (thus termed "data-driven") which can be performed rapidly in linear complexity without accessing the kernel function or matrix. The proposed data-driven low-rank approximation framework applies to general kernel functions and high-dimensional datasets. For datasets in R^d, the proposed two-sided and one-sided factorization schemes can compute a rank-r approximation to an m × n kernel matrix with O(dr(m + n)) and O(dr2(m + n)), respectively. Detailed theoretical analysis has been performed to reveal the impact of compressed data on the low-rank approximation accuracy and suggest the use of data compression schemes that could generate evenly distributed samples. The efficiency and robustness of the proposed framework is demonstrated against a few state-of-the-art methods on various kernel functions and datasets ranging from three to over a hundred dimensions.


 

Efficient Construction of An HSS Preconditioner for Symmetric Positive Definite H2 Matrices

Xin Xing, University of California, Berkeley

In an iterative approach for solving linear systems with dense, ill-conditioned, symmetric positive definite (SPD) kernel matrices, both fast matrix-vector products and fast preconditioning operations are required. Fast (linear-scaling) matrix-vector products are available by expressing the kernel matrix in an H2 representation or an equivalent fast multipole method representation. This talk is concerned with preconditioning such matrices using the hierarchically semiseparable (HSS) matrix representation. Previously, an algorithm was presented to construct an HSS approximation to an SPD kernel matrix that is guaranteed to be SPD. However, this algorithm has quadratic cost and was only designed for recursive binary partitioning of the points defining the kernel matrix. We present a general algorithm for constructing an SPD HSS approximation. Importantly, the algorithm uses the H2 representation of the SPD matrix to reduce its computational complexity from quadratic to quasilinear. Numerical experiments illustrate how this SPD HSS approximation performs as a preconditioner for solving linear systems arising from a range of kernel functions.


 

A Second Order Globally Convergent Method for Computing Top Eigenpair of Irreducible Non-negative Matrices and its Applications

Jiang Yang, Southern University of Science and Technology

In this work, a shifted inverse power method is proposed to compute top eigenpairs of irreducible non-negative matrices. We rigorously prove that this method is second order globally convergent, and the total number of iterations is independent of the matrix size. The effectiveness of the method is due to a convergent upper bound of the dominant eigenvalue. Then, we extend this method to symmetric tridiagonal matrices. Moreover, we propose a novel deflation technique to remove the known eigenpair but still preserve the tridiagonal structure. Thanks to the fact that any symmetric matrix can be tridiagonalized, the proposed shifted inverse power method combining with the novel deflation enables us to compute all eigenpairs of any symmetric matrix efficiently. Several numerical experiments are carried out to demonstrate the efficiency of the proposed method.


 

Structured Preconditioning for Neural Network Training

Qiang Ye, University of Kentucky

Batch normalization (BN) is a popular and ubiquitous method in deep neural network training that has been shown to decrease training time and improve generalization performance. Despite its success, BN is not theoretically well understood. It is not suitable for training with very small mini-batch sizes or online learning. In this talk, we present a structure based preconditioning method called Batch Normalization Preconditioning (BNP) to accelerate neural network training. We will analyze the structure of the Hessian matrix of a loss function, from which the effects of mini-batch statistics of a hidden variable can be deduced. We will propose a parameter transformation that is equivalent to normalizing the hidden variables to improve the conditioning of the Hessian. Compared with BN, one benefit of BNP is that it is not constrained on the mini-batch size and works in the online learning setting. We will present several experiments demonstrating competitiveness of BNP. Furthermore, we will discuss a connection to BN which provides theoretical insights on how BN improves training. The talk is based on a joint work with Susanna Lange and Kyle Helfrich.


 

Solving Inverse Problems with Deep Learning

Lexing Ying, Stanford University

This talk is about some recent progress on solving inverse problems using deep learning. Compared to traditional machine learning problems, inverse problems are often limited by the size of the training data set. We show how to overcome this issue by incorporating mathematical analysis and physics into the design of neural network architectures. We first describe neural network representations of pseudodifferential operators and Fourier integral operators. We then continue to discuss applications including electric impedance tomography, optical tomography, inverse acoustic/EM scattering, seismic imaging, and travel-time tomography.


 

A Spectrally Accurate Numerical Method for Computing the Bogoliubov-de Gennes Excitations of Dipolar Bose-Einstein Condensates

Yong Zhang, Tianjin University

In this paper, we propose an efficient and robust numerical method to study the elementary excitation of dipolar Bose-Einstein condensates (BEC), which is governed by the Bogoliubov-de Gennes equations (BdGEs) with nonlocal dipole-dipole interaction, around the mean field ground state. Analytical properties of the BdGEs are investigated, which could serve as benchmarks for the numerical methods. To evaluate the nonlocal interactions accurately and efficiently, we propose a new Simple Fourier Spectral Convolution method (SFSC). Then, integrating SFSC with the standard Fourier spectral method for spatial discretization and Implicitly Restarted Arnoldi Methods (IRAM) for the eigenvalue problem, we derive an efficient and spectrally accurate method, named as SFSC-IRAM method, for the BdGEs. Ample numerical tests are provided to illustrate the accuracy and efficiency. Finally, we apply the new method to study systematically the excitation spectrum and Bogoliubov amplitudes around the ground state with different parameters in different spatial dimensions.


 

A Fourth Order Finite Difference Method for Solving Elliptic Interface Problems with the FFT Acceleration

Shan Zhao, University of Alabama

In this talk, we will introduce an augmented matched interface and boundary (AMIB) method for solving elliptic interface problems. The AMIB method provides a uniform fictitious domain approach to correct the fourth order central difference near interfaces and boundaries. In treating a smoothly curved interface, zeroth and first order jump conditions are enforced repeatedly to generate fictitious values surrounding the interface. Different types of boundary conditions, including Dirichlet, Neumann, Robin and their mix combinations, can be imposed to generate fictitious values outside boundaries. Based on fictitious values at interfaces and boundaries, the AMIB method reconstructs Cartesian derivative jumps as additional unknowns and forms an enlarged linear system. In the Schur complement solution of such system, the FFT algorithm will not sense the solution discontinuities, so that the discrete Laplacian can be efficiently inverted. In our 2D tests, the FFT-based AMIB not only achieves a fourth order convergence in dealing with interfaces and boundaries, but also produces an overall complexity of O(n^2 log n) for a n-by-n uniform grid. Moreover, the AMIB method can provide fourth order accurate approximations to solution gradients and fluxes.