# QR factorization of ill-conditioned tall-and-skinny matrices on distributed-
memory systems

Nenad Mijić, Abhiram Kaushik, and Davor Davidović N. Mijić, A. Kaushik and D.
Davidović was with the Centre for Informatics and Computing, Ruđer Bošković
Institute, Zagreb, Croatia.
E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In this paper we present a novel algorithm developed for computing the QR
factorisation of extremely ill-conditioned tall-and-skinny matrices on
distributed memory systems. The algorithm is based on the communication-
avoiding CholeskyQR2 algorithm and its block Gram-Schmidt variant. The latter
improves the numerical stability of the CholeskyQR2 algorithm and
significantly reduces the loss of orthogonality even for matrices with
condition numbers up to $10^{15}$. Currently, there is no distributed GPU
version of this algorithm available in the literature which prevents the
application of this method to very large matrices. In our work we provide a
distributed implementation of this algorithm and also introduce a modified
version that improves the performance, especially in the case of extremely
ill-conditioned matrices. The main innovation of our approach lies in the
interleaving of the CholeskyQR steps with the Gram-Schmidt orthogonalisation,
which ensures that update steps are performed with fully orthogonalised
panels. The obtained orthogonality and numerical stability of our modified
algorithm is equivalent to CholeskyQR2 with Gram-Schmidt and other state-of-
the-art methods. Weak scaling tests performed with our test matrices show
significant performance improvements. In particular, our algorithm outperforms
state-of-the-art Householder-based QR factorisation algorithms available in
ScaLAPACK by a factor of $6$ on CPU-only systems and up to $80\times$ on GPU-
based systems with distributed memory.

###### Index Terms:

QR factorisation, CholeskyQR, Gram-Schmidt, Tall-and-skinny matrices, Ill-
conditioned matrices, Graphic processing units, Distributed systems, Parallel
algorithms

##  1 Introduction

In this paper, we develop and analyse scalable algorithms tailored to large-
scale distributed-memory systems that can compute the QR decomposition [1] of
a rectangular matrix $A\in\mathbb{R}^{m\times n}$ with $m\geq n$:

$A=QR,$ (1)

where $Q\in R^{m\times n}$ is an orthogonal matrix and
$R\in\mathbb{R}^{n\times n}$ is an upper triangular matrix. Finding the QR
decomposition of large rectangular matrices is a crucial step in various
numerical methods such as block methods (e.g. in the solution of linear
systems with multiple right-hand sides [2]), the Krylov subspace methods [3],
eigenvalue solvers (e.g. in the reduction to the band form of the multi-stage
eigensolver [4]), subspace iteration [5] and solving a dense least squares
problem of an overdetermined system. In some extreme cases, where $m\gg n$,
i.e. the matrix has many more rows than columns (so-called tall and skinny
matrix - TS), the calculation of the QR factorisation becomes a critical path
and requires special methods. One example is a subspace projection iterative
eigensolver [5] to calculate a small fraction of the extreme eigenvalues.

The most stable and accurate solution to compute QR factorisation for general
matrices is based on Householder reflectors [1], usually referred to as
Householder QR. The algorithm finds a set of orthogonal Householder matrices
that annihilate the entries below the main diagonal of $A$ column by column,
while the rest of the matrix is updated (trailing update). This type of
algorithm, where $R$ is constructed by applying an orthogonal matrix from the
left, $Q^{T}A=R$, is called orthogonal triangularisation and provides good
numerical stability. The highly efficient implementations can be found in many
numerical libraries such as ScaLAPACK [6] for distributed memory systems or
MAGMA [7] for heterogeneous and accelerator-based shared memory systems.
However, the Householder QR cannot achieve high performance for tall and
skinny matrices [8] because the panel factorisation of columns with many more
rows is performed with much slower Level-1 and Level-2 BLAS kernels, which
cannot be compensated by the highly parallel and optimised Level-3 BLAS
kernels used in the trailing update [9].

To avoid working with very tall and skinny panels, a Tall Skinny QR (TSQR)
algorithm [8, 10] was developed, which offers better parallelisation on
systems with distributed memory and reduces communication. The main idea is to
split the input matrix into a cyclic 1-D block row layout so that QR
factorisations can be performed on local blocks concurrently. In the following
reduction steps, the intermediate R-factors are grouped into pairs and
orthogonalised. This step is repeated until the final upper triangular $R$ is
reached. Although the degree of parallelism is significantly increased
compared to traditional Householder QR, a large number of flops is still
performed in terms of Level-1 and Level-2 BLAS kernels when computing the QR
of the local blocks. The distributed version of TSQR is implemented in the
SLATE library [11] within the Communication Avoiding QR (CAQR) algorithm,
where it is used for QR factorisation of tall-and-skinny panels. CAQR was
developed for general matrices and not for TS matrices. It is more performant
than the standard distributed QR as it replaces the communication volume with
more flops.

An alternative to the Householder-based QR algorithms for tall and skinny
matrices is the CholeskyQR algorithm [1]. CholeskyQR is a simple algorithm
with very low communication overhead and about half the arithmetic cost of
TSQR. The algorithm is based on Cholesky factorisation and triangular
orthogonalisation ($Q=AR^{-1}$) and is generally faster than TSQR. The main
drawback is that the algorithm suffers from a loss of orthogonality and is
numerically unstable. To compensate for this shortcoming, the authors in [12,
13] have proposed a new algorithm called CholeskyQR2 (see section 4 for more
details). The main idea is to repeat the CholeskyQR algorithm twice to
increase orthogonality. The new algorithm achieves accuracy and orthogonality
comparable to TSQR, but also requires twice as much communication as
CholeskyQR, with arithmetic cost equivalent to those of TSQR. However, the
algorithm fails for input matrices with very high condition numbers
($\kappa(A)>u^{-1/2}$, where $u$ is the unit roundoff). Distributed parallel
CholeskyQR is available in the SLATE library, where it is used as one of the
preferred methods for solving least squares and QR factorisation problems.
However, there is no report on how numerical instabilities have been tackled
with ill-conditioned matrices.

The CholeskyQR2 has proven to be the method of choice for solving least
squares and eigenvalue problems on highly parallel computer clusters. In [14],
the CholeskyQR2 is extended over a 3D processor grid, which results in less
communication between processors and outperforms the QR factorisation of
ScaLAPACK by up to $3.3\times$ for strong scaling tests and $1.9\times$ for
weak scaling tests. However, the main issue with the numerical stability of
the algorithm remains. If the condition number of the input matrix is too
large, the constructed Gram matrix fails to be positive semidefinite, causing
the Cholesky factorisation step to fail [12].

This issue is directly addressed in Shifted CholeskyQR algorithm [15]. The
idea is to construct a shifted Gram matrix $\hat{W}:=A^{T}A+\sigma I$ in which
the shift factor $\sigma$ ensures the numerical stability of the Cholesky
factorisation. However, the proposed solution does not work for matrices with
condition numbers close to or greater than $u^{-1/2}$. In the case of
extremely ill-conditioned matrices, the Shifted CholeskyQR is used as a
preconditioner for CholeskyQR2 and the new algorithm is called Shifted
CholeskyQR3. Although this approach significantly improves stability, it leads
to $50\%$ more flops compared to CholeskyQR2.

Another approach to improve the applicability of CholeskyQR for ill-condition
matrices is to use the LU decomposition as preconditioning for the Cholesky
factorisation [16]. After the LU factorisation $PA:=LU$ has been calculated,
the Cholesky decomposition can be applied to the matrix $L^{T}L$, since $L$ is
usually better conditioned than $A$. After the first call of CholeskyQR with
LU, the $Q_{1}$ obtained is refined by the second call of the standard
CholeskyQR algorithm. The achieved orthogonality and the residuals of the
algorithm are comparable to those of the Householder QR even for very ill-
conditioned matrices, i.e. $\kappa(A)>u^{-1/2}$. However, since the algorithm
computes the LU decomposition with partial pivoting, it exhibits less
parallelism compared to CholeskyQR2 and is about $1.5$ times slower on shared
memory systems and between 3 and 5 times on distributed memory systems.

The most recent improvements to the orthogonality and stability of the
CholeskyQR model was proposed in [17] with the so-called Randomised
Householder-Cholesky QR Factorization with Multisketching. The proposed
solution introduces up to two randomised sketch matrices (multisketching),
which prove that the orthogonality error is bounded by a constant of the order
of unit roundoff for matrices of arbitrary condition number. The first step is
to compute the randomised Householder QR that generates the matrix $Q_{1}$
orthogonal to the given sketch matrices. As in other similar approaches, the
obtained matrix is reorthogonalised in the second step by calling CholeskyQR2
to produce a fully orthogonal matrix $Q$. The authors reported that their
approach is applicable to extremely tall-and-skinny matrices ($n\leq 0.01\%m$
) and was negligibly faster than CholeskyQR2, but more stable than the shifted
CholeskyQR3.

Since the orthogonality error of CholeskyQR depends quadratically on the
condition number, the authors in [18] have proposed a mixed-precision
approach. In this approach, the input and output matrices are retained in
their required precision, while certain intermediate results are calculated at
doubled precision. The analysis has shown that the orthogonality error of the
mixed-precision CholeskyQR approach has a linear dependence on the condition
number of the input matrix when doubled precision is used. The main drawback
of the proposed algorithm is that the number of floating point instructions
increases significantly when doubled precision is used (especially in the case
when the target precision is 64-bit).

In [19], the stability of the Cholesky decomposition and orthogonality for
very ill-conditioned matrices is restored by combining the CholeskyQR with the
Gram-Schmidt method of re-orthogonalisation. A detailed overview can be found
in the subsection 5.2. The proposed algorithm has shown that the obtained
orthogonality is equivalent to the Householder-based QR factorisations
(general and TSQR), but cannot exploit massively parallel systems as the
method uses very tall and skinny column panels with full column size.

The original scientific contributions in this paper are:

  * •

A novel algorithm for computing the QR factorisation of ill-conditioned tall-
and-skinny matrices which combines the classical CholeskyQR2 algorithm with
Gram-Schmidt reorthogonalisation.

  * •

Orthogonality and numerical stability of the new algorithm are comparable to
those of Householder-based QR factorisation algorithms.

  * •

QR factorisation optimised for distributed memory GPU systems.

  * •

The algorithm outperforms other methods, including the Shifted CholeskyQR3
algorithm and traditional Householder-based approaches (such as ScaLAPACK’s
PDGEQRF).

  * •

An analysis and synthesis of the state-of-the-art in CholeskyQR-based QR
factorisation algorithms for tall-and-skinny matrices.

When defining the workflow in this paper, our main goal was to guide the
readers through all stages of the algorithm development, starting from the
basic version and gradually improving it to an advanced version. At each step,
we provide a numerical and a performance analysis which helps the reader to
better understand the improvements over the previous version. We believe that
this approach provides a deeper understanding of our algorithm. The paper
starts with the testing environment, state-of-the-art algorithms and test
matrices described in section 2. The CholeskyQR algorithm, a basic building
block of our new algorithm, is described in section 3, followed by the
analysis of the CholeskyQR2 algorithm in section 4. The main part of the paper
is in section 5, where we describe step by step how the algorithm evolves from
the first improvement, shifted CholeskyQR3, through the introduction of the
Gram-Schmidt method to our final version, the modified CholeskyQR2 with Gram-
Schmidt. The scalability analysis of the new algorithm and the comparison with
ScaLAPACK are shown and discussed in section 6. The paper is concluded in
section 7 with a summary of what has been achieved and an overview of future
work.

##  2 Testing environment

###  2.1 Testing platform

The numerical stability and performance were tested on the Supek supercomputer
at the University of Zagreb, University Computing Centre (SRCE). The tests
were performed on two partitions: GPU and CPU. The GPU partition consists of
20 nodes connected to the Cray Slingshot interconnect. Each node has an AMD
EPYC 7763 CPU with 64 cores and 512 GB of main memory, supported by four
NVIDIA A100 Tensor Core GPUs with 40 GB of device memory each. The CPU
partition consists of 52 nodes similar to the GPU partition, but with 2 AMD
EPYC processors (128 cores in total) and without GPU accelerators. The codes
were compiled with the GCC 12.1.0 compiler and the CUDA 11.6 library, which
provides BLAS and LAPACK functionality. Cray MPICH 8.1.20 and NVIDIA NCCL
2.12.12 were used for communication. All experiments were performed in double
precision arithmetic and an average execution time of 10 runs is reported for
each experiment.

The numerical accuracy is obtained by analysing the orthogonality of the
obtained matrix $Q$ using the formula $||Q^{T}Q-I||_{F}/\sqrt{n}$, where $I$
is the identity matrix, and the residual $||QR-A||_{F}/||A||_{F}$. Both the
orthogonality and the residual should be of the order of $O(u)$, where $u$ is
a machine precision of a certain numerical type precision (e.g. double
precision).

###  2.2 Test matrix suite

The input data are artificially generated matrices with a condition number
$\kappa(A)$ in the range of $\\{10^{0},10^{1},\ldots,10^{15}\\}$. The matrices
are generated using the SVD ($U\hat{\Sigma}V$), where U and V are left and
right singular vector matrices obtained from the SVD of a random input matrix.
The new diagonal matrix $\Sigma$ is constructed so that the diagonal elements
are $(1,\sigma^{\frac{1}{n-1}},\dots,\sigma^{\frac{n-2}{n-1}},\sigma)$, where
the parameter $\sigma$ controls the condition number of the generated matrix
$\kappa(A)\approx\sigma$.

For the numerical stability tests, we used matrices of size $30000\times 3000$
and the condition number in the range $\\{1,10^{1},\ldots,10^{15}\\}$ in
double precision. The strong tests were performed for matrices with $120000$
rows and the number of columns equal to $1\%$, $5\%$ and $10\%$ of the number
of rows. The weak scalability analysis was performed on matrices with
$40k,80k,120k,\ldots,480k$ rows and the number of columns fixed to $3000$,
resulting in blocks of size $10k\times 3k$ per process (MPI or NCCL rank).

###  2.3 Software

In our tests and analyses, we compare the performance of our novel
CholeskyQR2-based algorithms and evaluate their scalability on systems with
distributed memory. To the best of our knowledge, ScaLAPACK is the only
publicly available library that can handle QR factorization on distributed
memory architectures. The downside is that it does not support execution on
GPUs and only implements the Householder-based methods. The only alternative,
the SLATE library, which implements CholeskyQR for distributed multi-GPU
architectures, exhibits the same numerical instabilities for matrices with
very high condition numbers (greater than $10^{8}$), even when CholeskyQR2 is
used. Since it is not possible to compare our solutions with SLATE for ill-
conditioned matrices, the algorithm is only compared with ScaLAPACK.

For simplicity and easier comparison with our solution, the total computation
and communication costs for ScaLAPACK PSGEQRF implementation are:
$2\frac{mn^{2}}{P}-\frac{2}{3}\frac{n^{3}}{P}$ (number of operations),
$\frac{n^{2}}{2}\log P$ (number of words transmitted), $2n\log P$ (number of
messages), where $m$ and $n$ are the number of rows and columns, respectively,
and $P$ is the number of processes.

##  3 CholeskyQR

CholeskyQR is a simple algorithm that calculates the QR factorisation of a
tall and skinny matrix. The pseudocode is described in Algorithm 1. The
algorithm starts with the construction of the Gram matrix $A^{T}A$ (line 1)
using matrix-matrix multiplication, after which the upper triangular matrix
$R$ is obtained via Cholesky factorization (line 2). Finally, the orthogonal
matrix $Q$ is constructed by right-multiplying $A$ with the upper triangular
matrix $R$ (line 3), e.g. using the routine triangular system solve (trsm)
from the LAPACK library. The main advantage of CholeskyQR over Householder-
based algorithms such as TSQR is that all steps in Algorithm 1 can be
implemented as Level-3 BLAS operations, which ensures significantly better
performance on large parallel systems.

Algorithm 1 CholeskyQR

1:$A\in\mathbb{R}^{m\times n}$

2:$Q\in\mathbb{R}^{m\times n}$ orthogonal and $R\in\mathbb{R}^{n\times n}$
upper triangular matrix

3:$W:=A^{T}A$ $\triangleright$ Construct Gram matrix

4:$W=R^{T}R$ $\triangleright$ Cholesky factorization

5:$Q:=AR^{-1}$

Parallelization. Since the CholeskyQR is intended for factorization of
matrices with many more rows than columns, parallelization is done via the row
dimension. The matrix $A$ is divided into one-dimensional row blocks (see Fig.
2 middle ), with each processor $i$ processing a row block $A_{i}$. The Gram
matrix is first constructed locally, with each processor calculating its local
matrix $W_{i}=A^{T}_{i}A_{i}$ ( lines 1–3). Then, a collective communication
is required to collect and sum all local components to obtain the final Gram
matrix $W=W_{1}+W_{2}+\ldots+W_{P}$ (line 4). The Cholesky factorization (line
5) is performed redundantly by each processor using a highly optimized
implementation designed for shared memory systems, e.g. from the MAGMA or
LAPACK library. Finally, the construction of the orthogonal matrix $Q$ (lines
6–7), can be done in parallel as each processor constructs its own part of
$Q$, so that no further communication is required. Overall, CholeskyQR
requires only a single collective communication.

If the matrix $Q$ is required in a single memory space (e.g. a compute node),
an additional global communication (e.g. MPI_gather) is required to collect
the blocks $Q_{i}$ distributed across the MPI ranks. In the remainder of the
paper, we assume that $A$ is already distributed across the processors and $Q$
does not need to be collected on a single MPI rank. This simulates the real-
world scenario where QR factorization is only a part of a larger computer code
where $A$ is already distributed and $Q$ needs to be distributed for further
processing (see [5]). Throughout the paper, the parallel version of CholeskyQR
(CQR) is used as the main building block for all other algorithms.

Algorithm 2 Parallel CholeskyQR (CQR)

1:$A\in\mathbb{R}^{m\times n}$, $P$ number of processors

2:$Q\in\mathbb{R}^{m\times n}$ orthogonal and $R\in\mathbb{R}^{n\times n}$
upper triangular matrix

3:for $i\leftarrow 1,P$ do $\triangleright$ Parallel for

4: $W_{i}:=A_{i}^{T}A_{i}$ $\triangleright$ Compute Gram locally

5:end for

6:$W:=\sum_{i}W_{i}$ $\triangleright$ Allreduce

7:$W=R^{T}R$ $\triangleright$ Cholesky factorization

8:for $i\leftarrow 1,P$ do $\triangleright$ Parallel for

9: $Q_{i}:=A_{i}R^{-1}$ $\triangleright$ Compute Q locally

10:end for

The computational cost (flops) of parallel CholeskyQR is
$\frac{1}{3}n^{3}+2\frac{m}{P}n^{2}+n^{2}\log_{2}P$ with the dominant factor
$\left(2\frac{m}{P}n^{2}\right)$, coming from the construction of the Gram
matrix (gemm/syrk) and the construction of $Q$ (trsm). The total cost of CQR
is about half of the flops required by TSQR. As described in Algorithm 2,
CholeskyQR can be implemented with only one call to a collective communication
routine, for example Allreduce (Algorithm 2, line 4), which makes CQR suitable
for execution on large parallel systems. The number of words transmitted and
the number of messages is $n^{2}\log_{2}P$ and $\log_{2}P$ respectively and
corresponds to the transmission (broadcasting and reduction) of the local
matrices $W_{i}$ and is the same (except for a constant factor) to that of
TSQR. More details on the communication and calculation complexities can be
found in Table I (under CQR).

The biggest disadvantage is that the CholeskyQR is numerically unstable. The
loss of orthogonality of the calculated $Q$ increases with the condition of
the matrix $A$ and is upper bounded by $\mathcal{O}(\epsilon\kappa(A)^{2})$
[20], where $\epsilon$ is the machine precision ($\approx 10^{-16}$). Even for
small condition numbers, e.g. $\kappa(A)=\mathcal{O}(10)$, the deviation from
orthogonality obtained with double precision arithmetic is
$\mathcal{O}(10^{-14})$. Furthermore, since the construction of the gram
matrix squares the condition number, the resulting $W$ is not positive
semidefinite for ill-conditioned matrices, so the Cholesky factorisation step
cannot be performed.

##  4 CholeskyQR2

As already shown, the CholeskyQR cannot always generate orthogonal vectors
($Q$), and as the condition number increases, it even becomes numerically
unstable, which often leads to the failure of the algorithm (reduced stability
of the Cholesky factorisation). To address the problem with orthogonality, a
simple and very effective idea was presented to reorthogonalise the obtained
$Q$ by repeating the CholeskyQR algorithm (see [12]). The proposed algorithm
is called CholeskyQR2.

The pseudocode for the CholeskyQR2 is given in Algorithm 3. The algorithm
starts with the calculation of the QR factorization of $A$ (line 1), for which
the Cholesky QR algorithm (CQR) is used. The result is the matrix $Q_{1}$,
which does not have to be orthogonal. To improve the numerical accuracy of the
obtained matrix, another orthogonalization step is performed on $Q_{1}$ ( line
2). The idea of repeating the orthogonalization was introduced in [21], where
the authors showed that the stability of the Gram-Schmidt algorithm and the
orthogonality of the computed vectors can be improved if the algorithm is
executed twice. The idea is also applicable to the Cholesky QR, since both
algorithms are of the triangular orthogonalization type, i.e. the factor $Q$
is calculated by right multiplication of the triangular matrix $R$ ( Algorithm
1 line 3). Finally, the triangular matrix $R$ is obtained by multiplying the
upper triangular factors (line 3), which are generated by two calls to
CholeskyQR.

Algorithm 3 CholeskyQR2 (CQR2)

1:$A\in\mathbb{R}^{m\times n}$

2:$Q\in\mathbb{R}^{m\times n}$ orthogonal and $R\in\mathbb{R}^{n\times n}$
upper triangular matrix

3:$[Q_{1},R_{1}]:=CQR(A)$

4:$[Q,R_{2}]:=CQR(Q_{1})$

5:$R:=R_{2}R_{1}$

Although repeating CholeskyQR twice significantly improves orthogonality,
CholeskyQR2 is still numerically unstable for ill-conditioned matrices with
condition numbers greater than $\mathcal{O}(u^{-1/2})$, where $\mathbf{u}$ is
the unit roundoff. In practice, this means that the algorithm is stable up to
condition numbers of $10^{8}$ [12, 13] and the deviation from orthogonality of
the computed $Q$ is of order $\mathcal{O}(\kappa(A)^{2}u)$.

TABLE I: Computational (comp) and communication (comm) costs of routines in
parallel CholeskyQR and the CholeskyQR2 algorithms. We assume that the routine
syrk is used to compute the Gram matrix, and not gemm.

algorithm | routine | comp | comm
---|---|---|---
CQR2 | CQR | Gram | $\frac{m}{P}n^{2}$ | -
Gram_reduce | $n^{2}\log_{2}P$ | $n^{2}\log_{2}P$
Cholesky | $\frac{n^{3}}{3}$ | -
Construct_Q | $\frac{m}{P}\ n^{2}$ | -
Total | $\frac{1}{3}n^{3}+2\frac{m}{P}\ n^{2}+n^{2}\log_{2}P$ | $n^{2}\log_{2}P$
Compute_R |  | $\frac{1}{3}\ n^{3}$ | -
| Total |  | $n^{3}+4\frac{m\ n^{2}}{P}+2n^{2}\log_{2}P$ | $2\ n^{2}\log_{2}P$

The parallelization of the algorithm is simple as it exploits the parallelism
within the CholeskyQR algorithm. Therefore, no additional communication is
required, except for two collective routines, one for each call to CholeskyQR.
Since the algorithm repeats CholeskyQR twice, the communication costs (see
Table I CQR2) are twice as high as for CholeskyQR and amount to
$2n^{2}\log_{2}P$. Although the communication costs are much higher than for
TSQR, they are still much lower than for Householder QR. The computational
cost is also twice that of CholeskyQR, plus an additional cost of
$\frac{1}{3}n^{3}$ incurred by the construction of the matrix $R$ (Algorithm 3
line 3). The total computational cost of CholeskyQR2 is comparable to that of
TSQR.

##  5 CholeskyQR variants for extremely ill-conditioned matrices

The source of the numerical instability in the CholeskyQR algorithm is the
Cholesky decomposition, which is sensitive to nearly singular or ill-
conditioned matrices. Such matrices can occur in the construction of the Gram
matrix whose condition number is a squared condition of the original matrix,
resulting in a matrix that is not semi-positive definite in finite precision
arithmetic. In the following subsections, we present recent developments that
address the critical challenge of the breakdown of Cholesky factorization
within CholeskyQR algorithms and introduce our new parallel CholeskyQR-based
algorithm for distributed memory systems. These novel algorithms provide
innovative strategies to circumvent the limitations imposed by ill-conditioned
matrices and ensure the stability and reliability of the Cholesky
factorization process even for extremely ill-conditioned matrices.

###  5.1 Shifted CholeskyQR3

A recent proposal to improve the numerical stability of CholeskyQR is to shift
the Gram matrix in the CholeskyQR routine [15]. The main idea is to decrease
the condition number of the computed Gram matrix thus improving the stability
of the Cholesky factorization. The algorithm is called Shifted Cholesky QR,
and we will refer to it as sCQR throughout the paper. The algorithm presented
in [15] is described in Algorithm 4 and is similar to CholeskyQR except for
the key steps of introducing the shift $s$ in lines 2–3. The shift is chosen
to force the Gram matrix to become positive definite so that the Cholesky
algorithm can be completed. With a well-chosen shift, sCQR is suitable for
matrices of condition number up to $\mathcal{O}(u^{-1})$. In our research, we
choose the conservative approach (as proposed in [22]), where the Frobenius
norm is used instead of norm-2, because of significantly less computational
cost, which ensures numerical stability for our test matrices (see Fig. 1).

Algorithm 4 Shifted Cholesky QR (sCQR)

1:$A\in\mathbb{R}^{m\times n}$

2:$Q\in\mathbb{R}^{m\times n}$ orthogonal and $R\in\mathbb{R}^{n\times n}$
upper triangular matrix

3:$G=A^{T}\ A$

4:$s=\sqrt{m}\,\mathbf{u}\,||A||^{2}_{F}$ $\triangleright$ calculate shift

5:$W=G+sI$ $\triangleright$ construct Gram matrix

6:$W=R^{T}\ R$ $\triangleright$ Cholesky factorization

7:$Q=AR^{-1}$

The numerical stability for extremely ill-conditioned matrices (e.g.
$\kappa(A)\geq 10^{15}$) can be improved by choosing a better shift (usually a
smaller shift of the order of $\mathcal{O}(u\|A\|_{2}$), which can ensure the
stability of the Cholesky factorization. However, experimenting with other
approaches to compute the shift is not the focus of this article as it does
not affect the overall computational complexity or execution time. How to
choose an optimal shift for a given matrix can be read in [15].

The condition number of the computed matrix $Q$ in Algorithm 4 (line 5) in the
shifted CholeskyQR is roughly upper bounded by $\mathcal{O}(u^{-1/2})$, which
is a condition number that ensures the numerical stability of the CholeskyQR2.
To obtain the orthogonality and the residual of order $u$, the authors in [15]
have proposed a solution where the shifted Cholesky QR is used as a
preconditioner for the CholeskyQR2 that ensures its numerical stability. The
combined application of shifted CholeskyQR and CholeskyQR2 constitutes shifted
CholeskyQR3 (sCQR3), which is described in Algorithm 5. The first step (line
1) is to compute the matrix $Q_{1}$ with the reduced condition number, which
is then further orthogonalized via the CholeskyQR2 (line 2). Finally, the
resulting matrix $R$ (line 4) by multiplying intermediate triangular factors
computed by the shifted CholeskyQR and CholeskyQR2.

Algorithm 5 Shifted Cholesky QR3 (sCQR3)

1:$A\in\mathbb{R}^{m\times n}$

2:$Q\in\mathbb{R}^{m\times n}$ orthogonal and $R\in\mathbb{R}^{n\times n}$
upper triangular matrix

3:$[Q_{1},R_{1}]=sCQR(A)$

4:$[Q,R_{2}]=CQR2(Q_{1})$

5:$R:=R_{2}R_{1}$

The parallelization of the shifted CholeskyQR routine is similar to that of
the CholeskyQR routine described in Section 3. The algorithms are the same
except for the computation of the shift factor and shifting of the Gram
matrix. The main additional cost is the calculation of the Frobenius norm of
the matrix in computing the shift (Algorithm 4 line 2). This is done by
summing over the squares of the elements of the partial matrix stored in each
rank and communicating the results to rank 0, where the norm is calculated.
Based on the obtained results, the shift is added to the partial Gram matrix
stored in rank 0. The total computational cost of the parallel shifted
CholeskyQR3 with P processors is:

$\frac{5}{3}n^{3}+6\frac{m\ n^{2}}{P}+3\ n^{2}\ log_{2}P+2\frac{mn}{P}$ (2)

where the last term accounts for the calculation of the Frobenius norm when
computing the shift. In the above expression, we have neglected inexpensive
operations such as the calculation and addition of the shift parameter from
the norm, which cost $\mathcal{O}(n)$ or less. The computational cost of the
shifted CholeskyQR3 is larger than that of CQR2, for an additional cost of
sCQR. Since the cost of sCQR is about the same as that of CQR, the sCQR3 is
about $1.5\times$ higher than CQR2 plus the additional cost ($n^{3}$) required
to compute the final $R$. The communication cost is $50\%$ higher than for
cQR2.

Figure 1: Orthogonality and residuals of sCQR3 and CQR2 as a function of the
condition number, for input matrices with $m=30000$, $n=3000$ and conservative
shift for sCQR3.

With the right choice of shift, the numerical stability of the Cholesky
factorization can be ensured, even for the extremely ill-condition matrices.
However, to compute the orthogonal matrix $Q$, comparable to the Householder-
based approaches, e.g. TSQR, additional computational steps are required
compared to CholeskyQR2. These additional steps significantly increase both
computational cost and total execution time, making the shifted CholeskyQR3
inferior in terms of performance.

###  5.2 CholeskyQR2 with Gram-Schmidt

Another approach to address the problem with instability of the Cholesky
factorization for matrices with a condition number greater than
$\mathcal{O}(u^{-1/2})$ ($10^{8}$) is by combining the CholeskyQR2 with the
modified Gram-Schmidt re-orthogonalisation [19]. The basic version of the
algorithm proposed in [19] (without additional look-ahead optimizations) is
shown in Algorithm 6.

To make understanding more easy, the following notation is used in the rest of
the paper. In the description of the algorithms, the $A_{j}$ denotes the
$j$-th panel of the matrix $A$ with full row rank, and $A_{i:j}$ stands for
the range of panels from the $i$-th to the $j$-th with full row rank. In
addition, the term $A_{i,j}$ denotes the tile of dimension $(b\times b$),
where $b$ is the width of the panel.

Algorithm 6 Cholesky QR with Gram-Schmidt (CQRGS)

1:$A\in\mathbb{R}^{m\times n}$, panel width $b$ and number of panels
$k=\frac{n}{b}$

2:$Q\in\mathbb{R}^{m\times n}$ orthogonal and $R\in\mathbb{R}^{n\times n}$
upper triangular matrix

3:for $j=1\ldots k$ do

4: $W_{j}:=A_{j}^{T}A_{j}$ $\triangleright$ Construct Gram matrix

5: $W_{j}=U^{T}U$ $\triangleright$ Cholesky factorization

6: $Q_{j}=A_{j}U^{-1}$

7: $R_{j,j}=U$

8: $Y:=Q^{T}_{j}A_{j+1:k}$

9: $A_{j+1:k}:=A_{j+1:k}-Q_{j}Y$ $\triangleright$ Update panels

10: $R_{j,j+1:k}:=Y$

11:end for

The matrix $A$ is processed by panels, where $k$ is the number of panels and
$b$ is the panel width (line 1). The first step of each iteration is to
compute the Gram matrix (line 2), and the Cholesky factorization (line 3) of
the current panel $A_{j}$, followed by the construction of the orthogonal
matrix $Q_{j}$ (line 4), as in the CQR. As defined by the Gram-Schmidt
process, the rest of the panels to the right of the current panel are re-
orthogonalized (lines 6–7) w.r.t the $Q_{j}$ by applying the orthogonal
projection $Q_{j}Q^{T}_{j}$ to $A_{j+1:k}$. The advantage of the algorithm is
that all steps, except the Cholesky factorization in line 3, can be realized
using efficient level-3 BLAS kernels. In its original version [19], the
algorithm was developed and tested for shared memory and hybrid CPU-GPU
systems, where the Cholesky decomposition is performed on the CPU while the
trailing updates of $A$ are performed on the GPU. Since the algorithm operates
on panels with full row rank, parallelism was exploited only in the column
(panel) direction, depending on the high-performance Level-3 BLAS kernels used
for updating the trailing sub-matrices.

In the rest of this section, we describe our distributed CholeskyQR2 algorithm
with the modified Gram-Schmidt re-orthogonalization based on the algorithm
presented in [19]. Our algorithm is similar to the original CholeskyQR2
algorithm (see Algorithm 3), except that instead of calling CholeskyQR (CQR)
twice, our approach involves calling the CholeskyQR with Gram-Schmidt (CQRGS)
twice (see Algorithm 7).

Algorithm 7 CholeskyQR2 with Gram-Schmidt (CQR2GS)

1:$A\in\mathbb{R}^{m\times n}$

2:$Q\in\mathbb{R}^{m\times n}$ orthogonal and $R\in\mathbb{R}^{n\times n}$
upper triangular matrix

3:$[Q_{1},R_{1}]:=CQRGS(A)$

4:$[Q,R_{2}]:=CQRGS(Q_{1})$

5:$R:=R_{2}R_{1}$

The main contributions are extending the CholeskyQR and Gram-Schmidt (as
presented in [19]) to distributed memory systems equipped with GPU
accelerators and making some additional optimization that will make the
algorithm not only numerically stable for extremely ill-condition problems but
also in some cases superior in terms of performance compared to the original
CholeskyQR2.

Similar to CholeskyQR, the input matrix $A$ is partitioned and distributed
among $P$ processors in a one-dimensional block row layout (see Fig. 2
middle). Each block row $A_{p}$ is locally partitioned into $k$ panels
$A_{p,j}$ with width $b=n/k$ columns and $j\in\\{1,\ldots,k\\}$. The proposed
block row partitioning allows for coarse-grained parallelization between
processors (e.g. compute nodes) while partitioning into panels exploits fine-
grained parallelization at the node level.

Figure 2: Distributing and slicing of a matrix. An example with 2 processors.
The assignments of blocks and panels with processors are indicated on the
vertical axis.

Our distributed CholeskyQR algorithm with the blocked Gram-Schmidt re-
orthogonalization is described in Algorithm 8. The algorithm processes the
input matrix $A$ by panels, starting from the leftmost one. The first step is
to calculate the Gram matrix. To do this, each processor computes a local Gram
matrix from its current panel $A_{p,j}$ (line 2) and sums over all the local
matrices to get the final Gram matrix (line 3). Summing and distributing
requires collective communication (e.g. MPI_Allreduce), after which all
processors have identical Gram matrix. The Cholesky factorization (line 4) is
computed redundantly on each processor, and to avoid unnecessary
communication, each processor updates only its block row $Q_{p,j}$ of the
orthogonal panel $Q_{j}$ (line 5). Updating the panels right to the current
panel requires communication between the processors since the panel $Q_{j}$ is
applied to the full row rank of the trailing matrix $A_{p,j+1:k}$ (lines 7–9).
Each processor first calculates its partial product $Y_{p}$ (line 7), and then
a collective communication (line 8) is required to sum and broadcast the
global $Y$. Once the intermediate matrix $Y$ has been transmitted, each block
row can be updated independently (line 9).

Algorithm 8 Distributed Cholesky QR with blocked Gram-Schmidt

1:Number of processors $P$, $A\in\mathbb{R}^{m\times n}$ partitioned into
block rows and distributed among processors, panel width $b$

2:$Q\in\mathbb{R}^{m\times n}$ orthogonal and $R\in\mathbb{R}^{n\times n}$
upper triangular matrix

3:for j=1,2, …, k do

4: $W_{p,j}:=A_{p,j}^{T}A_{p,j}$

5: $W_{j}:=\text{MPI\\_Allreduce}(W_{p,j})$ $\triangleright$ Communication

6: $W_{j}=U^{T}U$

7: $Q_{p,j}:=A_{p,j}U^{-1}$

8: $R_{j,j}:=U$

9: $Y_{p}:=Q_{p,j}^{T}[A_{p,j+1},A_{p,j+2},\dots,A_{p,k}]$

10: $Y:=\text{MPI\\_Allreduce}(Y_{p})$ $\triangleright$ Communication

11: $[A_{p,j+1},\dots,A_{p,k}]:=[A_{p,j+1},\dots,A_{p,k}]-Q_{p,j}Y$

12: $[R_{j,j+1},R_{j,j+2},\dots,R_{j,k}]:=Y$

13:end for

A special case in which $b:=n$, leads to a single large panel eliminating the
iteration over $j$ and the panel updates, effectively removing the Gram-
Schmidt from the algorithm. In this case, the CQR2GS falls back to
CholeskyQR2, resulting in computational and communication costs equivalent to
those of CholeskyQR2.

A detailed breakdown of the computational and communication complexity for
each subroutine can be found in Table II. As shown in the table, the total
complexity cost of the parallel CholeskyQR2 with Gram-Schmidt (CQR2GS) is
dominated by expression $4\frac{m\ n^{2}}{P}$ that arises from the total cost
of routines syrk, trsm and update, that correspond to the construction of Gram
matrix, computing orthogonal matrix and updating part of Gram-Schmidt process,
respectively. Since it does not depend on the panel width $b$, the
computational complexity for these parts can easily be decreased by using more
processors and exploiting more parallelism. Opposite of that, the first and
last operands in the CQR2GS total computational costs, which correspond to the
Cholesky factorization and the matrix addition performed in the reduction
calls, depend on the panel width $b$. Therefore, $b$ can be considered as an
optimization parameter that can be adjusted to achieve either improved
numerical results, such as orthogonalization of the matrix, or optimized
computational performance by reducing the computational complexity. This
approach, characterized by a parameter $b$, not only increases the numerical
stability but also reduces the overall complexity compared to the conventional
CholeskyQR2 algorithm.

The reason why the numerical stability depends on the value of $b$ lies in the
fact that when input matrix $A$ is divided into panels, the condition number
of the first panel can be significantly decreased. Since the condition number
is the ratio between the largest and the smallest singular value, we have to
assess the upper and lower bounds of the singular values of the submatrix
(i.e. panel) $B\in\mathbb{R}^{m\times r}$ of the starting matrix
$A\in\mathbb{R}^{m\times n}$ with $n-r$ columns removed and $r<n$.

Let $\sigma_{1}\geq\sigma_{2}\geq\ldots\geq\sigma_{n-1}\geq\sigma_{n}$ be the
singular values of $A$ and
$\gamma_{1}\geq\gamma_{2}\geq\ldots\geq\gamma_{r-1}\geq\gamma_{r}$ singular
values of $B$, then the following inequalities hold:

$\displaystyle\sigma_{i}$ $\displaystyle\geq\gamma_{i},$ $\displaystyle
i=1,2,\ldots,r$ (3) $\displaystyle\gamma_{i}$
$\displaystyle\geq\sigma_{i+(n-r)},$ $\displaystyle i\leq r.$ (4)

For a proof see Colloraly 3.8.6 in [1]. From the above inequality, the largest
and the smallest singular values of $B$ are bounded by:

$\displaystyle\sigma_{1}$ $\displaystyle\geq\gamma_{1}\geq\sigma_{1+(n-r)}$
(5) $\displaystyle\sigma_{r}$ $\displaystyle\geq\gamma_{r}\geq\sigma_{n},$ (6)

and the condition number of $B$ is then:

$cond(A)=\frac{\sigma_{1}}{\sigma_{n}}\geq
cond(B)\geq\frac{\sigma_{1+(n-r)}}{\sigma_{r}}.$ (7)

In the worst case, the left equality holds and the condition number of the
panel is equal to the condition of the original matrix. This is the case of
very clustered singular values with an extremely large singular value where
the CQRGS (CQR2GS) cannot secure numerical stability. However, such extreme
use-cases are not the focus of this research and are the topic for future
work. In the general use-case, the condition number of $B$ is much smaller
than the upper bound, especially when the singular values of the original
matrix $A$ are equidistantly distributed. On this basis, the split of the
matrix, i.e. the panel width $b$ can be chosen so that the condition number of
the first panel $A_{1}$ is small enough so that its Gram matrix
$A_{1}^{T}A_{1}$ has a condition number less than $\mathcal{O}(u^{-1})$
($10^{15}$ in the case of double precision arithmetic).

TABLE II: Computational (up) and communication (down) complexity of CholeskyQR2 with Gram-Schmidt (GS). algorithm | routine | total comp |
---|---|---|---
CQR2GS | CQRGS | Gram | $b\ n\ \frac{m}{P}$ |
Gram_reduce | $b\ n\ log_{2}P$ |
Cholesky | $\frac{b^{2}\ n}{3}$ |
Construct_Q | $\frac{b\ m\ n}{P}$ |
GS | $2\frac{m\ n}{P}(n-b)$ |
GS_reduce | $\frac{n}{2}(n-b)\ log_{2}P$ |
| Total | $\frac{b^{2}\ n}{3}+2\frac{m\ n^{2}}{P}+\frac{n}{2}(n+b)\ log_{2}P$ |
Compute_R |  | $\frac{n^{3}}{3}$ |
Total |  | $2\frac{b^{2}\ n}{3}+\frac{n^{3}}{3}+4\frac{m\ n^{2}}{P}+n(n+b)\ log_{2}P$ |
|  |  | total comm | # of calls
CQR2GS | CQRGS | Gram_reduce | $b\ n\ log_{2}P$ | $\frac{n\ (n+b)}{2b^{2}}$
GS_reduce | $\frac{1}{2}n(n-b)\ log_{2}P$ | $\frac{n\ (n-b)}{2b^{2}}$
Total |  | $n\ (n+b)\ log_{2}P$ | $2\frac{n^{2}}{b^{2}}$

Too large $b$ can lead to a large condition number of the Gram matrix
$W_{j}=A_{j}^{T}A_{j}$ and to numerical instability of the Cholesky
factorization (see Fig. 3). For an ill-conditioned matrix with
$\kappa(A)=10^{\textbf{15}}$, the block size should be very small ($b=300$) to
achieve the orthogonality of order of ${u}$ resulting in dividing the matrix
into a larger number of panels (concretely 10) and iterating over lines 2–10
(Algorithm 8, line 1). On the other hand, for smaller condition numbers, the
panel width is increased and the number of iterations is decreased. Note that
in the case of matrices with condition numbers smaller than $O(10^{8})$, only
one panel is needed and CholeskyQR with Gram-Schmidt becomes a standard
CholeskyQR.

Figure 3: CQR2GS: Orthogonality of $Q$ as a function of panel size, for ill-
condition input matrices with $m=30000$, $n=3000$.

The connection between the number of panels and the performance (total
execution time) is shown in Fig. 4. It is observed that the panel width $b$
has a significant impact on the performance of the parallel CQR2GS and that
for a larger panel width (i.e. a smaller number of panels) the execution time
decreases. By decreasing the panel width (i.e. increasing the number of
panels), both the number of flops and the communication (volume of data)
decreases (see total flops, comm and #calls in Table II). However, with a
larger number of panels, the number of communication calls significantly
increases (with the factor $n^{3}/2)$, which results in a significantly longer
total execution time (the left side in Fig 4.

Figure 4: Time to solution of CQR2GS on 4 GPUs as a function of panel size for
well-condition input matrices ($\kappa(A)=10^{4}$) with $m=\\{30000,300000\\}$
and a fixed number of columns $n=3000$.

Another advantage of the CholeskyQR2 with Gram-Schmidt is the reduced amount
of communicated data, resulting in lower communication intensity compared to
CholeskyQR2. The total communication cost of the distributed CQRGS is:

$n\ (n+b)\operatorname{log}_{2}P,$ (8)

compared to the $2n^{2}\operatorname{log}_{2}P$ communication cost of
CholeskyQR2. Although additional communication is introduced by Gram-Schmidt
process (see Algorithm 8, line 8), a much significant communication decrease
is achieved in constructing Gram matrix with smaller $b$ resulting in a
smaller Gram matrix that has to be reduced and broadcasted among all MPI
ranks. This expression highlights that the volume of communicated data remains
consistently lower than that of the CholeskyQR2 algorithm as far as $b<<n$.

It is worth noting that although a smaller $b$ contributes to a reduction in
total costs, it leads to a higher number of memory operation calls
$\left(2\frac{n^{2}}{b^{2}}\right)$ if a larger number of panels is used,
which results in higher overhead. If the panel width $b$ is halved, the number
of communication messages increases by $4$ times. Despite the reduced flops
count for a smaller $b$, the influence of the communication on the total
execution time is more dominant than the improvements achieved by decreasing
the number of flops as illustrated in Fig. 4. Although the flops count
decreases with the smaller panel width (e.g. $b=100$), the total execution
time is much longer than if two panels are used (e.g. $b=1500).$

This highlights the nuanced optimization required in the selection of $b$ to
balance computational efficiency, communication overhead, and numerical
results. According to Fig. 3, for our specific input matrices, the algorithm
requires 10 panels to maintain orthogonality at a satisfactory level for input
matrices with condition number $10^{15}$. As a result, the algorithm itself
does not work in optimal mode. In this context, it would be beneficial to
investigate a modification that still works with good numerical behavior, but
with slightly larger values of $b$.

###  5.3 Modified CholeskyQR2 with Gram-Schmidt

The modified Cholesky-QR2 algorithm introduces a modification of CholeskyQR2
with Gram-Schmidt (mCQR2GS) by rearranging the order of computational tasks to
reduce the larger number of panels required in CQRGS. The computational
complexity and communication are therefore equivalent to the CQRGS algorithm
with the same number of panels. This idea provides an adaptive paneling
strategy that is tailored to the inherent properties of each matrix and
optimizes the factorization process for both ill- and well-conditioned
matrices. An example is shown in Fig. 5, where the matrix is divided into 3
panels.

Algorithm 9 Modified Cholesky QR with Gram-Schmidt (mCQRGS)

1:$A\in\mathbb{R}^{m\times n}$, number of panels $k$

2:$Q\in\mathbb{R}^{m\times n}$ orthogonal and $R\in\mathbb{R}^{n\times n}$
upper triangular matrix

3:$Q_{1},R_{1,1}=CQR2(A_{:,1})$ $\triangleright$ Orthogonalize first panel

4:for $j=2\ldots k$ do

5: $Y:=Q^{T}_{j-1}\ A_{:,j:k}$ $\triangleright$ Projections of orthogonal
panels on non-orthogonal

6: $A_{:,j:k}:=A_{:,j:k}-Q_{j-1}\ Y$ $\triangleright$ Update panels
$A_{j},\ldots A_{k}$

7: $[R_{j-1,j},R_{j-1,j+1},\ldots,R_{j-1,k}]:=Y$

8: $CQR(A_{:,j})$

9: $A_{:,j}:=A_{:,j}-Q_{1:j-1}\ Q^{T}_{1:j-1}\ A_{:,j}$ $\triangleright$
Reorthogonalize current panel

10: $Q_{j}=CQR(A_{:,j})$

11:end for

The input matrix $A$ is partitioned and distributed among $P$ processors in
the same manner as in CQRGS algorithm (see Fig. 2) with partitioned of local
block matrix $A_{p}$ into panels (Fig. 5). The pseudocode for the modified
algorithm is listed in Algorithm 9. Each rank starts with orthogonalizing the
first panel by applying full CQR2 (line 1) and then reorthogonalize the
trailing panels via the Gram-Schmidt process (line 3 i 4) in parallel. The
algorithm then moves by computing the QR of the second panel using CQR (line
6). However, the orthogonalization of this second panel is only partially
completed at this stage. Therefore, another step of reorthogonalization with
respect to the already orthogonal panels is required (line 7), followed by a
call to CQR to compute fully orthogonal $Q_{j}$ (line 8). This careful process
ensures the complete orthogonalization of the second panel and establishes its
orthogonality to the first panel. The algorithm then proceeds to the third
panel with the same routine as for the second panel.

Figure 5: Graphical overview of matrix distribution on 4 ranks and
computational operations on local matrix data.

Our analysis showed that the 3-panel strategy is optimal for our artificially
generated use cases as illustrated in Fig. 6. As expected, the mCQR2GS
algorithm still breaks down on our use-cases when the 2-panel strategy is used
on matrices with a very high condition number ($\geq 10^{15}$). The reason is
that the condition number of the first panel, according to the Eq. 7, is upper
bounded by the condition of the input matrix ($10^{15}$) and lower-bounded by
half of the condition number ($\approx 10^{8}$) which, in the best scenario of
reaching the lower bound, results in a Gram matrix with condition number
greater than $10^{15}$ and not fully semi-positive definite. Furthermore, the
mCQR2GS requires less flops compared to CQR2GS as illustrated in Fig. 7, as it
does not require the explicit construction of factor $R$ at the end (see Alg.
3 step 3). Note that for up to the condition number $10^{11}$ both algorithms
use 2-panel strategy as the optimal one. With the increasing condition number
CQR2GS requires a larger number of panels to secure the stability and
orthogonality, while mCQR2GS still requires only 2 panels, except for the last
use-case in which 3 panels are needed to improve the orthogonality.

Figure 6: Orthogonality of the $Q$ factor of mCQR2GS with 2 and 3 panels
w.r.t. to the condition number on one node and 4 GPUs. Matrix size $30k\times
3k$. Figure 7: Total execution time of CQR2GS and mCQR2GS when using the
optimal number of panels on one node and 4 GPUs. Matrix size $30k\times 3k$.

##  6 Scalability analysis

In the last part, we analyse the strong and weak scaling performance of the
Modified CholeskyQR2 with Gram-Schmidt (mCQR2QR) and compare it with
ScaLAPACK. In our test cases, the 3-panel strategy for mCQR2GS is used for
extremely ill-conditioned matrices. All tests achieve the required numerical
stability and orthogonality close to the machine precision $u$. In the case of
the GPU version, the NCCL communicator was used for collective communication
instead of the CUDA-aware MPI, as NCCL achieves much better performance by
significantly reducing the communication overhead.

The strong scaling was tested on the 3 artificial matrices (see subsection
2.2) with the condition number $\kappa(10^{4})$ and the dimensions $120k\times
1.2k$, $120k\times 6k$ and $120k\times 12k$ (the panels are $400,2000$ and
$4000$ wide respectively). Fig. 8 illustrates the strong scaling behaviour of
the CPU-only version of mCQR2GS and ScaLAPACK. Note that the scalability of
the mCQR2GS algorithm decreases with the number of nodes. The reason for this
is the communication (operation Allreduce), which is performed when building
the Gram matrix (Algorithm 8, line 3) and does not scale with the number of
nodes. In strong scaling tests, the width of the panel is fixed, which leads
to a constant load in Gram reduction operations, while the communication load
increases with the number of nodes, as shown in Fig. 9 (yellow line). Although
the computation parts scale close to the ideal line, the communication time
remains constant or slightly increases, resulting in an overall lower
scalability of the code as communication becomes more and more dominant, up to
$50\%$ of the total execution time (purple line). The observed peaks in
communication performance (Fig. 9, yellow line) on 4 and 8 nodes (16 and 32
MPI processes respectively) are the result of the internal NCCL optimisation
of allreduce operations implemented as a binary tree. For all tested matrices,
the mCQR2GS CPU outperforms ScaLAPACK by up to $4.7\times$.

Figure 8: Strong scaling of mCQR2GS-CPU and ScaLAPACK with artificial matrices
$m=120k$, $n=\\{1.2k,6k,12k\\}$ w.r.t the number of nodes. The dashed line is
the ideal scaling. Figure 9: Total execution time (green), execution times for
calculation (red) and communication (yellow) and share of communication in the
total execution time (purple) of mCQR2GS. Matrix size $120k\times 12k$.

Weak scaling experiments show the potential of the novel algorithm to perform
computations on large matrices while keeping the load per process or node
constant. The tests were performed on both CPU and GPU partitions of the Supek
supercomputer. The ScaLAPACK configuration is set to 16 tasks (MPIs) per node
and 8 threads per task, with the block height set to the number of rows
divided by the number of MPI processes and the block width set to 32. This
configuration achieved the best performance in our sweet spot analysis.

Fig. 10 shows that the weak scaling is nearly optimal for both the CPU and GPU
versions of mCQR2GS, with the total execution time ranging from $71.78$ msec
on one node to $124$ msec on 12 nodes for the GPU version. A steep time jump
from one to two nodes ($1.53\times$) is the NCCL communication overhead, which
was not observed when all NCCL processes are on the same node. Although the
introduction of inter-node communication is significant when switching to two
nodes, it remains almost constant when we introduce more nodes. Note that the
performance of ScaLAPACK is decreasing as the number of nodes increases. Since
the number of columns is fixed, the matrices become thinner and thinner as the
number of nodes increases (i.e. fixed number of rows per process/node),
resulting in a significant performance degradation for ScaLAPACK, which is
tailored for general and square matrices instead of tall and-skinny matrices.
Although the communication in mCQR2GS is significant compared to the total
execution time, both the GPU and CPU versions of the code achieve a
significant speedup compared to ScaLAPACK, with the CPU variant gaining more
than $6\times$ and the GPU variant $80\times$ in speedup.

Figure 10: Weak scaling of mCQR2GS-GPU, mCQR2GS-CPU and SCALAPACK with fixed
block size per node $40k\times 3k$ and 4 MPI processes per node. Dashed lines
are the ideal scaling.

##  7 Conclusion

In this paper, we have discussed the advantages and competitiveness of several
introduced algorithms, namely CQR2, sCQR3, CQR2GS and Householder QR, in terms
of key aspects including: suitability to perform in distributed environments,
numerical stability and the handling of tall-and-skinny matrices. Each of the
above algorithms exhibits strengths in one or more of these areas. CQR2, for
example, shows robust distributed parallelisation capabilities, but lack in
numerical stability. The sCQR3 algorithm has good numerical stability, but
with a tradeoff in terms of extra flops. On the other hand, while the
Householder-based QR provides good numerical stability, it performs sub-
optimally on tall and skinny matrices. We present a novel algorithm, which we
call modified CholeskyQR2 with Gram-Schmidt (mCQR2GS), to compute the QR
factorization of tall and skinny matrices on distributed multi-GPU
architectures. Our approach attempts to find a balance in terms of all the
above-mentioned features.

The novelty of the algorithm lies in complete orthogonalisation of the working
panel before performing the Gram-Schmidt step, which is not the case in the
original CholeskyQR2 with Gram-Schmidt. The novel algorithm achieves much
better numerical stability compared to other solutions tailored to tall-and-
skinny matrices, such as CholeskyQR2, especially in the case of extremely ill-
conditioned matrices (condition number of up to $10^{16}$). Our algorithm
proved to be faster than CQR2GS and outperforms ScaLAPACK on both distributed
CPU and GPU systems by $6\times$ and $80\times$, respectively.

The most important tuning parameter for performance and stability is the panel
width. By slicing the input matrix into smaller panels, the condition number
of each panel decreases. However, there is a tradeoff in terms of performance
when doing so. It is essential to tune this parameter to balance between
performance and stability and for our testing architecture and matrix dataset
(with equally distributed singular values), the optimal number of panels was
3. However, in the case of clustered singular values, our approach can not
improve the loss of orthogonality or numerical stability. To address this
issue, we plan to extend our solution with a shifting strategy to ensure
stability in such extreme use-cases.

The main bottleneck of our algorithm is the collective communication, whose
ratio to the total execution time increases by up to $50\%$ with an increasing
number of processes. An ongoing effort is the implementation of a look-ahead
approach to overlap the update of panels with computing the CholeskyQR of the
next panel (Algorithm 9 lines 4 and 6). The communication in CholeskyQR-based
algorithms is expensive compared to the computation, especially if CholeskyQR
has to be repeated multiple times to improve the orthogonality of the factor
$Q$, and increases with the number of processors. Moreover, the condition
number steeply decreases as we proceed with the panel processing, opening up a
space for further optimisation in reducing the number of flops by applying a
runtime decision on how many repetitions of CholeskyQR to perform.

The source code of the Shifted CholeskyQR3, CQR2GS and mCQR2GS versions is
available on the GitHub pages 111https://github.com/HybridScale/CholeskyQR2-IM
and Zenodo [23].

## Acknowledgments

This research was supported by the Croatian Science Foundation under grant
number HRZZ-UIP-2020-02-4559 and the European Regional Development Fund under
grant KK.01.1.1.01.009 - DATACROSS. All the computations were performed using
the Advanced computing service provided by the University of Zagreb University
Computing Centre - SRCE.

## References

  * [1] G. H. Golub and C. F. Van Loan, _Matrix Computations_ , 4th ed. Baltimore: The Johns Hopkinks University Press, 2013, vol. 37, no. 13.
  * [2] M. Heyouni and A. Essai, “Matrix Krylov subspace methods for linear systems with multiple right-hand sides,” _Numerical Algorithms_ , vol. 40, no. 2, pp. 137–156, 10 2005. [Online]. Available: https://link.springer.com/article/10.1007/s11075-005-1526-2http://link.springer.com/10.1007/s11075-005-1526-2
  * [3] M. H. Gutknecht, “Block Krylov Space Methods for Linear Systems With Multiple Right-hand Sides: An Introduction,” 2006.
  * [4] D. Davidovic and E. S. Quintana-Orti, “Applying OOC Techniques in the Reduction to Condensed Form for Very Large Symmetric Eigenproblems on GPUs,” in _2012 20th Euromicro International Conference on Parallel, Distributed and Network-based Processing_. IEEE, 2 2012, pp. 442–449. [Online]. Available: http://www.scopus.com/inward/record.url?eid=2-s2.0-84862136726&partnerID=tZOtx3y1http://ieeexplore.ieee.org/document/6169620/
  * [5] J. Winkelmann, P. Springer, and E. D. Napoli, “ChASE: Chebyshev Accelerated Subspace iteration Eigensolver for sequences of Hermitian eigenvalue problems,” _ACM Transactions on Mathematical Software_ , vol. 45, no. 2, pp. 1–34, 6 2019. [Online]. Available: http://dl.acm.org/citation.cfm?doid=3326465.3313828https://dl.acm.org/doi/10.1145/3313828
  * [6] J. Choi, J. Demmel, I. Dhillon, J. Dongarra, S. Ostrouchov, A. Petitet, K. Stanley, D. Walker, and R. C. Whaley, “ScaLAPACK: A portable linear algebra library for distributed memory computers - Design issues and performance,” _Computer Physics Communications_ , vol. 97, no. 1-2, pp. 1–15, 8 1996.
  * [7] A. Haidar, S. Tomov, P. Luszczek, and J. Dongarra, “MAGMA Embedded: Towards a Dense Linear Algebra Library for Energy Efficient Extreme Computing,” in _High Performance Extreme Computing Conference (HPEC)_ , 2015.
  * [8] J. Demmel, L. Grigori, M. Hoemmen, and J. Langou, “Communication-optimal parallel and sequential QR and LU factorizations,” _SIAM Journal on Scientific Computing_ , vol. 34, no. 1, pp. A206–A239, 1 2012. [Online]. Available: https://epubs.siam.org/doi/abs/10.1137/080731992http://arxiv.org/abs/0808.2664http://epubs.siam.org/doi/10.1137/080731992
  * [9] L. S. Blackford, A. Petitet, R. Pozo, K. Remington, R. C. Whaley, J. Demmel, J. Dongarra, I. Duff, S. Hammarling, G. Henry, and others, “An updated set of basic linear algebra subprograms (BLAS),” _ACM Transactions on Mathematical Software_ , vol. 28, no. 2, pp. 135–151, 2002.
  * [10] J. Demmel, L. Grigori, M. Hoemmen, and J. Langou, “Communication-avoiding parallel and sequential QR factorizations,” _CoRR abs/0806.2159_ , pp. 1–96, 2008. [Online]. Available: http://scholar.google.com/scholar?hl=en&btnG=Search&q=intitle:Communication-avoiding+parallel+and+sequential+QR+factorizations#0
  * [11] M. Gates, J. Kurzak, A. Charara, A. Yarkhan, and J. Dongarra, “SLATE: Design of a modern distributed and accelerated linear algebra library,” in _International Conference for High Performance Computing, Networking, Storage and Analysis, SC_. New York, NY, USA: IEEE Computer Society, 11 2019, pp. 1–18. [Online]. Available: https://dl.acm.org/doi/10.1145/3295500.3356223
  * [12] T. Fukaya, Y. Nakatsukasa, Y. Yanagisawa, and Y. Yamamoto, “CholeskyQR2: A Simple and Communication-Avoiding Algorithm for Computing a Tall-Skinny QR Factorization on a Large-Scale Parallel System,” in _2014 5th Workshop on Latest Advances in Scalable Algorithms for Large-Scale Systems_. IEEE, 11 2014, pp. 31–38. [Online]. Available: http://ieeexplore.ieee.org/document/7016731/
  * [13] Y. Yamamoto, Y. Nakatsukasa, Y. Yanagisawa, and T. Fukaya, “Roundoff error analysis of the Cholesky QR2 algorithm,” _Electronic Transactions on Numerical Analysis_ , vol. 44, pp. 306–326, 2015.
  * [14] E. Hutter and E. Solomonik, “Communication-avoiding CholeskyQR2 for rectangular matrices,” in _Proceedings - 2019 IEEE 33rd International Parallel and Distributed Processing Symposium, IPDPS 2019_. Institute of Electrical and Electronics Engineers Inc., 5 2019, pp. 89–100.
  * [15] T. Fukaya, R. Kannan, Y. Nakatsukasa, Y. Yamamoto, and Y. Yanagisawa, “Shifted cholesky qr for computing the qr factorization of ill-conditioned matrices,” _SIAM Journal on Scientific Computing_ , vol. 42, no. 1, pp. A477–A503, 2 2020. [Online]. Available: https://epubs.siam.org/doi/10.1137/18M1218212
  * [16] T. Terao, K. Ozaki, and T. Ogita, “LU-Cholesky QR algorithms for thin QR decomposition,” _Parallel Computing_ , vol. 92, p. 102571, 4 2020.
  * [17] A. J. Higgins, D. B. Szyld, E. G. Boman, and I. Yamazaki, “Analysis of Randomized Householder-Cholesky QR Factorization with Multisketching,” 9 2023. [Online]. Available: http://arxiv.org/abs/2309.05868
  * [18] I. Yamazaki, S. Tomov, and J. Dongarra, “Mixed-Precision Cholesky QR Factorization and Its Case Studies on Multicore CPU with Multiple GPUs,” _SIAM Journal on Scientific Computing_ , vol. 37, no. 3, pp. C307–C330, 1 2015. [Online]. Available: http://epubs.siam.org/doi/10.1137/14M0973773
  * [19] A. E. Tomás and E. S. Quintana-Ortí, “Cholesky and Gram-Schmidt Orthogonalization for Tall-and-Skinny QR Factorizations on Graphics Processors,” in _Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics)_ , vol. 11725 LNCS. Springer, 8 2019, pp. 469–480.
  * [20] A. Stathopoulos and K. Wu, “A Block Orthogonalization Procedure with Constant Synchronization Requirements,” _https://doi.org/10.1137/S1064827500370883_ , vol. 23, no. 6, pp. 2165–2182, 7 2006. [Online]. Available: https://epubs.siam.org/doi/10.1137/S1064827500370883
  * [21] L. Giraud, J. Langou, M. Rozložník, and J. V. D. Eshof, “Rounding error analysis of the classical Gram-Schmidt orthogonalization process,” _Numerische Mathematik_ , vol. 101, no. 1, pp. 87–100, 5 2005. [Online]. Available: http://link.springer.com/10.1007/s00211-005-0615-4
  * [22] T. Fukuya, R. Kannan, Y. Nakatsukasa, Y. Yamamoto, and Y. Yanagisawa, “Performance evaluation of the shifted Cholesky QR algorithm for ill-conditioned matrices,” in _sc18.supercomputing.org_. IEEE Computer Society, 2018.
  * [23] D. Davidović, N. Mijić, and A. K. Badrinarayanan, “CholeskyQR2-IM: CholeskyQR2 for ill-conditioned matrices v1.0.0,” 2024. [Online]. Available: https://zenodo.org/records/10888693

Nenad Mijić received his M.S. degree in Physics from the University of Zagreb
in 2017. He is currently working on his PhD under the supervision of Davor
Davidović. His research interests include code portability and eigenvalue
solvers on heterogeneous distributed systems.
---
Abhiram Kaushik is a postdoctoral researcher at the Ruđer Bošković Institute.
He received his PhD in Physics from the Indian Institute of Science in 2019.
His research interests include the application of high-performance computing
techniques to problems in high-energy physics and numerical linear algebra. He
has published 11 papers in peer-reviewed international scientific journals.
---
Davor Davidović is a senior research associate at the Ruđer Bošković
Institute. He received his PhD in computer science from the University of
Zagreb in 2014. His research interests are parallel and distributed computing,
code optimisation, hybrid and GPU programming, and scalable algorithms for
linear algebra. He has participated in more than 15 international research
projects (FP7, H2020, COST). He has published 34 original scientific papers in
international journals and in conference proceedings.
---