Reproducible Tall-Skinny QR

We explored half-precision implementation of blocked QR algorithms with the following motivations: 1. New GPUs perform fast half-precision arithmetic (4 to 16 times as fast as doubleprecision). 2. QR factorization is a basic linear algebra tool useful for many physics and data analysis applications. 3. Communication-avoiding, parallelizable QR algorithms already exist for talland- skinny matrices. While the standard QR algorithms are highly unstable in half-precision, our numerical simulations show that the Tall-and-Skinny QR (TSQR) algorithm can improve the backward error of QR factorization. When using subspace iteration for graph clustering applications, half-precision accuracy in forming the eigenspace is sufficient for clustering with high precision and recall for some medium-scale benchmark problems. Note that all half-precision arithmetic were simulated with conversions to single-precision floats. Therefore, the results from this work are somewhat optimistic.

1. Linear systems: Introduction to Linear Systems. Gauss Elimination. Numerical Solutions. 2. Vectors: Vector Operations. Dot Product. Span. Linear Independence. The Product [Ax]. Cross Product. Lines, Planes, and Hyperplanes. 3. Matrices: Matrix Operations. Matrix Inverse. Elementary and Invertible Matrices. LU Factorization. 4. Vector spaces: Subspaces of [R to the n Power]. Vector Spaces. Linear Independence Bases. Dimension. Coordinate Vectors and Change of Basis. Rank and Nullity. Applications to Coding Theory. 5. Linear transformations: Matrix Transformations. Linear Transformations. Kernel and Range. The Matrix of a Linear Transformation. The Algebra of Linear Transformations. 6. Determinants: Determinants Cofactor Expansion. Properties of Determinants. The Adjoint Cramer's Rule. Determinants with Permutations. 7. Eigenvalues and eigenvectors: Eigenvalues and Eigenvectors. Diagonalization. Approximations of Eigenvalues and Eigenvectors. Applications to Dynamical Systems. Applications to Markov Chains. 8. Dot and inner products: Orthogonal Sets and Matrices. Orthogonal Projections Gram-Schmidt Process. The QR Factorization. Least Squares. Orthogonalization of Symmetric Matrices. Quadratic Forms and Conic Sections. The Singular Value Decomposition (SVD). Inner Products. Appendices: Answers to selected exercises Each chapter concludes with Applications, Mini-Projects, and Computer Exercises.

The AllReduce algorithm is a promising new algorithm for parallelizing the Householder QR decomposition A = QR of a tall and skinny matrix. It divides the input matrix A vertically in a recursive manner, computes the QR decompositions of each submatrix independently, and merges the results to obtain the QR decomposition of A. While this algorithm has been shown to achieve excellent speedup in various parallel environments, its rounding error properties have not been elucidated yet. In this paper, we present theoretical error analysis of the AllReduce algorithm. Specifically, we derive bounds for the backward error of A and deviation from orthogonality of the computed Q factor. Our analysis shows that both of these bounds are smaller than their counterparts for the conventional Householder QR algorithm. Moreover, the bounds decrease as the number of submatrices increases. These results are supported by numerical experiments. Thus we can conclude that the AllReduce algorithm can be used as a reliable method of orthogonalization in parallel environments.

Chapter 18 - The Algebraic Eigenvalue Problem

A multiple-input-multiple-output (MIMO) system using Space-Time Block Coding (STBC) techniques can be implemented to improve the capacity of a wireless link. The optimal decoder used the maximum likelihood(ML). But as the number of the antennas in the system and the data rates increase, the ML detector becomes too complex to use. The detectors the least complex and efficient are Zero-Forcing(ZF), MMSE and especially the Sphere Decoding(SD) that is based on QR Decomposition(QRD). In this paper, we study the Gram-Schmidt, Householder and Givens Rotation decomposition methods and compare their efficiency in terms of computational complexity and error rate performance. The channel coding is added to the chain of transmission in order to improve performance. SD receivers are optimized for coded MIMO. Simulations show that the SD detector based on Gram-Schmidt decomposition with channel encoding gives the best results.

Multishift QR algorithms are efficient for solving the symmetric tridiagonal eigenvalue problem on a parallel computer. In this paper, we focus on three variants of the multishift QR algorithm, namely, the conventional multishift QR algorithm, the deferred shift QR algorithm and the fully pipelined multishift QR algorithm, and construct performance models for them. Our models are designed for shared-memory parallel machines, and given the basic performance characteristics of the target machine and the problem size, predict the execution time of these algorithms. Experimental results show that our models can predict the relative performance of these algorithms to the accuracy of 10% in many cases. Thus our models are useful for choosing the best algorithm to solve a given problem in a specified computational environment, as well as for finding the best value of the performance parameters.

Chapter 17 - Implementing the QR Decomposition

SVD and QR Decomposition have attracted much attention recently for solving linear algebra and digital signal processing problems. The House-holder QR Decomposition requires a smaller number of operations then the Givens QR Decomposition and therefore is strongly recommended for sequential computer implementations (as found in nearly all known libraries like NAG, LINPACK, EISPACK etc.) However for parallel implementation it has not been widely used, because the Givens Rotations possess inherent parallelism and the Householder Reflections do not. This is true if one analyzes the original algorithm. Four independent loops and three bottlenecks between the loops are constraint for pipelining the computations. This leads to an inefficient solution, i.e. 4n time moments for each iteration and since there are n iterations then the total time to execute the algorithm is O(4n 2 ). Compared to 3n m 2 time steps for the Givens QR Decomposition it is uneconomic. In this paper we have used some recent results for the elimination of the computational and data broadcast, and data synchronization to derive a fully localized form of the Householder QR Decomposition algorithm. We have succeeded in reorganizing and transforming the algorithm from three bottlenecks to one and four loops to one. A linear and a double pipeline array of n+1 processors are presented to solve the problem in O ( n 2 /2) time steps. It is also shown that the bottlenecks for the bidiagonalization and SVD computation cannot be eliminated.

QR and LU decompositions are the most important matrix decomposition algorithms. Many studies work on accelerating these algorithms by FPGA or ASIC in a case by case style. In this paper, we propose a unified framework for the matrix decomposition algorithms, combining three QR decomposition algorithms and LU algorithm with pivoting into a unified linear array structure. The QR and LU decomposition algorithms exhibit the same two-level loop structure and the same data dependency. Utilizing the similarities in loop structure and data dependency of matrix decomposition, we unify a fine-grained algorithm for all four matrix decomposition algorithms. Furthermore, we present a unified coprocessor structure with a scalable linear array of processing elements (PEs), in which four types of PEs are same in the structure of memory channels and PE connections, but the only difference exists in the internal structure of data path. Our unified co-processor, which is IEEE 32-bit floating-point precision, is implemented and mapped onto a Xilinx Virtex5 FPGA chip. Experimental results show that our co-processors can achieve speedup of 2.3 to 14.9 factors compared to a Pentium Dual CPU with double SSE threads.

Out-of-core processing aims to handle large amount of data when the memory is limited. There exists several out-of-core applications including disk-memory and CPU-GPU processing. Ideally, these out-of-core applications can be expected to be close to the peak performance of the in-core computations, if the data movement between different memory hierarchies can be overlapped by the in-core computations effectively. However, with the emergence of matrix accelerators such as TensorCore GPU, the imbalance between the speed of computations and data movement is further exacerbated, such that even high computation intensity kernels can be dominated by data movement cost. In such cases, the algorithms need to be redesigned to reduce communication volume and overlap the data movement by pipelines. In this paper, we select classic Gram-Schmidt QR factorization as an example to illustrate our recursive strategy, which shows smaller amount of data movement and higher overlapping ratio than the conventional blocking QR factorization algorithm. The results suggest this technique can potentially be applied to broader matrix computations kernels.

Compared with the single tree detector, the layered orthogonal lattice detector (LORD), developed by Siti et al. , is a well-known soft-output multiple input multiple output detector to exploit ${M}$ parallel tree traversals to deliver data with ${M}$ times of detection throughput rate. The preprocessing QR-decomposition (QRD) of the ${M}$ -by- ${M}$ channel matrix for the single tree detector is of complexity proportional to ${M^{3}}$ . However, the preprocessing QRD for the LORD needs to compute the ${M}$ permuted channel matrices that are constructed from the original ${M}$ -by- ${M}$ channel matrix through the root conditioning criterion. The original LORD algorithm for this root conditioning QRD (RC-QRD) relies on the Gram-Schmidt orthogonalization and is of complexity proportional to ${M^{4}}$ for large ${M}$ . In this brief, we apply the Givens rotation and take advantage of the relationships among the ${M}$ permuted matrices to develop an RC-QRD algorithm with complexity proportional to ${M^{3}}$ . Furthermore, when ${M}$ is large, our proposed RC-QRD algorithm requires the number of real Givens rotations about 1.8 times necessary for computing a conventional matrix QRD. Also, for ${M=4}$ , our proposed RC-QRD hardware architecture requires gate count 2.1 times that required by the conventional triangular systolic array to compute a matrix QRD. Accordingly, with only about two times of complexity for the preprocessing RC-QRD, the LORD is able to perform ${M}$ tree traversals to deliver data with ${M}$ times of throughput rate.

Signal processing in encrypted domain has become an important mean to protect privacy in an untrusted network environment. Due to the limitations of the underlying encryption methods, many useful algorithms that are sophisticated are not well implemented. Considering that QR decomposition is widely used in many fields, in this paper, we propose to implement QR decomposition in homomorphic encrypted domain. We firstly realize some necessary primitive operations in homomorphic encrypted domain, including division and open square operation. Gram-Schmidt process is then studied in the encrypted domain. We propose the implementation of QR decomposition in the encrypted domain by using the secure implementation of Gram-Schmidt process. We conduct experiments to demonstrate the effectiveness and analyze the performance of the proposed outsourced QR decomposition.

Assume that two subspaces F and G of a unitary space are defined as the ranges (or null spaces) of given rectangular matrices A and B . Accurate numerical methods are developed for computing the principal angles θ k ( F , G ) {\theta _k}(F,G) and orthogonal sets of principal vectors u k ∈ F {u_k} \in F and v k ∈ G , k = 1 , 2 , ⋯ , q = dim ⁡ ( G ) ≦ dim ⁡ ( F ) {v_k} \in G,k = 1,2, \cdots ,q = \dim (G) \leqq \dim (F) . An important application in statistics is computing the canonical correlations σ k = cos ⁡ θ k {\sigma _k} = \cos {\theta _k} between two sets of variates. A perturbation analysis shows that the condition number for θ k {\theta _k} essentially is max ( κ ( A ) , κ ( B ) ) \max (\kappa (A),\kappa (B)) , where κ \kappa denotes the condition number of a matrix. The algorithms are based on a preliminary QR -factorization of A and B (or A H {A^H} and B H {B^H} ), for which either the method of Householder transformations (HT) or the modified Gram-Schmidt method (MGS) is used. Then cos θ k \cos \;{\theta _k} and sin θ k \sin \;{\theta _k} are computed as the singular values of certain related matrices. Experimental results are given, which indicates that MGS gives θ k {\theta _k} with equal precision and fewer arithmetic operations than HT. However, HT gives principal vectors, which are orthogonal to working accuracy, which is not generally true for MGS. Finally, the case when A and/or B are rank deficient is discussed.

This paper presents algorithms for quintic time scaling trajectory generation of an end- effector, while discretizing a Bézier curve satisfying a finite jerk constraint. The well-known Gaussian elimination and three versions of QR algorithms, i.e., classical, Givens, and House- holder, are utilized to identify a smooth trajectory of an end-effector. Based on these algorithms, a least-time complexity algorithm will be identified for the quintic time-scaling trajectory generation. Imposing a finite jerk constraint reduces the vibrations within a servicing robotic arm or any other manipulator. Moreover, a Bézier curve is chosen to define the path that an end-effector should follow. Subsequently, the path is discretized and constraints are imposed. The coefficients of time-scaling for each discretized segment are calculated by solving systems of linear equations resembling a confluent Vandermonde matrix. Since there are multiple time segments, we iteratively use Gaussian elimination, Gram-Schmidt process, Householder reflections, and Givens rotations to solve a multiple system of equations and determine the lowest time-scaling, opposed to conventional approaches, to determine a smooth quintic time scaling for end-effector trajectory generation of a servicing robotic arm.

The QR factorisation is a cornerstone of numerical linear algebra, essential for solving overdetermined linear systems, eigenvalue problems, and various scientific computing tasks. However, computing it for ill-conditioned tall-and-skinny (TS) matrices on large-scale distributed-memory systems, particularly those with multiple GPUs, presents significant challenges in balancing numerical stability, high performance, and efficient communication. Traditional Householder-based QR methods provide numerical stability but perform poorly on TS matrices due to their reliance on memory-bound kernels. This paper introduces a novel algorithm for computing the QR factorisation of ill-conditioned TS matrices based on CholeskyQR methods. Although CholeskyQR is fast, it typically fails due to severe loss of orthogonality for ill-conditioned inputs. To solve this, our new algorithm, mCQRGSI+, combines the speed of CholeskyQR with stabilising techniques from the Gram–Schmidt process. It is specifically optimised for distributed multi-GPU systems, using adaptive strategies to balance computation and communication. Our analysis shows the method achieves accuracy comparable to Householder QR, even for extremely ill-conditioned matrices (condition numbers up to 1016). Scaling experiments demonstrate speedups of up to 12× over ScaLAPACK and 16× over SLATE’s CholeskyQR2. This work delivers a method that is both robust and highly parallel, advancing the state-of-the-art for this challenging class of problems.