Conditions for Digit Stability in Iterative Methods Using the Redundant Number Representation

Iterative methods play an important role in science and engineering applications, with uses ranging from linear system solvers in finite element methods to optimization solvers in model predictive control. Recently, a new computational strategy for iterative methods called ARCHITECT was proposed by Li et al. in [1] that uses the redundant number representation to create "stable digits" in the Most-significant Digits (MSDs) of an iterate, allowing the future iterations to assume the stable MSDs have not changed their value, eliminating the need to recompute them. In this work, we present a theoretical analysis of how these "stable digits" arise in iterative methods by showing that a Fejér monotone sequence in the redundant number representation can develop stable MSDs in the elements of the sequence as the sequence grows in length. This property of Fejér monotone sequences allows us to expand the class of iterative methods known to have MSD stability when using the redundant number representation to include any fixed-point iteration of a contractive Lipschitz continuous function. We then show that this allows for the theoretical guarantee of digit stability not just in the Jacobi method that was previously analyzed by Li et al. in [2], but also in other commonly used methods such as Newton's method

Matrix Multiplication in Multiword Arithmetic: Error Analysis and Application to GPU Tensor Cores

In multiword arithmetic, a matrix is represented as the unevaluated sum of two or more lower-precision matrices, and a matrix product is formed by multiplying the constituents in low precision. We investigate the use of multiword arithmetic for improving the performance-accuracy tradeoff of matrix multiplication with mixed precision block fused multiply-add (FMA) hardware, focusing especially on the tensor cores available on NVIDIA GPUs. Building on a general block FMA framework, we develop a comprehensive error analysis of multiword matrix multiplication. After confirming the theoretical error bounds experimentally by simulating low precision in software, we use the cuBLAS and CUTLASS libraries to implement a number of matrix multiplication algorithms using double-fp16 (double-binary16) arithmetic. When running the algorithms on NVIDIA V100 and A100 GPUs, we find that double-fp16 is not as accurate as fp32 (binary32) arithmetic despite satisfying the same worst-case error bound. Using probabilistic error analysis, we explain why this issue is likely to be caused by the rounding mode used by the NVIDIA tensor cores, and propose a parameterized blocked summation algorithm that alleviates the problem and significantly improves the performance-accuracy tradeoff

A Parametrization of Structure-Preserving Transformations for Matrix Polynomials

Given a matrix polynomial $A(\lambda)$ of degree $d$ and the associated vector space of pencils \DLP(A) described in Mackey, Mackey, Mehl, and Mehrmann [SIAM J. Matrix Anal. Appl., 28 (2006), pp. 971-1004], we construct a parametrization for the set of left and right transformations that preserve the block structure of such pencils, and hence produce a new matrix polynomial \At(\lambda) that is still of degree $d$ and is unimodularly equivalent to $A(\lambda)$. We refer to such left and right transformations as structure-preserving transformations (SPTs). Unlike previous work on SPTs, we do not require the leading matrix coefficient of $A(\lambda)$ to be nonsingular. We show that additional constraints on the parametrization lead to SPTs that also preserve extra structures in $A(\lambda)$ such as symmetric, alternating, and $T$-palindromic structures. Our parametrization allows easy construction of SPTs that are low-rank modifications of the identity matrix. The latter transform $A(\lambda)$ into an equivalent matrix polynomial \At(\lambda) whose $j$th matrix coefficient \At_j is a low-rank modification of $A_j$. We expect such SPTs to be one of the key tools for developing algorithms that reduce a matrix polynomial to Hessenberg form or tridiagonal form in a finite number of steps and without the use of a linearization

Probabilistic Rounding Error Analysis of Householder QR Factorization

When an $m\times n$ matrix is premultiplied by a product of $n$ Householder matrices the worst-case normwise rounding error bound is proportional to $mnu$, where $u$ is the unit roundoff. We prove that this bound can be replaced by one proportional to $\sqrt{mn}u$ that holds with high probability if the rounding errors are mean independent and of mean zero. The proof makes use of a matrix concentration inequality. In particular, this result applies to Householder QR factorization. The same square rooting of the error constant applies to two-sided transformations by Householder matrices and hence to standard QR-type algorithms for computing eigenvalues and singular values. It also applies to Givens QR factorization. These results complement recent probabilistic rounding error analysis results for inner-product based algorithms and show that the square rooting effect is widespread in numerical linear algebra. Our numerical experiments, which make use of a new backward error formula for QR factorization, show that the probabilistic bounds give a much better indicator of the actual backward errors and their rate of growth than the worst-case bounds

Computing the square root of a low-rank perturbation of the scaled identity matrix

We consider the problem of computing the square root of a perturbation of the scaled identity matrix, A = α Iₙ + UVᴴ, where U and V are n × k matrices with k ≤ n. This problem arises in various applications, including computer vision and optimization methods for machine learning. We derive a new formula for the pth root of A that involves a weighted sum of powers of the pth root of the k × k matrix α Iₖ + VᴴU. This formula is particularly attractive for the square root, since the sum has just one term when p = 2. We also derive a new class of Newton iterations for computing the square root that exploit the low-rank structure. We test these new methods on random matrices and on positive definite matrices arising in applications. Numerical experiments show that the new approaches can yield much smaller residual than existing alternatives and can be significantly faster when the perturbation UVᴴ has low rank

An approach to protein structure using information geometry

In the light of recent structural developments in DNA structural diversity crystallographic studies and the Protein Data Bank*, this note is intended to draw attention to an interesting feature of the ordering of amino acids along protein chains. They all exhibited clustering compared to a random distribution, so there is a stable long range ordering that is unexpected. To date we have no clear explanation of why this should be the case. * https://doi.org/10.1016/j.jbc.2021.100553

Randomized Low Rank Matrix Approximation: Rounding Error Analysis and a Mixed Precision Algorithm

The available error bounds for randomized algorithms for computing a low rank approximation to a matrix assume exact arithmetic. Rounding errors potentially dominate the approximation error, though, especially when the algorithms are run in low precision arithmetic. We give a rounding error analysis of the method that computes a randomized rangefinder and then computes an approximate singular value decomposition approximation. Our analysis covers the basic method and the power iteration for the fixed-rank problem, as well as the power iteration for the fixed-precision problem. We see that for the fixed-rank problem, the bound for the power iteration is favourable in terms of simplicity and rounding error contribution. We give both worst-case and probabilistic rounding error bounds as functions of the problem dimensions and the rank. The worst-case bounds are pessimistic, but the probabilistic bounds are reasonably tight and still reliably bound the error in practice. We also propose a mixed precision version of the algorithm that offers potential speedups by gradually decreasing the precision during the execution of the algorithm

Randomized Low Rank Matrix Approximation: Rounding Error Analysis and a Mixed Precision Algorithm

The available error bounds for randomized algorithms for computing a low rank approximation to a matrix assume exact arithmetic. Rounding errors potentially dominate the approximation error, though, especially when the algorithms are run in low precision arithmetic. We give a rounding error analysis of the method that computes a randomized rangefinder and then computes an approximate singular value decomposition approximation. Our analysis covers the basic method and the power iteration for the fixed-rank problem, as well as the power iteration for the fixed-precision problem. We give both worst-case and probabilistic rounding error bounds as functions of the problem dimensions and the rank. The worst-case bounds are pessimistic, but the probabilistic bounds are reasonably tight and still reliably bound the error in practice. We also propose a mixed precision version of the algorithm that offers potential speedups by gradually decreasing the precision during the execution of the algorithm

Quasi-triangularization of matrix polynomials over arbitrary fields

In \cite{TasTisZab}, Taslaman, Tisseur, and Zaballa show that any regular matrix polynomial P(\la) over an algebraically closed field is spectrally equivalent to a triangular matrix polynomial of the same degree. When P(\la) is real and regular, they also show that there is a real quasi-triangular matrix polynomial of the same degree that is spectrally equivalent to P(\la), in which the diagonal blocks are of size at most $2 \times 2$. This paper generalizes these results to regular matrix polynomials P(\la) over arbitrary fields \bF, showing that any such P(\la) can be quasi-triangularized to a spectrally equivalent matrix polynomial over \bF of the same degree, in which the largest diagonal block size is bounded by the highest degree appearing among all of the \bF-irreducible factors in the Smith form for P(\la)