75,981 research outputs found
Fast algorithms for hyperspectral Diffuse Optical Tomography
The image reconstruction of chromophore concentrations using Diffuse Optical
Tomography (DOT) data can be described mathematically as an ill-posed inverse
problem. Recent work has shown that the use of hyperspectral DOT data, as
opposed to data sets comprising of a single or, at most, a dozen wavelengths,
has the potential for improving the quality of the reconstructions. The use of
hyperspectral diffuse optical data in the formulation and solution of the
inverse problem poses a significant computational burden. The forward operator
is, in actuality, nonlinear. However, under certain assumptions, a linear
approximation, called the Born approximation, provides a suitable surrogate for
the forward operator, and we assume this to be true in the present work.
Computation of the Born matrix requires the solution of thousands of large
scale discrete PDEs and the reconstruction problem, requires matrix-vector
products with the (dense) Born matrix. In this paper, we address both of these
difficulties, thus making the Born approach a computational viable approach for
hyDOT reconstruction. In this paper, we assume that the images we wish to
reconstruct are anomalies of unknown shape and constant value, described using
a parametric level set approach, (PaLS) on a constant background. Specifically,
to address the issue of the PDE solves, we develop a novel recycling-based
Krylov subspace approach that leverages certain system similarities across
wavelengths. To address expense of using the Born operator in the inversion, we
present a fast algorithm for compressing the Born operator that locally
compresses across wavelengths for a given source-detector set and then
recursively combines the low-rank factors to provide a global low-rank
approximation. This low-rank approximation can be used implicitly to speed up
the recovery of the shape parameters and the chromophore concentrations
Fast Direct Methods for Gaussian Processes
A number of problems in probability and statistics can be addressed using the
multivariate normal (Gaussian) distribution. In the one-dimensional case,
computing the probability for a given mean and variance simply requires the
evaluation of the corresponding Gaussian density. In the -dimensional
setting, however, it requires the inversion of an covariance
matrix, , as well as the evaluation of its determinant, . In many
cases, such as regression using Gaussian processes, the covariance matrix is of
the form , where is computed using a specified
covariance kernel which depends on the data and additional parameters
(hyperparameters). The matrix is typically dense, causing standard direct
methods for inversion and determinant evaluation to require
work. This cost is prohibitive for large-scale modeling. Here, we show that for
the most commonly used covariance functions, the matrix can be
hierarchically factored into a product of block low-rank updates of the
identity matrix, yielding an algorithm for inversion.
More importantly, we show that this factorization enables the evaluation of the
determinant , permitting the direct calculation of probabilities in
high dimensions under fairly broad assumptions on the kernel defining . Our
fast algorithm brings many problems in marginalization and the adaptation of
hyperparameters within practical reach using a single CPU core. The combination
of nearly optimal scaling in terms of problem size with high-performance
computing resources will permit the modeling of previously intractable
problems. We illustrate the performance of the scheme on standard covariance
kernels
HLIBCov: Parallel Hierarchical Matrix Approximation of Large Covariance Matrices and Likelihoods with Applications in Parameter Identification
We provide more technical details about the HLIBCov package, which is using
parallel hierarchical (\H-) matrices to identify unknown parameters of the
covariance function (variance, smoothness, and covariance length). These
parameters are estimated by maximizing the joint Gaussian log-likelihood
function. The HLIBCov package approximates large dense inhomogeneous covariance
matrices with a log-linear computational cost and storage requirement. We
explain how to compute the Cholesky factorization, determinant, inverse and
quadratic form in the H-matrix format. To demonstrate the numerical
performance, we identify three unknown parameters in an example with 2,000,000
locations on a PC-desktop.Comment: 26 pages, 11 figure
A Fast Block Low-Rank Dense Solver with Applications to Finite-Element Matrices
This article presents a fast solver for the dense "frontal" matrices that
arise from the multifrontal sparse elimination process of 3D elliptic PDEs. The
solver relies on the fact that these matrices can be efficiently represented as
a hierarchically off-diagonal low-rank (HODLR) matrix. To construct the
low-rank approximation of the off-diagonal blocks, we propose a new
pseudo-skeleton scheme, the boundary distance low-rank approximation, that
picks rows and columns based on the location of their corresponding vertices in
the sparse matrix graph. We compare this new low-rank approximation method to
the adaptive cross approximation (ACA) algorithm and show that it achieves
betters speedup specially for unstructured meshes. Using the HODLR direct
solver as a preconditioner (with a low tolerance) to the GMRES iterative
scheme, we can reach machine accuracy much faster than a conventional LU
solver. Numerical benchmarks are provided for frontal matrices arising from 3D
finite element problems corresponding to a wide range of applications
Fast hierarchical solvers for sparse matrices using extended sparsification and low-rank approximation
Inversion of sparse matrices with standard direct solve schemes is robust,
but computationally expensive. Iterative solvers, on the other hand,
demonstrate better scalability; but, need to be used with an appropriate
preconditioner (e.g., ILU, AMG, Gauss-Seidel, etc.) for proper convergence. The
choice of an effective preconditioner is highly problem dependent. We propose a
novel fully algebraic sparse matrix solve algorithm, which has linear
complexity with the problem size. Our scheme is based on the Gauss elimination.
For a given matrix, we approximate the LU factorization with a tunable accuracy
determined a priori. This method can be used as a stand-alone direct solver
with linear complexity and tunable accuracy, or it can be used as a black-box
preconditioner in conjunction with iterative methods such as GMRES. The
proposed solver is based on the low-rank approximation of fill-ins generated
during the elimination. Similar to H-matrices, fill-ins corresponding to blocks
that are well-separated in the adjacency graph are represented via a
hierarchical structure. The linear complexity of the algorithm is guaranteed if
the blocks corresponding to well-separated clusters of variables are
numerically low-rank
Likelihood Approximation With Hierarchical Matrices For Large Spatial Datasets
We use available measurements to estimate the unknown parameters (variance,
smoothness parameter, and covariance length) of a covariance function by
maximizing the joint Gaussian log-likelihood function. To overcome cubic
complexity in the linear algebra, we approximate the discretized covariance
function in the hierarchical (H-) matrix format. The H-matrix format has a
log-linear computational cost and storage , where the rank is
a small integer and is the number of locations. The H-matrix technique
allows us to work with general covariance matrices in an efficient way, since
H-matrices can approximate inhomogeneous covariance functions, with a fairly
general mesh that is not necessarily axes-parallel, and neither the covariance
matrix itself nor its inverse have to be sparse. We demonstrate our method with
Monte Carlo simulations and an application to soil moisture data. The C, C++
codes and data are freely available.Comment: 23 pages, 24 figures, 3 tables, version after the second major
revisio
An efficient method for block low-rank approximations for kernel matrix systems
In the iterative solution of dense linear systems from boundary integral
equations or systems involving kernel matrices, the main challenges are the
expensive matrix-vector multiplication and the storage cost which are usually
tackled by hierarchical matrix techniques such as and
matrices. However, hierarchical matrices also have a high
construction cost that is dominated by the low-rank approximations of the
sub-blocks of the kernel matrix. In this paper, an efficient method is proposed
to give a low-rank approximation of the kernel matrix block in
the form of an interpolative decomposition (ID) for a kernel function
and two properly located point sets . The proposed method combines
the ID using strong rank-revealing QR (sRRQR), which is purely algebraic, with
analytic kernel information to reduce the construction cost of a rank-
approximation from , for ID using sRRQR alone, to
which is not related to . Numerical experiments show that
matrix construction with the proposed algorithm only requires a
computational cost linear in the matrix dimension
Tucker Tensor analysis of Matern functions in spatial statistics
In this work, we describe advanced numerical tools for working with
multivariate functions and for the analysis of large data sets. These tools
will drastically reduce the required computing time and the storage cost, and,
therefore, will allow us to consider much larger data sets or finer meshes.
Covariance matrices are crucial in spatio-temporal statistical tasks, but are
often very expensive to compute and store, especially in 3D. Therefore, we
approximate covariance functions by cheap surrogates in a low-rank tensor
format. We apply the Tucker and canonical tensor decompositions to a family of
Matern- and Slater-type functions with varying parameters and demonstrate
numerically that their approximations exhibit exponentially fast convergence.
We prove the exponential convergence of the Tucker and canonical approximations
in tensor rank parameters. Several statistical operations are performed in this
low-rank tensor format, including evaluating the conditional covariance matrix,
spatially averaged estimation variance, computing a quadratic form,
determinant, trace, loglikelihood, inverse, and Cholesky decomposition of a
large covariance matrix. Low-rank tensor approximations reduce the computing
and storage costs essentially. For example, the storage cost is reduced from an
exponential to a linear scaling , where
is the spatial dimension, is the number of mesh points in one
direction, and is the tensor rank. Prerequisites for applicability of the
proposed techniques are the assumptions that the data, locations, and
measurements lie on a tensor (axes-parallel) grid and that the covariance
function depends on a distance, .Comment: 23 pages, 2 diagrams, 2 tables, 9 figure
Linear Multiple Low-Rank Kernel Based Stationary Gaussian Processes Regression for Time Series
Gaussian processes (GP) for machine learning have been studied systematically
over the past two decades and they are by now widely used in a number of
diverse applications. However, GP kernel design and the associated
hyper-parameter optimization are still hard and to a large extend open
problems. In this paper, we consider the task of GP regression for time series
modeling and analysis. The underlying stationary kernel can be approximated
arbitrarily close by a new proposed grid spectral mixture (GSM) kernel, which
turns out to be a linear combination of low-rank sub-kernels. In the case where
a large number of the sub-kernels are used, either the Nystr\"{o}m or the
random Fourier feature approximations can be adopted to deal efficiently with
the computational demands. The unknown GP hyper-parameters consist of the
non-negative weights of all sub-kernels as well as the noise variance; their
estimation is performed via the maximum-likelihood (ML) estimation framework.
Two efficient numerical optimization methods for solving the unknown
hyper-parameters are derived, including a sequential majorization-minimization
(MM) method and a non-linearly constrained alternating direction of multiplier
method (ADMM). The MM matches perfectly with the proven low-rank property of
the proposed GSM sub-kernels and turns out to be a part of efficiency, stable,
and efficient solver, while the ADMM has the potential to generate better local
minimum in terms of the test MSE. Experimental results, based on various
classic time series data sets, corroborate that the proposed GSM kernel-based
GP regression model outperforms several salient competitors of similar kind in
terms of prediction mean-squared-error and numerical stability.Comment: 15 pages, 5 figures, submitte
Krylov Subspace Recycling for Fast Iterative Least-Squares in Machine Learning
Solving symmetric positive definite linear problems is a fundamental
computational task in machine learning. The exact solution, famously, is
cubicly expensive in the size of the matrix. To alleviate this problem, several
linear-time approximations, such as spectral and inducing-point methods, have
been suggested and are now in wide use. These are low-rank approximations that
choose the low-rank space a priori and do not refine it over time. While this
allows linear cost in the data-set size, it also causes a finite, uncorrected
approximation error. Authors from numerical linear algebra have explored ways
to iteratively refine such low-rank approximations, at a cost of a small number
of matrix-vector multiplications. This idea is particularly interesting in the
many situations in machine learning where one has to solve a sequence of
related symmetric positive definite linear problems. From the machine learning
perspective, such deflation methods can be interpreted as transfer learning of
a low-rank approximation across a time-series of numerical tasks. We study the
use of such methods for our field. Our empirical results show that, on
regression and classification problems of intermediate size, this approach can
interpolate between low computational cost and numerical precision
- …