Topic overview

Computation

1468 works3673 researchers

Open map Browse papers

Map preview

Start with the graph, then narrow the list

1468works

3673researchers

Next steps

Use the topic as a working map

Open the full map for clusters, then return here to scan ranked papers and people.

Inspect nearby papers, researchers, institutions and communities without opening a separate graph page.

Building this map preview

BZPEER is loading the nearby papers, people, topics and institutions for this page.

preprint2016arXiv

Kalman-based Stochastic Gradient Method with Stop Condition and Insensitivity to Conditioning

Modern proximal and stochastic gradient descent (SGD) methods are believed to efficiently minimize large composite objective functions, but such methods have two algorithmic challenges: (1) a lack of fast or justified stop conditions, and (2) sensitivity to the objective function's conditioning. In response to the first challenge, modern proximal and SGD methods guarantee convergence only after multiple epochs, but such a guarantee renders proximal and SGD methods infeasible when the number of component functions is very large or infinite. In response to the second challenge, second order SGD methods have been developed, but they are marred by the complexity of their analysis. In this work, we address these challenges on the limited, but important, linear regression problem by introducing and analyzing a second order proximal/SGD method based on Kalman Filtering (kSGD). Through our analysis, we show kSGD is asymptotically optimal, develop a fast algorithm for very large, infinite or streaming data sources with a justified stop condition, prove that kSGD is insensitive to the problem's conditioning, and develop a unique approach for analyzing the complex second order dynamic

preprint2016arXiv

Bayesian inverse problems with $l_1$ priors: a Randomize-then-Optimize approach

Prior distributions for Bayesian inference that rely on the $l_1$-norm of the parameters are of considerable interest, in part because they promote parameter fields with less regularity than Gaussian priors (e.g., discontinuities and blockiness). These $l_1$-type priors include the total variation (TV) prior and the Besov $B^s_{1,1}$ space prior, and in general yield non-Gaussian posterior distributions. Sampling from these posteriors is challenging, particularly in the inverse problem setting where the parameter space is high-dimensional and the forward problem may be nonlinear. This paper extends the randomize-then-optimize (RTO) method, an optimization-based sampling algorithm developed for Bayesian inverse problems with Gaussian priors, to inverse problems with $l_1$-type priors. We use a variable transformation to convert an $l_1$-type prior to a standard Gaussian prior, such that the posterior distribution of the transformed parameters is amenable to Metropolized sampling via RTO. We demonstrate this approach on several deconvolution problems and an elliptic PDE inverse problem, using TV or Besov $B^s_{1,1}$ space priors. Our results show that the transformed RTO algorithm ch

preprint2016arXiv

Fast symmetric factorization of hierarchical matrices with applications

We present a fast direct algorithm for computing symmetric factorizations, i.e. $A = WW^T$, of symmetric positive-definite hierarchical matrices with weak-admissibility conditions. The computational cost for the symmetric factorization scales as $\mathcal{O}(n \log^2 n)$ for hierarchically off-diagonal low-rank matrices. Once this factorization is obtained, the cost for inversion, application, and determinant computation scales as $\mathcal{O}(n \log n)$. In particular, this allows for the near optimal generation of correlated random variates in the case where $A$ is a covariance matrix. This symmetric factorization algorithm depends on two key ingredients. First, we present a novel symmetric factorization formula for low-rank updates to the identity of the form $I+UKU^T$. This factorization can be computed in $\mathcal{O}(n)$ time if the rank of the perturbation is sufficiently small. Second, combining this formula with a recursive divide-and-conquer strategy, near linear complexity symmetric factorizations for hierarchically structured matrices can be obtained. We present numerical results for matrices relevant to problems in probability \& statistics (Gaussian processes), interp

preprint2016arXiv

Fast Bayesian whole-brain fMRI analysis with spatial 3D priors

Spatial whole-brain Bayesian modeling of task-related functional magnetic resonance imaging (fMRI) is a great computational challenge. Most of the currently proposed methods therefore do inference in subregions of the brain separately or do approximate inference without comparison to the true posterior distribution. A popular such method, which is now the standard method for Bayesian single subject analysis in the SPM software, is introduced in Penny et al. (2005b). The method processes the data slice-by-slice and uses an approximate variational Bayes (VB) estimation algorithm that enforces posterior independence between activity coefficients in different voxels. We introduce a fast and practical Markov chain Monte Carlo (MCMC) scheme for exact inference in the same model, both slice-wise and for the whole brain using a 3D prior on activity coefficients. The algorithm exploits sparsity and uses modern techniques for efficient sampling from high-dimensional Gaussian distributions, leading to speed-ups without which MCMC would not be a practical option. Using MCMC, we are for the first time able to evaluate the approximate VB posterior against the exact MCMC posterior, and show that

preprint2014arXiv

An implementation of a randomized algorithm for principal component analysis

Recent years have witnessed intense development of randomized methods for low-rank approximation. These methods target principal component analysis (PCA) and the calculation of truncated singular value decompositions (SVD). The present paper presents an essentially black-box, fool-proof implementation for Mathworks' MATLAB, a popular software platform for numerical computation. As illustrated via several tests, the randomized algorithms for low-rank approximation outperform or at least match the classical techniques (such as Lanczos iterations) in basically all respects: accuracy, computational efficiency (both speed and memory usage), ease-of-use, parallelizability, and reliability. However, the classical procedures remain the methods of choice for estimating spectral norms, and are far superior for calculating the least singular values and corresponding singular vectors (or singular subspaces).

preprint2017arXiv

Scalable Bayesian variable selection and model averaging under block orthogonal design

We propose a scalable algorithmic framework for exact Bayesian variable selection and model averaging in linear models under the assumption that the Gram matrix is block-diagonal, and as a heuristic for exploring the model space for general designs. In block-diagonal designs our approach returns the most probable model of any given size without resorting to numerical integration. The algorithm also provides a novel and efficient solution to the frequentist best subset selection problem for block-diagonal designs. Posterior probabilities for any number of models are obtained by evaluating a single one-dimensional integral that can be computed upfront, and other quantities of interest such as variable inclusion probabilities and model averaged regression estimates by carrying out an adaptive, deterministic one-dimensional numerical integration. The overall computational cost scales linearly with the number of blocks, which can be processed in parallel, and exponentially with the block size, rendering it most adequate in situations where predictors are organized in many moderately-sized blocks. For general designs, we approximate the Gram matrix by a block-diagonal using spectral clus

preprint2016arXiv

Multivariate distributions with fixed marginals and correlations

Consider the problem of drawing random variates $(X_1,\ldots,X_n)$ from a distribution where the marginal of each $X_i$ is specified, as well as the correlation between every pair $X_i$ and $X_j$. For given marginals, the Fréchet-Hoeffding bounds put a lower and upper bound on the correlation between $X_i$ and $X_j$. Any achievable correlation between $X_i$ and $X_j$ is a convex combinations of these bounds. The value $λ(X_i,X_j) \in [0,1]$ of this convex combination is called here the convexity parameter of $(X_i,X_j),$ with $λ(X_i,X_j) = 1$ corresponding to the upper bound and maximal correlation. For given marginal distributions functions $F_1,\ldots,F_n$ of $(X_1,\ldots,X_n)$ we show that $λ(X_i,X_j) = λ_{ij}$ if and only if there exist symmetric Bernoulli random variables $(B_1,\ldots,B_n)$ (that is $\{0,1\}$ random variables with mean 1/2) such that $λ(B_i,B_j) = λ_{ij}$. In addition, we characterize completely the set of convexity parameters for symmetric Bernoulli marginals in two, three and four dimensions.

preprint2016arXiv

High-dimensional Filtering using Nested Sequential Monte Carlo

Sequential Monte Carlo (SMC) methods comprise one of the most successful approaches to approximate Bayesian filtering. However, SMC without good proposal distributions struggle in high dimensions. We propose nested sequential Monte Carlo (NSMC), a methodology that generalises the SMC framework by requiring only approximate, properly weighted, samples from the SMC proposal distribution, while still resulting in a correct SMC algorithm. This way we can exactly approximate the locally optimal proposal, and extend the class of models for which we can perform efficient inference using SMC. We show improved accuracy over other state-of-the-art methods on several spatio-temporal state space models.

preprint2016arXiv

Error bounds for gradient density estimation computed from a finite sample set using the method of stationary phase

For a twice continuously differentiable function $S$, we define the density function of its gradient (derivative in one dimension) $s = S^{\prime}$ as a random variable transformation of a uniformly distributed random variable using $s$ as the transformation function. Given $N$ values of $S$ sampled at equally spaced locations, we demonstrate using the method of stationary phase that the approximation error between the integral of the scaled, discrete power spectrum of the wave function $ϕ^{D}_τ=\frac{1}{\sqrt{L}}\exp\left(\frac{iS}τ\right)$ and the integral of the true density function of $s$ over an arbitrarily small interval is bounded above by $O(1/N)$ as $N \rightarrow \infty$ ($τ\rightarrow 0$). In addition to its easy implementation and fast computability in $O(N \log N)$ that only requires computing the discrete Fourier transform, our framework for obtaining the derivative density does not involve any parameter selection like the number of histogram bins, width of the histogram bins, width of the kernel parameter, number of mixture components etc. as required by other widely applied methods like histograms and Parzen windows.

preprint2016arXiv

ManifoldOptim: An R Interface to the ROPTLIB Library for Riemannian Manifold Optimization

Manifold optimization appears in a wide variety of computational problems in the applied sciences. In recent statistical methodologies such as sufficient dimension reduction and regression envelopes, estimation relies on the optimization of likelihood functions over spaces of matrices such as the Stiefel or Grassmann manifolds. Recently, Huang, Absil, Gallivan, and Hand (2016) have introduced the library ROPTLIB, which provides a framework and state of the art algorithms to optimize real-valued objective functions over commonly used matrix-valued Riemannian manifolds. This article presents ManifoldOptim, an R package that wraps the C++ library ROPTLIB. ManifoldOptim enables users to access functionality in ROPTLIB through R so that optimization problems can easily be constructed, solved, and integrated into larger R codes. Computationally intensive problems can be programmed with Rcpp and RcppArmadillo, and otherwise accessed through R. We illustrate the practical use of ManifoldOptim through several motivating examples involving dimension reduction and envelope methods in regression.

preprint2016arXiv

ggRandomForests: Exploring Random Forest Survival

Random forest (Leo Breiman 2001a) (RF) is a non-parametric statistical method requiring no distributional assumptions on covariate relation to the response. RF is a robust, nonlinear technique that optimizes predictive accuracy by fitting an ensemble of trees to stabilize model estimates. Random survival forests (RSF) (Ishwaran and Kogalur 2007; Ishwaran et al. 2008) are an extension of Breimans RF techniques allowing efficient nonparametric analysis of time to event data. The randomForestSRC package (Ishwaran and Kogalur 2014) is a unified treatment of Breimans random forest for survival, regression and classification problems. Predictive accuracy makes RF an attractive alternative to parametric models, though complexity and interpretability of the forest hinder wider application of the method. We introduce the ggRandomForests package, tools for visually understand random forest models grown in R (R Core Team 2014) with the randomForestSRC package. The ggRandomForests package is structured to extract intermediate data objects from randomForestSRC objects and generate figures using the ggplot2 (Wickham 2009) graphics package. This document is structured as a tutorial for building r

preprint2016arXiv

Improved Laplace Approximation for Marginal Likelihoods

Statistical applications often involve the calculation of intractable multidimensional integrals. The Laplace formula is widely used to approximate such integrals. However, in high-dimensional or small sample size problems, the shape of the integrand function may be far from that of the Gaussian density, and thus the standard Laplace approximation can be inaccurate. We propose an improved Laplace approximation that reduces the asymptotic error of the standard Laplace formula by one order of magnitude, thus leading to third-order accuracy. We also show, by means of practical examples of various complexity, that the proposed method is extremely accurate, even in high dimensions, improving over the standard Laplace formula. Such examples also demonstrate that the accuracy of the proposed method is comparable with that of other existing methods, which are computationally more demanding. An R implementation of the improved Laplace approximation is also provided through the R package iLaplace available on CRAN.

preprint2015arXiv

ggRandomForests: Visually Exploring a Random Forest for Regression

Random Forests [Breiman:2001] (RF) are a fully non-parametric statistical method requiring no distributional assumptions on covariate relation to the response. RF are a robust, nonlinear technique that optimizes predictive accuracy by fitting an ensemble of trees to stabilize model estimates. The randomForestSRC package (http://cran.r-project.org/package=randomForestSRC) is a unified treatment of Breiman's random forests for survival, regression and classification problems. Predictive accuracy make RF an attractive alternative to parametric models, though complexity and interpretability of the forest hinder wider application of the method. We introduce the ggRandomForests package (http://cran.r-project.org/package=ggRandomForests), for visually understand random forest models grown in R with the randomForestSRC package. The vignette is a tutorial for using the ggRandomForests package with the randomForestSRC package for building and post-processing a regression random forest. In this tutorial, we explore a random forest model for the Boston Housing Data, available in the MASS package. We grow a random forest for regression and demonstrate how ggRandomForests can be used when de

preprint2016arXiv

Geodesic Lagrangian Monte Carlo over the space of positive definite matrices: with application to Bayesian spectral density estimation

We extend the application of Hamiltonian Monte Carlo to allow for sampling from probability distributions defined over symmetric or Hermitian positive definite matrices. To do so, we exploit the Riemannian structure induced by Cartan's century-old canonical metric. The geodesics that correspond to this metric are available in closed-form and---within the context of Lagrangian Monte Carlo---provide a principled way to travel around the space of positive definite matrices. Our method improves Bayesian inference on such matrices by allowing for a broad range of priors, so we are not limited to conjugate priors only. In the context of spectral density estimation, we use the (non-conjugate) complex reference prior as an example modeling option made available by the algorithm. Results based on simulated and real-world multivariate time series are presented in this context, and future directions are outlined.

preprint2016arXiv

Efficient Bayesian computation by proximal Markov chain Monte Carlo: when Langevin meets Moreau

Modern imaging methods rely strongly on Bayesian inference techniques to solve challenging imaging problems. Currently, the predominant Bayesian computation approach is convex optimisation, which scales very efficiently to high dimensional image models and delivers accurate point estimation results. However, in order to perform more complex analyses, for example image uncertainty quantification or model selection, it is necessary to use more computationally intensive Bayesian computation techniques such as Markov chain Monte Carlo methods. This paper presents a new and highly efficient Markov chain Monte Carlo methodology to perform Bayesian computation for high dimensional models that are log-concave and non-smooth, a class of models that is central in imaging sciences. The methodology is based on a regularised unadjusted Langevin algorithm that exploits tools from convex analysis, namely Moreau-Yoshida envelopes and proximal operators, to construct Markov chains with favourable convergence properties. In addition to scaling efficiently to high dimensions, the method is straightforward to apply to models that are currently solved by using proximal optimisation algorithms. We provi

preprint2016arXiv

Variable Selection with Scalable Bootstrap in Generalized Linear Model for Massive Data

Bootstrap is commonly used as a tool for non-parametric statistical inference to estimate meaningful parameters in Variable Selection Models. However, for massive dataset that has exponential growth rate, the computation of Bootstrap Variable Selection (BootVS) can be a crucial issue. In this paper, we propose the method of Variable Selection with Bag of Little Bootstraps (BLBVS) on General Linear Regression and extend it to Generalized Linear Model for selecting important parameters and assessing the quality of estimators' computation efficiency by analyzing results of multiple bootstrap sub-samples. The introduced method best suits large datasets which have parallel and distributed computing structures. To test the performance of BLBVS, we compare it with BootVS from different aspects via empirical studies. The results of simulations show our method has excellent performance. A real data analysis, Risk Forecast of Credit Cards, is also presented to illustrate the computational superiority of BLBVS on large scale datasets, and the result demonstrates the usefulness and validity of our proposed method.

preprint2016arXiv

A Point-process Response Model for Spike Trains from Single Neurons in Neural Circuits under Optogenetic Stimulation

Optogenetics is a new tool to study neuronal circuits that have been genetically modified to allow stimulation by flashes of light. We study recordings from single neurons within neural circuits under optogenetic stimulation. The data from these experiments present a statistical challenge of modeling a high frequency point process (neuronal spikes) while the input is another high frequency point process (light flashes). We further develop a generalized linear model approach to model the relationships between two point processes, employing additive point-process response functions. The resulting model, Point-process Responses for Optogenetics (PRO), provides explicit nonlinear transformations to link the input point process with the output one. Such response functions may provide important and interpretable scientific insights into the properties of the biophysical process that governs neural spiking in response to optogenetic stimulation. We validate and compare the PRO model using a real dataset and simulations, and our model yields a superior area-under-the- curve value as high as 93% for predicting every future spike. For our experiment on the recurrent layer V circuit in the pr

preprint2016arXiv

Localization in High-Dimensional Monte Carlo Filtering

The high dimensionality and computational constraints associated with filtering problems in large-scale geophysical applications are particularly challenging for the Particle Filter (PF). Approximate but efficient methods such as the Ensemble Kalman Filter (EnKF) are therefore usually preferred. A key element of these approximate methods is localization, which is in principle a general technique to avoid the curse of dimensionality and consists in limiting the influence of observations to neighboring sites. However, while it works effectively with the EnKF, localization introduces harmful discontinuities in the estimated physical fields when applied blindly to the PF. In the present paper, we explore two possible local algorithms based on the EnKPF, a hybrid method combining the EnKF and the PF. A simulation study in a conjugate normal setup allows to highlight the trade-offs involved when applying localization to PF type of algorithms in the high-dimensional setting. Experiments with the Lorenz96 model demonstrate the ability of the local EnKPF algorithms to perform well even with a small number of particles compared to the problem size.

preprint2016arXiv

Simple, Scalable and Accurate Posterior Interval Estimation

There is a lack of simple and scalable algorithms for uncertainty quantification. Bayesian methods quantify uncertainty through posterior and predictive distributions, but it is difficult to rapidly estimate summaries of these distributions, such as quantiles and intervals. Variational Bayes approximations are widely used, but may badly underestimate posterior covariance. Typically, the focus of Bayesian inference is on point and interval estimates for one-dimensional functionals of interest. In small scale problems, Markov chain Monte Carlo algorithms remain the gold standard, but such algorithms face major problems in scaling up to big data. Various modifications have been proposed based on parallelization and approximations based on subsamples, but such approaches are either highly complex or lack theoretical support and/or good performance outside of narrow settings. We propose a very simple and general posterior interval estimation algorithm, which is based on running Markov chain Monte Carlo in parallel for subsets of the data and averaging quantiles estimated from each subset. We provide strong theoretical guarantees and illustrate performance in several applications.

preprint2016arXiv

A Note on a Sum of Lognormals

This note considers the applicability of Gauss-Hermite quadrature and direct numerical quadrature for computation of moment generating function (mgf) and the derivatives. A preprocessing using the asymptotic technique is employed while computing the characteristic function (chf) using Gauss Hermite quadrature while this is optional for mgf. The mgf of the low and high amplitude regions of a single lognormal variable and the derivatives is examined and attention is drawn to the effect of variance. The problem of inversion of the mgf/chf of a sum of lognormals to obtain the CDF/pdf is considered with special reference to methods related to Post Widder technique, Gaussian quadrature and the Fourier series method. The method based on the complex exponential integral which makes use of the derivative of the cumulant is an alternative. Segmentation of the mgf/chf on the basis of the derivative structure which indicates activity rate is shown to be useful.

preprint2016arXiv

Fast and Adaptive Sparse Precision Matrix Estimation in High Dimensions

This paper proposes a new method for estimating sparse precision matrices in the high dimensional setting. It has been popular to study fast computation and adaptive procedures for this problem. We propose a novel approach, called Sparse Column-wise Inverse Operator, to address these two issues. We analyze an adaptive procedure based on cross validation, and establish its convergence rate under the Frobenius norm. The convergence rates under other matrix norms are also established. This method also enjoys the advantage of fast computation for large-scale problems, via a coordinate descent algorithm. Numerical merits are illustrated using both simulated and real datasets. In particular, it performs favorably on an HIV brain tissue dataset and an ADHD resting-state fMRI dataset.

preprint2016arXiv

High-Dimensional Bayesian Regularised Regression with the BayesReg Package

Bayesian penalized regression techniques, such as the Bayesian lasso and the Bayesian horseshoe estimator, have recently received a significant amount of attention in the statistics literature. However, software implementing state-of-the-art Bayesian penalized regression, outside of general purpose Markov chain Monte Carlo platforms such as STAN, is relatively rare. This paper introduces bayesreg, a new toolbox for fitting Bayesian penalized regression models with continuous shrinkage prior densities. The toolbox features Bayesian linear regression with Gaussian or heavy-tailed error models and Bayesian logistic regression with ridge, lasso, horseshoe and horseshoe$+$ estimators. The toolbox is free, open-source and available for use with the MATLAB and R numerical platforms.

preprint2016arXiv

CoinCalc -- A new R package for quantifying simultaneities of event series

We present the new R package CoinCalc for performing event coincidence analysis (ECA), a novel statistical method to quantify the simultaneity of events contained in two series of observations, either as simultaneous or lagged coincidences within a user-specific temporal tolerance window. The package also provides different analytical as well as surrogate-based significance tests (valid under different assumptions about the nature of the observed event series) as well as an intuitive visualization of the identified coincidences. We demonstrate the usage of CoinCalc based on two typical geoscientific example problems addressing the relationship between meteorological extremes and plant phenology as well as that between soil properties and land cover.

preprint2016arXiv

The Reduced-Order Hybrid Monte Carlo Sampling Smoother

Hybrid Monte-Carlo (HMC) sampling smoother is a fully non-Gaussian four-dimensional data assimilation algorithm that works by directly sampling the posterior distribution formulated in the Bayesian framework. The smoother in its original formulation is computationally expensive due to the intrinsic requirement of running the forward and adjoint models repeatedly. Here we present computationally efficient versions of the HMC sampling smoother based on reduced-order approximations of the underlying model dynamics. The schemes developed herein are tested numerically using the shallow-water equations model on Cartesian coordinates. The results reveal that the reduced-order versions of the smoother are capable of accurately capturing the posterior probability density, while being significantly faster than the original full order formulation.

339 works