151 research outputs found
Low-rank updates and a divide-and-conquer method for linear matrix equations
Linear matrix equations, such as the Sylvester and Lyapunov equations, play
an important role in various applications, including the stability analysis and
dimensionality reduction of linear dynamical control systems and the solution
of partial differential equations. In this work, we present and analyze a new
algorithm, based on tensorized Krylov subspaces, for quickly updating the
solution of such a matrix equation when its coefficients undergo low-rank
changes. We demonstrate how our algorithm can be utilized to accelerate the
Newton method for solving continuous-time algebraic Riccati equations. Our
algorithm also forms the basis of a new divide-and-conquer approach for linear
matrix equations with coefficients that feature hierarchical low-rank
structure, such as HODLR, HSS, and banded matrices. Numerical experiments
demonstrate the advantages of divide-and-conquer over existing approaches, in
terms of computational time and memory consumption
TMB: Automatic Differentiation and Laplace Approximation
TMB is an open source R package that enables quick implementation of complex
nonlinear random effect (latent variable) models in a manner similar to the
established AD Model Builder package (ADMB, admb-project.org). In addition, it
offers easy access to parallel computations. The user defines the joint
likelihood for the data and the random effects as a C++ template function,
while all the other operations are done in R; e.g., reading in the data. The
package evaluates and maximizes the Laplace approximation of the marginal
likelihood where the random effects are automatically integrated out. This
approximation, and its derivatives, are obtained using automatic
differentiation (up to order three) of the joint likelihood. The computations
are designed to be fast for problems with many random effects (~10^6) and
parameters (~10^3). Computation times using ADMB and TMB are compared on a
suite of examples ranging from simple models to large spatial models where the
random effects are a Gaussian random field. Speedups ranging from 1.5 to about
100 are obtained with increasing gains for large problems. The package and
examples are available at http://tmb-project.org
Efficient approximation of functions of some large matrices by partial fraction expansions
Some important applicative problems require the evaluation of functions
of large and sparse and/or \emph{localized} matrices . Popular and
interesting techniques for computing and , where
is a vector, are based on partial fraction expansions. However,
some of these techniques require solving several linear systems whose matrices
differ from by a complex multiple of the identity matrix for computing
or require inverting sequences of matrices with the same
characteristics for computing . Here we study the use and the
convergence of a recent technique for generating sequences of incomplete
factorizations of matrices in order to face with both these issues. The
solution of the sequences of linear systems and approximate matrix inversions
above can be computed efficiently provided that shows certain decay
properties. These strategies have good parallel potentialities. Our claims are
confirmed by numerical tests
Solving rank structured Sylvester and Lyapunov equations
We consider the problem of efficiently solving Sylvester and Lyapunov
equations of medium and large scale, in case of rank-structured data, i.e.,
when the coefficient matrices and the right-hand side have low-rank
off-diagonal blocks. This comprises problems with banded data, recently studied
by Haber and Verhaegen in "Sparse solution of the Lyapunov equation for
large-scale interconnected systems", Automatica, 2016, and by Palitta and
Simoncini in "Numerical methods for large-scale Lyapunov equations with
symmetric banded data", SISC, 2018, which often arise in the discretization of
elliptic PDEs.
We show that, under suitable assumptions, the quasiseparable structure is
guaranteed to be numerically present in the solution, and explicit novel
estimates of the numerical rank of the off-diagonal blocks are provided.
Efficient solution schemes that rely on the technology of hierarchical
matrices are described, and several numerical experiments confirm the
applicability and efficiency of the approaches. We develop a MATLAB toolbox
that allows easy replication of the experiments and a ready-to-use interface
for the solvers. The performances of the different approaches are compared, and
we show that the new methods described are efficient on several classes of
relevant problems
SlabLU: A Two-Level Sparse Direct Solver for Elliptic PDEs
The paper describes a sparse direct solver for the linear systems that arise
from the discretization of an elliptic PDE on a two dimensional domain. The
solver is designed to reduce communication costs and perform well on GPUs; it
uses a two-level framework, which is easier to implement and optimize than
traditional multi-frontal schemes based on hierarchical nested dissection
orderings. The scheme decomposes the domain into thin subdomains, or "slabs".
Within each slab, a local factorization is executed that exploits the geometry
of the local domain. A global factorization is then obtained through the LU
factorization of a block-tridiagonal reduced coefficient matrix. The solver has
complexity for the factorization step, and for each
solve once the factorization is completed.
The solver described is compatible with a range of different local
discretizations, and numerical experiments demonstrate its performance for
regular discretizations of rectangular and curved geometries. The technique
becomes particularly efficient when combined with very high-order convergent
multi-domain spectral collocation schemes. With this discretization, a
Helmholtz problem on a domain of size (for
which N=100 \mbox{M}) is solved in 15 minutes to 6 correct digits on a
high-powered desktop with GPU acceleration
Integrated Nested Laplace Approximations for Large-Scale Spatial-Temporal Bayesian Modeling
Bayesian inference tasks continue to pose a computational challenge. This
especially holds for spatial-temporal modeling where high-dimensional latent
parameter spaces are ubiquitous. The methodology of integrated nested Laplace
approximations (INLA) provides a framework for performing Bayesian inference
applicable to a large subclass of additive Bayesian hierarchical models. In
combination with the stochastic partial differential equations (SPDE) approach
it gives rise to an efficient method for spatial-temporal modeling. In this
work we build on the INLA-SPDE approach, by putting forward a performant
distributed memory variant, INLA-DIST, for large-scale applications. To perform
the arising computational kernel operations, consisting of Cholesky
factorizations, solving linear systems, and selected matrix inversions, we
present two numerical solver options, a sparse CPU-based library and a novel
blocked GPU-accelerated approach which we propose. We leverage the recurring
nonzero block structure in the arising precision (inverse covariance) matrices,
which allows us to employ dense subroutines within a sparse setting. Both
versions of INLA-DIST are highly scalable, capable of performing inference on
models with millions of latent parameters. We demonstrate their accuracy and
performance on synthetic as well as real-world climate dataset applications.Comment: 22 pages, 14 figure
Solution of partial differential equations on vector and parallel computers
The present status of numerical methods for partial differential equations on vector and parallel computers was reviewed. The relevant aspects of these computers are discussed and a brief review of their development is included, with particular attention paid to those characteristics that influence algorithm selection. Both direct and iterative methods are given for elliptic equations as well as explicit and implicit methods for initial boundary value problems. The intent is to point out attractive methods as well as areas where this class of computer architecture cannot be fully utilized because of either hardware restrictions or the lack of adequate algorithms. Application areas utilizing these computers are briefly discussed
The automatic solution of partial differential equations using a global spectral method
A spectral method for solving linear partial differential equations (PDEs)
with variable coefficients and general boundary conditions defined on
rectangular domains is described, based on separable representations of partial
differential operators and the one-dimensional ultraspherical spectral method.
If a partial differential operator is of splitting rank , such as the
operator associated with Poisson or Helmholtz, the corresponding PDE is solved
via a generalized Sylvester matrix equation, and a bivariate polynomial
approximation of the solution of degree is computed in
operations. Partial differential operators of
splitting rank are solved via a linear system involving a block-banded
matrix in operations. Numerical
examples demonstrate the applicability of our 2D spectral method to a broad
class of PDEs, which includes elliptic and dispersive time-evolution equations.
The resulting PDE solver is written in MATLAB and is publicly available as part
of CHEBFUN. It can resolve solutions requiring over a million degrees of
freedom in under seconds. An experimental implementation in the Julia
language can currently perform the same solve in seconds.Comment: 22 page
- …