Search

Stochastic Rounding is a probabilistic rounding mode that is surprisingly effective in large-scale computations and low-precision arithmetic. Its random nature promotes error cancellation rather than error accumulation, resulting in slower growth of roundoff errors as the problem size increases, especially when compared to traditional deterministic rounding methods, such as rounding-to-nearest. We advocate for SR as a foundational tool in the complexity analysis of algorithms, and suggest several research directions., Comment: 5 pages

The notion of `stable rank' of a matrix is central to the analysis of randomized matrix algorithms, covariance estimation, deep neural networks, and recommender systems. We compare the properties of the stable rank and intrinsic dimension of real and complex matrices to those of the classical rank. Basic proofs and examples illustrate that the stable rank does not satisfy any of the fundamental rank properties, while the intrinsic dimension satisfies a few. In particular, the stable rank and intrinsic dimension of a submatrix can exceed those of the original matrix; adding a Hermitian positive semi-definite matrix can lower the intrinsic dimension of the sum; and multiplication by a nonsingular matrix can drastically change the stable rank and the intrinsic dimension. We generalize the concept of stable rank to the p-stable in any Schatten p-norm, thereby unifying the concepts of stable rank and intrinsic dimension: The stable rank is the 2-stable rank, while the intrinsic dimension is the 1-stable rank of a Hermitian positive semi-definite matrix. We derive sum and product inequalities for the pth root of the p-stable rank, and show that it is well-conditioned in the norm-wise absolute sense. The conditioning improves if the matrix and the perturbation are Hermitian positive semi-definite.

We a present and analyze rpCholesky-QR, a randomized preconditioned Cholesky-QR algorithm for computing the thin QR factorization of real mxn matrices with rank n. rpCholesky-QR has a low orthogonalization error, a residual on the order of machine precision, and does not break down for highly singular matrices. We derive rigorous and interpretable two-norm perturbation bounds for rpCholesky-QR that require a minimum of assumptions. Numerical experiments corroborate the accuracy of rpCholesky-QR for preconditioners sampled from as few as 3n rows, and illustrate that the two-norm deviation from orthonormality increases with only the condition number of the preconditioned matrix, rather than its square -- even if the original matrix is numerically singular.

Motivated by the popularity of stochastic rounding in the context of machine learning and the training of large-scale deep neural network models, we consider stochastic nearness rounding of real matrices $\mathbf{A}$ with many more rows than columns. We provide novel theoretical evidence, supported by extensive experimental evaluation that, with high probability, the smallest singular value of a stochastically rounded matrix is well bounded away from zero -- regardless of how close $\mathbf{A}$ is to being rank deficient and even if $\mathbf{A}$ is rank-deficient. In other words, stochastic rounding \textit{implicitly regularizes} tall and skinny matrices $\mathbf{A}$ so that the rounded version has full column rank. Our proofs leverage powerful results in random matrix theory, and the idea that stochastic rounding errors do not concentrate in low-dimensional column spaces.

We perturb a real matrix $A$ of full column rank, and derive lower bounds for the smallest singular values of the perturbed matrix, in terms of normwise absolute perturbations. Our bounds, which extend existing lower-order expressions, demonstrate the potential increase in the smallest singular values, and represent a qualitative model for the increase in the small singular values after a matrix has been downcast to a lower arithmetic precision. Numerical experiments confirm the qualitative validity of this model and its ability to predict singular values changes in the presence of decreased arithmetic precision.

We analyse the calibration of BayesCG under the Krylov prior, a probabilistic numeric extension of the Conjugate Gradient (CG) method for solving systems of linear equations with symmetric positive definite coefficient matrix. Calibration refers to the statistical quality of the posterior covariances produced by a solver. Since BayesCG is not calibrated in the strict existing notion, we propose instead two test statistics that are necessary but not sufficient for calibration: the Z-statistic and the new S-statistic. We show analytically and experimentally that under low-rank approximate Krylov posteriors, BayesCG exhibits desirable properties of a calibrated solver, is only slightly optimistic, and is computationally competitive with CG., Comment: 40 Pages

We advocate a numerically reliable and accurate approach for practical parameter identifiability analysis: Applying column subset selection (CSS) to the sensitivity matrix, instead of computing an eigenvalue decomposition of the Fischer information matrix. Identifiability analysis via CSS has three advantages: (i) It quantifies reliability of the subsets of parameters selected as identifiable and unidentifiable. (ii) It establishes criteria for comparing the accuracy of different algorithms. (iii) The implementations are numerically more accurate and reliable than eigenvalue methods applied to the Fischer matrix, yet without an increase in computational cost. The effectiveness of the CSS methods is illustrated with extensive numerical experiments on sensitivity matrices from six physical models, as well as on adversarial synthetic matrices. Among the CSS methods, we recommend an implementation based on the strong rank-revealing QR algorithm because of its rigorous accuracy guarantees for both identifiable and non-identifiable parameters.

We analyze the forward error in the floating point summation of real numbers, for computations in low precision or extreme-scale problem dimensions that push the limits of the precision. We present a systematic recurrence for a martingale on a computational tree, which leads to explicit and interpretable bounds without asymptotic big-O terms. Two probability parameters strengthen the precision-awareness of our bounds: one parameter controls the first order terms in the summation error, while the second one is designed for controlling higher order terms in low precision or extreme-scale problem dimensions. Our systematic approach yields new deterministic and probabilistic error bounds for three classes of mono-precision algorithms: general summation, shifted general summation, and compensated (sequential) summation. Extension of our systematic error analysis to mixed-precision summation algorithms that allow any number of precisions yields the first probabilistic bounds for the mixed-precision FABsum algorithm. Numerical experiments illustrate that the probabilistic bounds are accurate, and that among the three classes of mono-precision algorithms, compensated summation is generally the most accurate. As for mixed precision algorithms, our recommendation is to minimize the magnitude of intermediate partial sums relative to the precision in which they are computed.

For real symmetric matrices that are accessible only through matrix vector products, we present Monte Carlo estimators for computing the diagonal elements. Our probabilistic bounds for normwise absolute and relative errors apply to Monte Carlo estimators based on random Rademacher, sparse Rademacher, normalized and unnormalized Gaussian vectors, and to vectors with bounded fourth moments. The novel use of matrix concentration inequalities in our proofs represents a systematic model for future analyses. Our bounds mostly do not depend on the matrix dimension, target different error measures than existing work, and imply that the accuracy of the estimators increases with the diagonal dominance of the matrix. An application to derivative-based global sensitivity metrics corroborates this, as do numerical experiments on synthetic test matrices. We recommend against the use in practice of sparse Rademacher vectors, which are the basis for many randomized sketching and sampling algorithms, because they tend to deliver barely a digit of accuracy even under large sampling amounts.

We analyse the forward error in the floating point summation of real numbers, from algorithms that do not require recourse to higher precision or better hardware. We derive informative explicit expressions, and new deterministic and probabilistic bounds for errors in three classes of algorithms: general summation,shifted general summation, and compensated (sequential) summation. Our probabilistic bounds for general and shifted general summation hold to all orders. For compensated summation, we also present deterministic and probabilistic first and second order bounds, with a first order bound that differs from existing ones. Numerical experiments illustrate that the bounds are informative and that among the three algorithm classes, compensated summation is generally the most accurate method.

The explosion of biobank data offers immediate opportunities for gene-environment (GxE) interaction studies of complex diseases because of the large sample sizes and the rich collection in genetic and non-genetic information. However, the extremely large sample size also introduces new computational challenges in GxE assessment, especially for set-based GxE variance component (VC) tests, which are a widely used strategy to boost overall GxE signals and to evaluate the joint GxE effect of multiple variants from a biologically meaningful unit (e.g., gene). In this work, we focus on continuous traits and present SEAGLE, a Scalable Exact AlGorithm for Large-scale set-based GxE tests, to permit GxE VC tests for biobank-scale data. SEAGLE employs modern matrix computations to achieve the same "exact" results as the original GxE VC tests without imposing additional assumptions or relying on approximations. SEAGLE can easily accommodate sample sizes in the order of $10^5$, is implementable on standard laptops, and does not require specialized computing equipment. We demonstrate SEAGLE's performance through extensive simulations. We illustrate its utility by conducting genome-wide gene-based GxE analysis on the Taiwan Biobank data to explore the interaction of gene and physical activity status on body mass index.

This paper presents a probabilistic perspective on iterative methods for approximating the solution $\mathbf{x}_* \in \mathbb{R}^d$ of a nonsingular linear system $\mathbf{A} \mathbf{x}_* = \mathbf{b}$. In the approach a standard iterative method on $\mathbb{R}^d$ is lifted to act on the space of probability distributions $\mathcal{P}(\mathbb{R}^d)$. Classically, an iterative method produces a sequence $\mathbf{x}_m$ of approximations that converge to $\mathbf{x}_*$. The output of the iterative methods proposed in this paper is, instead, a sequence of probability distributions $\mu_m \in \mathcal{P}(\mathbb{R}^d)$. The distributional output both provides a "best guess" for $\mathbf{x}_*$, for example as the mean of $\mu_m$, and also probabilistic uncertainty quantification for the value of $\mathbf{x}_*$ when it has not been exactly determined. Theoretical analysis is provided in the prototypical case of a stationary linear iterative method. In this setting we characterise both the rate of contraction of $\mu_m$ to an atomic measure on $\mathbf{x}_*$ and the nature of the uncertainty quantification being provided. We conclude with an empirical illustration that highlights the insight into solution uncertainty that can be provided by probabilistic iterative methods.

The Bayesian Conjugate Gradient method (BayesCG) is a probabilistic generalization of the Conjugate Gradient method (CG) for solving linear systems with real symmetric positive definite coefficient matrices. Our CG-based implementation of BayesCG under a structure-exploiting prior distribution represents an 'uncertainty-aware' version of CG. Its output consists of CG iterates and posterior covariances that can be propagated to subsequent computations. The covariances have low-rank and are maintained in factored form. This allows easy generation of accurate samples to probe uncertainty in downstream computations. Numerical experiments confirm the effectiveness of the low-rank posterior covariances., Comment: 34 Pages including supplementary material (main paper is 23 pages, supplement is 11 pages). Computer codes are available at https://github.com/treid5/ProbNumCG_Supp

Multivariate multiple linear regression (MMLR), which occurs in a number of practical applications, generalizes traditional least squares (multivariate linear regression) to multiple right-hand sides. We extend recent MLR analyses to sketched MMLR in general Schatten $p$-norms by interpreting the sketched problem as a multiplicative perturbation. Our work represents an extension of Maher's results on Schatten $p$-norms. We derive expressions for the exact and perturbed solutions in terms of projectors for easy geometric interpretation. We also present a geometric interpretation of the action of the sketching matrix in terms of relevant subspaces. We show that a key term in assessing the accuracy of the sketched MMLR solution can be viewed as a tangent of a largest principal angle between subspaces under some assumptions. Our results enable additional interpretation of the difference between an orthogonal and oblique projector with the same range.

Probabilistic models are proposed for bounding the forward error in the numerically computed inner product (dot product, scalar product) between of two real $n$-vectors. We derive probabilistic perturbation bounds, as well as probabilistic roundoff error bounds for the sequential accumulation of the inner product. These bounds are non-asymptotic, explicit, and make minimal assumptions on perturbations and roundoffs. The perturbations are represented as independent, bounded, zero-mean random variables, and the probabilistic perturbation bound is based on Azuma's inequality. The roundoffs are also represented as bounded, zero-mean random variables. The first probabilistic bound assumes that the roundoffs are independent, while the second one does not. For the latter, we construct a Martingale that mirrors the sequential order of computations. Numerical experiments confirm that our bounds are more informative, often by several orders of magnitude, than traditional deterministic bounds -- even for small vector dimensions~$n$ and very stringent success probabilities. In particular the probabilistic roundoff error bounds are functions of $\sqrt{n}$ rather than~$n$, thus giving a quantitative confirmation of Wilkinson's intuition. The paper concludes with a critical assessment of the probabilistic approach.

Several recent works have developed a new, probabilistic interpretation for numerical algorithms solving linear systems in which the solution is inferred in a Bayesian framework, either directly or by inferring the unknown action of the matrix inverse. These approaches have typically focused on replicating the behavior of the conjugate gradient method as a prototypical iterative method. In this work surprisingly general conditions for equivalence of these disparate methods are presented. We also describe connections between probabilistic linear solvers and projection methods for linear systems, providing a probabilistic interpretation of a far more general class of iterative methods. In particular, this provides such an interpretation of the generalised minimum residual method. A probabilistic view of preconditioning is also introduced. These developments unify the literature on probabilistic linear solvers, and provide foundational connections to the literature on iterative solvers for linear systems.

Linear regression is a classic method of data analysis. In recent years, sketching -- a method of dimension reduction using random sampling, random projections, or both -- has gained popularity as an effective computational approximation when the number of observations greatly exceeds the number of variables. In this paper, we address the following question: How does sketching affect the statistical properties of the solution and key quantities derived from it? To answer this question, we present a projector-based approach to sketched linear regression that is exact and that requires minimal assumptions on the sketching matrix. Therefore, downstream analyses hold exactly and generally for all sketching schemes. Additionally, a projector-based approach enables derivation of key quantities from classic linear regression that account for the combined model- and algorithm-induced uncertainties. We demonstrate the usefulness of a projector-based approach in quantifying and enabling insight on excess uncertainties and bias-variance decompositions for sketched linear regression. Finally, we demonstrate how the insights from our projector-based analyses can be used to produce practical sketching diagnostics to aid the design of judicious sketching schemes.

Given a real symmetric positive semi-definite matrix E, and an approximation S that is a sum of n independent matrix-valued random variables, we present bounds on the relative error in S due to randomization. The bounds do not depend on the matrix dimensions but only on the numerical rank (intrinsic dimension) of E. Our approach resembles the low-rank approximation of kernel matrices from random features, but our accuracy measures are more stringent. In the context of parameter selection based on active subspaces, where S is computed via Monte Carlo sampling, we present a bound on the number of samples so that with high probability the angle between the dominant subspaces of E and S is less than a user-specified tolerance. This is a substantial improvement over existing work, as it is a non-asymptotic and fully explicit bound on the sampling amount n, and it allows the user to tune the success probability. It also suggests that Monte Carlo sampling can be efficient in the presence of many parameters, as long as the underlying function f is sufficiently smooth.

This is a systematic investigation into the sensitivity of low-rank approximations of real matrices. We show that the low-rank approximation errors, in the two-norm, Frobenius norm and more generally, any Schatten p-norm, are insensitive to additive rank-preserving perturbations in the projector basis; and to matrix perturbations that are additive or change the number of columns (including multiplicative perturbations). Thus, low-rank matrix approximations are always well-posed and do not require a singular value gap. In the presence of a singular value gap, connections are established between low-rank approximations and subspace angles.

We present randomized algorithms for estimating the trace and deter- minant of Hermitian positive semi-definite matrices. The algorithms are based on subspace iteration, and access the matrix only through matrix vector products. We analyse the error due to randomization, for starting guesses whose elements are Gaussian or Rademacher random variables. The analysis is cleanly separated into a structural (deterministic) part followed by a probabilistic part. Our absolute bounds for the expectation and concentration of the estimators are non-asymptotic and informative even for matrices of low dimension. For the trace estimators, we also present asymptotic bounds on the number of samples (columns of the starting guess) required to achieve a user-specified relative error. Numerical experiments illustrate the performance of the estimators and the tightness of the bounds on low-dimensional matrices; and on a challenging application in uncertainty quantification arising from Bayesian optimal experimental design., Comment: 38 pages. Accepted for publication in Numerische Mathematik

Rue and Held (2005) proposed a method for efficiently computing the Gaussian likelihood for stationary Markov random field models, when the data locations fall on a complete regular grid, and the model has no additive error term. The calculations rely on the availability of the covariances. We prove a theorem giving the rate of convergence of a spectral method of computing the covariances, establishing that the error decays faster than any polynomial in the size of the computing grid. We extend the exact likelihood calculations to the case of non-rectangular domains and missing values on the interior of the grid and to the case when an additive uncorrelated error term (nugget) is present in the model. We also give an alternative formulation of the likelihood that has a smaller memory burden, parts of which can be computed in parallel. We show in simulations that using the exact likelihood can give far better parameter estimates than using standard Markov random field approximations. Having access to the exact likelihood allows for model comparisons via likelihood ratios on large datasets, so as an application of the methods, we compare several state-of-the-art methods for large spatial datasets on an aerosol optical thickness dataset. We find that simple block independent likelihood and composite likelihood methods outperform stochastic partial differential equation approximations in terms of computation time and returning parameter estimates that nearly maximize the likelihood.

The leverage scores of a full-column rank matrix A are the squared row norms of any orthonormal basis for range(A). We show that corresponding leverage scores of two matrices A and A + \Delta A are close in the relative sense, if they have large magnitude and if all principal angles between the column spaces of A and A + \Delta A are small. We also show three classes of bounds that are based on perturbation results of QR decompositions. They demonstrate that relative differences between individual leverage scores strongly depend on the particular type of perturbation \Delta A. The bounds imply that the relative accuracy of an individual leverage score depends on: its magnitude and the two-norm condition of A, if \Delta A is a general perturbation; the two-norm condition number of A, if \Delta A is a perturbation with the same norm-wise row-scaling as A; (to first order) neither condition number nor leverage score magnitude, if \Delta A is a component-wise row-scaled perturbation. Numerical experiments confirm the qualitative and quantitative accuracy of our bounds., Comment: This version has been accepted to SIMAX but has not yet gone through copy editing

Given a real matrix A with n columns, the problem is to approximate the Gram product AA^T by c << n weighted outer products of columns of A. Necessary and sufficient conditions for the exact computation of AA^T (in exact arithmetic) from c >= rank(A) columns depend on the right singular vector matrix of A. For a Monte-Carlo matrix multiplication algorithm by Drineas et al. that samples outer products, we present probabilistic bounds for the 2-norm relative error due to randomization. The bounds depend on the stable rank or the rank of A, but not on the matrix dimensions. Numerical experiments illustrate that the bounds are informative, even for stringent success probabilities and matrices of small dimension. We also derive bounds for the smallest singular value and the condition number of matrices obtained by sampling rows from orthonormal matrices., Comment: Update to title in third version. Major revisions in second version including new bounds and a more detailed experimental section. Submitted to SIMAX

Motivated by the least squares solver Blendenpik, we investigate three strategies for uniform sampling of rows from m x n matrices Q with orthonormal columns. The goal is to determine, with high probability, how many rows are required so that the sampled matrices have full rank and are well-conditioned with respect to inversion. Extensive numerical experiments illustrate that the three sampling strategies (without replacement, with replacement, and Bernoulli sampling) behave almost identically, for small to moderate amounts of sampling. In particular, sampled matrices of full rank tend to have two-norm condition numbers of at most 10. We derive a bound on the condition number of the sampled matrices in terms of the coherence {\mu} of Q. This bound applies to all three different sampling strategies; it implies a, not necessarily tight, lower bound of O(m {\mu} ln(n)) for the number of sampled rows; and it is realistic and informative even for matrices of small dimension and the stringent requirement of a 99 percent success probability. For uniform sampling with replacement we derive a potentially tighter condition number bound in terms of the leverage scores of Q. To obtain a more easily computable version of this bound, in terms of just the largest leverage scores, we first derive a general bound on the two-norm of diagonally scaled matrices. To facilitate the numerical experiments and test the tightness of the bounds, we present algorithms to generate matrices with user-specified coherence and leverage scores. These algorithms, the three sampling strategies, and a large variety of condition number bounds are implemented in the Matlab toolbox kappa_SQ_v3., Comment: Under revision for SIMAX. Matlab code can be downloaded from http://www4.ncsu.edu/~tawentwo/kappaSQ_upToDate.zip

A sequence of approximations for the determinant and its logarithm of a complex matrixis derived, along with relative error bounds. The determinant approximations are derived from expansions of det(X)=exp(trace(log(X))), and they apply to non-Hermitian matrices. Examples illustrate that these determinant approximations are efficient for lattice simulations of finite temperature nuclear matter, and that they use significantly less space than Gaussian elimination. The first approximation in the sequence is a block diagonal approximation; it represents an extension of Fischer's and Hadamard's inequalities to non-Hermitian matrices. In the special case of Hermitian positive-definite matrices, block diagonal approximations can be competitive with sparse inverse approximations. At last, a different representation of sparse inverse approximations is given and it is shown that their accuracy increases as more matrix elements are included., Comment: 14 pages, 6 figures, 3 tables

La Budde's method computes the characteristic polynomial of a real matrix A in two stages: first it applies orthogonal similarity transformations to reduce A to upper Hessenberg form H, and second it computes the characteristic polynomial of H from characteristic polynomials of leading principal submatrices of H. If A is symmetric, then H is symmetric tridiagonal, and La Budde's method simplifies to the Sturm sequence method. If A is diagonal then La Budde's method reduces to the Summation Algorithm, a Horner-like scheme used by the MATLAB function POLY to compute characteristic polynomials from eigenvalues. We present recursions to compute the individual coefficients of the characteristic polynomial in the second stage of La Budde's method, and derive running error bounds for symmetric and nonsymmetric matrices. We also show that La Budde's method can be more accurate than POLY, especially for indefinite and nonsymmetric matrices A. Unlike POLY, La Budde's method is not affected by illconditioning of eigenvalues, requires only real arithmetic, and allows the computation of individual coefficients., Comment: 24 pages, 10 figures

We introduce a new approximation to nucleon matrix determinants that is physically motivated by chiral effective theory. The method involves breaking the lattice into spatial zones and expanding the determinant in powers of the boundary hopping parameter., Comment: 20 pages, 6 figures, revtex4 (version to appear in PRC)