<figure><img src="image_1.png" /></figure>

HAL
open science

# Mixed precision iterative refinement for low-rank matrix and tensor approximations

Marc Baboulin, Oguz Kaya, Théo Mary, Matthieu Robeyns

<figure><img src="image_2.png" /></figure>

## To cite this version:

Marc Baboulin, Oguz Kaya, Théo Mary, Matthieu Robeyns. Mixed precision iterative refinement for low-rank matrix and tensor approximations. 2023. hal-04115337

HAL Id: hal-04115337

https://inria.hal.science/hal-04115337

Preprint submitted on 2 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

HAL
open science

# Mixed precision iterative refinement for low-rank matrix and tensor approximations

Marc Baboulin, Oguz Kaya, Theo Mary, Matthieu Robeyns

► To cite this version:

Marc Baboulin, Oguz Kaya, Theo Mary, Matthieu Robeyns. Mixed precision iterative refinement for low-rank matrix and tensor approximations. SIAM Journal on Scientific Computing, 2025, 47 (5), pp.A2906-A2935. 10.1137/23M1576992. hal-04115337v2

HAL Id: hal-04115337

https://hal.science/hal-04115337v2

Submitted on 9 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# MIXED PRECISION ITERATIVE REFINEMENT FOR LOW-RANK MATRIX AND TENSOR APPROXIMATIONS*

MARC BABOULIN†, OGUZ KAYA†, THEO MARY‡, AND MATTHIEU ROBEYNS†,§

**Abstract.** We present a new mixed precision algorithm to compute low-rank matrix and tensor approximations, a fundamental task in numerous applications in scientific computing and data analysis. Our algorithm is reminiscent of the iterative refinement framework for linear systems: we first compute a low-rank approximation in low precision and then refine its accuracy by iteratively updating it. We carry out an error analysis of our algorithm which proves that we can reach a high accuracy while performing most of the operations in low precision. We measure the computational cost of the algorithm, which depends on the numerical rank of the input (matrix or tensor) as well as the speed ratio between low and high precision arithmetic. We identify two situations where our method has a strong potential: when the hardware provides fast low precision matrix multiply-accumulate units, and when the numerical rank of the input is small at low accuracy levels. We confirm experimentally the potential of our algorithm for computing various low-rank matrix and tensor decompositions such as SVD, QR, Tucker, hierarchical Tucker, and tensor-train.

**Key words.** Matrix and tensor computations, low-rank approximations, mixed precision algorithms, iterative refinement, randomized SVD, tensor decompositions.

**AMS subject classifications.** 65F55, 65G50, 65Y20, 15A69

**1. Introduction.** Low-rank approximations (LRA) are a powerful tool used in many scientific applications to reduce the dimension of large scale data. For example, an $n \times n$ matrix $X$ may be approximated by a low-rank product $UV^T$ of $n \times r$ matrices $U$ and $V$, reducing the storage cost from $O(n^2)$ to $Onr)$. This storage cost is even more critical for tensors [20, 6, 14], the higher dimensional generalization of matrices : a $d$th order tensor $X$ requires storage in $O(n^d)$ that is exponential in the number of dimensions $d$. This *curse of dimensionality* can be tackled via LRA by decomposing the full tensor as a product of tensors of lower order and lower rank. Various such decompositions exist, such as the Tucker [28, 9], the tensor-train (TT) [27], and the hierarchical Tucker (HT) [15, 13] formats.

However, computing matrix or tensor low-rank decompositions at a controlled approximation accuracy $\varepsilon$ is computationally expensive; it represents the bottleneck of many LRA-based applications. Therefore, developing high performance algorithms for computing LRA is an important problem that has been the object of many studies; see, for instance, these surveys on matrices [19] and tensors [20, 14, 6].

Our work is specifically motivated by the recent and rapid emergence of *low precision arithmetics*, in particular half precision floating-point arithmetics such as the IEEE fp16 and bfloat16 formats. These low precision arithmetics provide great computational benefits, in terms of storage, speed, and energy [17]. On modern hardware, and in particular on Graphics Processing Unit (GPU) accelerators, low precision arithmetics can be orders of magnitude faster than the standard single (fp32) or double (fp64) precision arithmetics. However, using low precision also degrades the accuracy of the computations; for example half precision arithmetics provide at best between 3 and 4 digits of accuracy, depending on the format. This motivates

<footnote>*Version of June 2, 2023.</footnote>

<footnote>†University Paris-Saclay, CNRS, LISN, Gif-sur-Yvette, France, (marc.baboulin,oguz.kaya,matthieu.robeyns@universite-paris-saclay.fr)</footnote>

<footnote>‡Sorbonne Université, CNRS, LIP6, Paris, France, (theo.mary@lip6.fr)</footnote>

<footnote>§Corresponding author.</footnote>
# MIXED PRECISION ITERATIVE REFINEMENT FOR LOW-RANK MATRIX AND TENSOR APPROXIMATIONS*

MARC BABOULIN†, OGUZ KAYA†, THEO MARY‡, AND MATTHIEU ROBEYNS†,§

**Abstract.** We present a new mixed precision algorithm to compute low-rank matrix and tensor approximations, a fundamental task in numerous applications in scientific computing and data analysis. Our algorithm is reminiscent of the iterative refinement framework for linear systems: we first compute a low-rank approximation in low precision and then refine its accuracy by iteratively updating it. We carry out an error analysis of our algorithm which proves that we can reach a high accuracy while performing most of the operations in low precision. We measure the computational cost of the algorithm, which depends on the numerical rank of the input (matrix or tensor) as well as the speed ratio between low and high precision arithmetic. We identify two situations where our method has a strong potential: when the hardware provides fast low precision matrix multiply- accumulate units, and when the numerical rank of the input is small at low accuracy levels. We confirm experimentally the potential of our algorithm for computing various low-rank matrix and tensor decompositions such as SVD, QR, Tucker, hierarchical Tucker, and tensor-train.

**Key words.** Matrix and tensor computations, low-rank approximations, mixed precision algorithms, iterative refinement, randomized SVD, tensor decompositions.

**AMS subject classifications.** 65F55, 65G50, 65Y20, 15A69

**1. Introduction.** Low-rank approximations (LRA) are a powerful tool used in many scientific applications to reduce the dimension of large scale data. For example, an $n \times n$ matrix $X$ may be approximated by a low-rank product $UV^T$ of $n \times r$ matrices $U$ and $V$, reducing the storage cost from $O(n^2)$ to $O nr)$. This storage cost is even more critical for tensors [22, 8, 16], the higher dimensional generalization of matrices: a $d$th order tensor $X$ requires storage in $O(n^d)$ that is exponential in the number of dimensions $d$. This *curse of dimensionality* can be tackled via LRA by decomposing the full tensor as a product of tensors of lower order and lower rank. Various such decompositions exist, such as the Tucker [31, 11], the tensor-train (TT) [29], and the hierarchical Tucker (HT) [17, 15] formats.

However, computing matrix or tensor low-rank decompositions at a controlled approximation accuracy $\varepsilon$ is computationally expensive; it represents the bottleneck of many LRA-based applications. Therefore, developing high performance algorithms for computing LRA is an important problem that has been the object of many studies; see, for instance, these surveys on matrices [21] and tensors [22, 16, 8].

Our work is specifically motivated by the recent and rapid emergence of *low precision arithmetics*, in particular half precision floating-point arithmetics such as the IEEE fp16 and bfloat16 formats. These low precision arithmetics provide great computational benefits, in terms of storage, speed, and energy [19]. On modern hardware, and in particular on Graphics Processing Unit (GPU) accelerators, low precision arithmetics can be orders of magnitude faster than the standard single (fp32) or double (fp64) precision arithmetics. However, using low precision also degrades the accuracy of the computations; for example half precision arithmetics provide at best between 3 and 4 digits of accuracy, depending on the format. This motivates

<footnote>*Version of May 19, 2025.</footnote>

<footnote>†University Paris-Saclay, CNRS, LISN, Gif-sur-Yvette, France, (marc.baboulin,oguz.kaya, matthieu.robeyns@universite-paris-saclay.fr)</footnote>

<footnote>‡Sorbonne Université, CNRS, LIP6, Paris, France, (theo.mary@lip6.fr)</footnote>

<footnote>§Corresponding author.</footnote>
43

the need for mixed precision algorithms, which combine multiple precision formats
with the goal of achieving the high performance of the low precision while preserving
the high accuracy of the higher precision. The emergence of low precision on modern
hardware has generated much recent interest in mixed precision algorithms, with many
successful examples in numerical linear algebra; see the survey [17] for an overview of
this field.

49

In this paper, we seek to develop mixed precision algorithms for computing LRA.
Contrary to other linear algebra routines such as the solution of linear systems, there
has been relatively little work on designing mixed (or even low) precision algorithms
for LRA. Amestoy et al. [2] describe a mixed precision matrix LRA that partitions the
low-rank factors into several block columns stored in different precisions depending
on the singular values of the matrix; this approach can make use of low precisions
for matrices with rapidly decaying singular values. A similar approach is proposed
by Ooi et al. [23] for H-matrices. Connolly et al. [7] perform the rounding error
analysis of the randomized singular value decomposition (SVD) for matrix LRA, in
particular for the adaptive algorithm of Martinsson and Voronin [22] which builds an
LRA by increasing its rank until the prescribed accuracy is reached. The analysis in
[7] reveals that lower precision can be used when the norm of the approximation error
becomes small enough; this once more corresponds to the case where the matrix has
rapidly decaying singular values. Finally, a recent paper of Ootomo and Yokota [25]
also proposes a mixed precision randomized SVD; their approach stores the random
projection matrix in half precision while keeping the input matrix in single precision,
and obtains speedup by emulating single precision computations with GPU tensor
cores [11, 24]. There have also been very few attempts to develop mixed precision
LRA for tensors. We can mention a recent work by Yang et al. [30] that proposes an
iterative CP decomposition using mixed precision stochastic gradient descent. In a
context different from LRA, the recent work of Agullo et al. [1] is worth mentioning:
they consider the solution of linear systems Ax = b via GMRES, where A can be
approximated under tensor format [8]. This approach can then benefit from a mixed
precision implementation of GMRES.

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

In this paper, we propose a new method for computing LRA in mixed precision
arithmetic. Our approach is applicable to basically any LRA algorithm, involving
either matrices or tensors. It is reminiscent of the iterative refinement framework
used for solving linear systems: the idea is to first compute an LRA in low precision,
then evaluate the error (or residual) from this first LRA, and re-apply the same LRA
kernel to this error term to obtain a correction term that is used to refine the accuracy
of the LRA. This can be repeated iteratively to reach any level of desired accuracy.
The refined LRA is obtained as the sum of the original low precision LRA and the
correction term, and is thus of larger yet still of low-rank. In order to contain the
rank growth and maintain the optimal rank throughout the iterations, our method
employs a “recompression” strategy that is performed in high precision but whose
cost stays asymptotically smaller than that of LRA.

85

86

87

88

89

90

91

92

We carry out an error analysis of our method which proves that it can reach a
high accuracy while performing most of the operations in low precision. In particular,
we show that the precision used for the LRA kernel—which is the computational
bottleneck of the whole method—only affects the convergence speed of the process,
but not its attainable accuracy. In order to assess under which conditions we can
expect our method to be beneficial, we perform a complexity analysis that measures
the cost of the method as a function of the numerical rank of the input as well as the
speed ratio between the low and high precision arithmetic. We identify two situations
the need for mixed *precision algorithms*, which combine multiple precision formats with the goal of achieving the high performance of the low precision while preserving the high accuracy of the higher precision. The emergence of low precision on modern hardware has generated much recent interest in mixed precision algorithms, with many successful examples in numerical linear algebra; see the survey [19] for an overview of this field.

In this paper, we seek to develop mixed precision algorithms for computing LRA. Contrary to other linear algebra routines such as the solution of linear systems, there has been relatively little work on designing mixed (or even low) precision algorithms for LRA. Amestoy et al. [2] describe a mixed precision matrix LRA that partitions the low-rank factors into several block columns stored in different precisions depending on the singular values of the matrix; this approach can make use of low precisions for matrices with rapidly decaying singular values. A similar approach is proposed by Ooi et al. [25] for *H*-matrices. Connolly et al. [9] perform the rounding error analysis of the randomized singular value decomposition (SVD) for matrix LRA, in particular for the adaptive algorithm of Martinsson and Voronin [24] which builds an LRA by increasing its rank until the prescribed accuracy is reached. The analysis in [9] reveals that lower precision can be used when the norm of the approximation error becomes small enough; this once more corresponds to the case where the matrix has rapidly decaying singular values. Finally, a recent paper of Ootomo and Yokota [27] also proposes a mixed precision randomized SVD; their approach stores the random projection matrix in half precision while keeping the input matrix in single precision, and obtains speedup by emulating single precision computations with GPU tensor cores [13, 26]. There have also been very few attempts to develop mixed precision LRA for tensors. We can mention a recent work by Yang et al. [33] that proposes an iterative CP decomposition using mixed precision stochastic gradient descent. In a context different from LRA, the recent work of Agullo et al. [1] is worth mentioning: they consider the solution of linear systems $Ax = b$ via GMRES, where A can be approximated under tensor format [10]. This approach can then benefit from a mixed precision implementation of GMRES.

In this paper, we propose a new method for computing LRA in mixed precision arithmetic. Our approach is applicable to basically any LRA algorithm, involving either matrices or tensors. It is reminiscent of the iterative refinement framework used for solving linear systems: the idea is to first compute an LRA in low precision, then evaluate the error (or residual) from this first LRA, and re-apply the same LRA kernel to this error term to obtain a correction term that is used to refine the accuracy of the LRA. This can be repeated iteratively to reach any level of desired accuracy. The refined LRA is obtained as the sum of the original low precision LRA and the correction term, and is thus of larger yet still of low-rank. In order to contain the rank growth and maintain the optimal rank throughout the iterations, our method employs a “recompression” strategy that is performed in high precision but whose cost stays asymptotically smaller than that of LRA.

We carry out an error analysis of our method based on a general parameterized error model that only assumes that we have numerically stable implementations of the basic kernels used in our algorithm (LRA, matrix multiplication, and recompression). For matrices, many instances of such kernels are known to be stable. For tensors, the numerical stability of most tensor computations is still an open question, but extensive experiments certainly support the idea that they are indeed stable. Under this model, we prove that the method can reach a high accuracy while performing most of the operations in low precision. In particular, we show that the precision used for the
where our method has a strong potential. The first is when the hardware provides fast low precision matrix multiply-accumulate units, also called block FMA units [4], which allow for computing the low precision LRA at very high speed. The second is when the numerical rank of the input is small at low accuracy levels, which means that the singular values of the matrix or tensor are rapidly decaying; in this case, the first iterations of our method becomes inexpensive.

We apply our method to various low-rank matrix and tensor decompositions: SVD, QR, Tucker, HT, and TT. We perform some MATLAB experiments to confirm that our method is able to compute an LRA at any desired level of accuracy, while mostly using the low precision arithmetic.

The rest of this paper is organized as follows. In section 2, we describe our main algorithm of iterative refinement for LRA and provide the corresponding error and complexity analysis. Then we apply this method to various LRA algorithms for matrices and tensors in section 3 and section 4, respectively. We validate our method experimentally in section 5. Finally, we provide concluding remarks in section 6.

**1.1. Notations.** We denote as $X$ the object of interest, a matrix or a tensor, and as $F$ its low-rank factors that we seek to compute. If $X \in \mathbb{R}^{m \times n}$ is a matrix, then $F = UV^T$ is usually implicitly represented by the product of two matrices $U \in \mathbb{R}^{m \times r}$, $V \in \mathbb{R}^{n \times r}$ of smaller dimensions, where $r$ is the rank of $F$. Similarly, if $X \in \mathbb{R}^{n_1 \times \dots \times n_d}$ is a tensor of order $d$, its low-rank factors $F$ are represented using a number of matrices and/or tensors of smaller order; the precise expression of $F$ depends on the chosen format. In particular, in this article we will consider three different low-rank tensor decompositions:

* the Tucker format [9], which uses $d$ matrices $U^{(1)} \in \mathbb{R}^{n_1 \times r_1}, \dots, U^{(d)} \in \mathbb{R}^{n_d \times r_d}$, where $r_i \le n_i$, and a $d$-th order core tensor $G \in \mathbb{R}^{r_1 \times \dots \times r_d}$;

* the TT format [27], which uses $d$ 3rd-order tensors $U^{(1)} \in \mathbb{R}^{r_0 \times n_1 \times r_1}, U^{(2)} \in \mathbb{R}^{r_1 \times n_2 \times r_2}, \dots, U^{(d)} \in \mathbb{R}^{r_{d-1} \times n_d \times r_d}$, where $r_0 = r_d = 1$;

* the HT format [15, 13], which recursively applies the Tucker format to the core tensor $G$ to obtain a binary tree of tensors. Leaf nodes of the tree are matrices $U^{(1)} \in \mathbb{R}^{n_1 \times r_{d-1}}, \dots, U^{(d)} \in \mathbb{R}^{n_d \times r_{2d-2}}$ and the rest are 3rd-order tensors $G^{(1)} \in \mathbb{R}^{r_0 \times r_1 \times r_2}, \dots, G^{(d-1)} \in \mathbb{R}^{r_d \times r_{2d-3} \times r_{2d-2}}$, where $G^{(1)}$ is the root and $r_0 = 1$.

Note that our main iterative refinement method described in the next section does not require any specific knowledge or property of the tensor decomposition, and simply denotes its low-rank factors as $F$.

To control the accuracy of the LRA, we will distinguish two types of parameters, denoted as $\varepsilon$ and $u$. The truncation threshold $\varepsilon$ controls the low-rank truncation error and thus the accuracy of the LRA in exact arithmetic. The unit roundoff $u$ controls the rounding errors and thus the additional loss of accuracy incurred by the use of finite precision arithmetic. Actually, in our iterative refinement method, we will have four parameters $\varepsilon, u, \varepsilon_\ell$, and $u_\ell$, where $\varepsilon_\ell$ and $u_\ell$ denote the truncation threshold and unit roundoff used by the low precision LRA, whereas $\varepsilon$ and $u$ denote the truncation threshold and unit roundoff used by the rest of the operations performed in high precision, and correspond to the final target accuracy.

For a matrix $X$, we denote by $\text{RANK}(X, \varepsilon)$ the numerical rank of $X$ at accuracy $\varepsilon$, which is the smallest integer $r_\varepsilon$ such that there exists low-rank factors $F$ of rank $r_\varepsilon$ satisfying a relative accuracy $\varepsilon$, that is, $\|X - F\| \le \varepsilon \|X\|$ for a given choice of norm. This definition of the numerical rank extends to the Tucker, TT, and HT representations of tensors by defining $\text{RANK}(X, \varepsilon)$ as the vector of length $\tilde{d}$ whose
LRA kernel—which is the computational bottleneck of the whole method—only affects the convergence speed of the process, but not its attainable accuracy. Moreover, the method’s convergence speed is independent of the condition number of the matrix or tensor, which makes it widely applicable. In order to assess under which conditions we can expect our method to be beneficial, we perform a complexity analysis that measures the cost of the method as a function of the numerical rank of the input as well as the speed ratio between the low and high precision arithmetic. We identify two situations where our method has a strong potential. The first is when the hardware provides fast low precision matrix multiply-accumulate units, also called block FMA units [6], which allow for computing the low precision LRA at very high speed. The second is when the numerical rank of the input is small at low accuracy levels, which means that the singular values of the matrix or tensor are rapidly decaying; in this case, the first iterations of our method becomes inexpensive.

We apply our method to various low-rank matrix and tensor decompositions: SVD, QR, Tucker, HT, and TT. We perform some MATLAB experiments to confirm that our method is able to compute an LRA at any desired level of accuracy, while mostly using the low precision arithmetic.

The rest of this paper is organized as follows. In section 2, we describe our main algorithm of iterative refinement for LRA and provide the corresponding error and complexity analysis. Then we apply this method to various LRA algorithms, namely QRCP and randomized SVD for matrices, and HOSVD, TTSVD, and HTSVD for tensors in section 3 and section 4, respectively. We validate our method experimentally in section 5. Finally, we provide concluding remarks in section 6.

**1.1. Notations.** We denote as $X$ the object of interest, a matrix or a tensor, and as $F$ its low-rank factors that we seek to compute. If $X \in \mathbb{R}^{m \times n}$ is a matrix, then $F = UV^T$ is usually implicitly represented by the product of two matrices $U \in \mathbb{R}^{m \times r}$, $V \in \mathbb{R}^{n \times r}$ of smaller dimensions, where $r$ is the rank of $F$. Similarly, if $X \in \mathbb{R}^{n_1 \times \dots \times n_d}$ is a tensor of order $d$, its low-rank factors $F$ are represented using a number of matrices and/or tensors of smaller order; the precise expression of $F$ depends on the chosen format. In particular, in this article we will consider three different low-rank tensor decompositions:

* the Tucker format [11], which uses $d$ matrices $U^{(1)} \in \mathbb{R}^{n_1 \times r_1}, \dots, U^{(d)} \in \mathbb{R}^{n_d \times r_d}$, where $r_i \le n_i$, and a $d$-th-order core tensor $G \in \mathbb{R}^{r_1 \times \dots \times r_d}$;

* the TT format [29], which uses $d$ 3rd-order tensors $U^{(1)} \in \mathbb{R}^{r_0 \times n_1 \times r_1}, U^{(2)} \in \mathbb{R}^{r_1 \times n_2 \times r_2}, \dots, U^{(d)} \in \mathbb{R}^{r_{d-1} \times n_d \times r_d}$, where $r_0 = r_d = 1$;

* the HT format [17, 15], which recursively applies the Tucker format to the core tensor $G$ to obtain a binary tree of tensors. Leaf nodes of the tree are matrices $U^{(1)} \in \mathbb{R}^{n_1 \times r_{d-1}}, \dots, U^{(d)} \in \mathbb{R}^{n_d \times r_{2d-2}}$ and the rest are 3rd-order tensors $G^{(1)} \in \mathbb{R}^{r_0 \times r_1 \times r_2}, \dots, G^{(d-1)} \in \mathbb{R}^{r_d \times r_{2d-3} \times r_{2d-2}}$, where $G^{(1)}$ is the root and $r_0 = 1$.

Note that our main iterative refinement method described in the next section does not require any specific knowledge or property of the tensor decomposition, and simply denotes its low-rank factors as $F$.

Note that the low-rank factors $F$ are a structured (compressed) object representing a matrix or tensor of the same order and outer dimensions as $X$. Therefore, the subtraction $X - F$ is well-defined; since $X$ is uncompressed, evaluating the result of this subtraction requires to decompress $F$. Given two factors $F_1$ and $F_2$ representing matrices or tensors of the same order and outer dimensions, the addition $F_1 + F_2$ is also well-defined; since both $F_1$ and $F_2$ are compressed, the result of this addition can
coefficients correspond to the “inner” dimensions connecting the underlying matrices and/or tensors in the decomposition. The length $\bar{d}$ depends on both $d$ and the type of decomposition:

$$ \bar{d} = \begin{cases} d & (\text{Tucker}) \\ d-1 & (\text{TT}) \\ 2d-2 & (\text{HT}) \end{cases} \qquad (1.1) $$

Technically, the rank vector is determined by a vector of error tolerances $\epsilon_i$, $i = 1: \bar{d}$, whose sum yields the target accuracy $\epsilon = \sum_{i=1}^{\bar{d}} \epsilon_i$. In practice, the truncation error is uniformly distributed by setting $\epsilon_i = \epsilon/\bar{d}$ [27, 29, 13].

## 2. Iterative refinement for low-rank approximation.
In this section, we propose a mixed precision algorithm to compute LRA at high accuracy while performing most of the computations in low precision.

We consider the problem of computing low-rank factors $F$ of a matrix or tensor $X$ satisfying a relative accuracy $\epsilon$, that is, $\|X - F\| \le \epsilon\|X\|$. In finite precision arithmetic, we must carefully select the precision in order to achieve this accuracy. In a uniform precision context (where we perform all computations in the same precision), we must use a precision whose unit roundoff $u$ is safely smaller than $\epsilon$, that is, $u = \theta\epsilon$, with $\theta \le 1$ some parameter. Indeed, we may then expect the computed factors $\hat{F}$ to satisfy

$$ \|X - \hat{F}\| \le (\epsilon + cu)\|X\| = (1 + c\theta)\epsilon\|X\| $$
for some constant $c$, where the term $c\theta$ accounts for the effect of rounding errors and can be made as small as needed by decreasing $\theta$ (decreasing $u$, that is, increasing the precision).

For applications requiring relatively high accuracy (small values of $\epsilon$), this uniform precision approach therefore cannot make use of lower precisions. This motivates us to propose a mixed precision method, outlined below, which uses low precision to compute a first approximate set of factors $F_0$, and then refines them into more accurate factors $F_1$ as follows.

1. Compute low-rank factors $F_0$ of $X$ at low accuracy $\epsilon_\ell$ and in low precision $u_\ell = \theta\epsilon_\ell$.

2. Compute the error $E = X - F_0$ in high precision $u$.

3. Compute low-rank factors $F_E$ of $E$ at low accuracy $\epsilon_\ell$ and in low precision $u_\ell = \theta\epsilon_\ell$.

4. Update the low-rank factors of $X$ to $F_1 = F_0 + F_E$ in high precision $u$.

This approach is based on the observation that, while the first factors $F_0$ will achieve a low accuracy $\epsilon_\ell$ relative to $\|X\|$,

$$ \|X - F_0\| \le \epsilon_\ell \|X\|, $$

the factors $F_E$ of the error $E$ will achieve a low accuracy $\epsilon_\ell$ relative to $\|E\|$,

$$ \|E - F_E\| \le \epsilon_\ell \|E\|. $$

Neglecting for now the effect of rounding errors, we can expect to have $\|E\| \approx \epsilon_\ell \|X\|$, and therefore the combined factors $F_0 + F_E$ will satisfy an accuracy of about $\epsilon_\ell^2$:

$$ \|X - F_1\| \lesssim \epsilon_\ell^2 \|X\|. $$

This idea is reminiscent of the iterative refinement for linear systems $Ax = b$, where step 1 corresponds to computing an initial solution $x_0$, step 2 corresponds to
be formed by using low-rank arithmetic, that is, can be represented by a low-rank factor $F = F_1 + F_2$ with larger ranks.

To control the accuracy of the LRA, we will distinguish two types of parameters, denoted as $\varepsilon$ and $u$. The truncation threshold $\varepsilon$ controls the low-rank truncation error and thus the accuracy of the LRA in exact arithmetic. The unit roundoff $u$ controls the rounding errors and thus the additional loss of accuracy incurred by the use of finite precision arithmetic. Actually, in our iterative refinement method, we will have four parameters $\varepsilon$, $u$, $\varepsilon_\ell$, and $u_\ell$, where $\varepsilon_\ell$ and $u_\ell$ denote the truncation threshold and unit roundoff used by the low precision LRA, whereas $\varepsilon$ and $u$ denote the truncation threshold and unit roundoff used by the rest of the operations performed in high precision, and correspond to the final target accuracy.

For a matrix $X$, we denote by $\text{RANK}(X, \varepsilon)$ the numerical rank of $X$ at accuracy $\varepsilon$, which is the smallest integer $r_\varepsilon$ such that there exists low-rank factors $F$ of rank $r_\varepsilon$ satisfying a relative accuracy $\varepsilon$, that is, $\|X - F\| \le \varepsilon\|X\|$ for a given choice of norm. We also define $\text{RANK}(X)$, the non-numerical (or exact) rank of a matrix $X$. This definition of the numerical rank extends to the Tucker, TT, and HT representations of tensors by defining $\text{RANK}(X, \varepsilon)$ and $\text{RANK}(X)$ as the vector of length $\bar{d}$ whose coefficients correspond to the “inner” dimensions connecting the underlying matrices and/or tensors in the decomposition. The length $\bar{d}$ depends on both $d$ and the type of decomposition:

$$ \bar{d} = \begin{cases} d & (\text{Tucker}) \\ d-1 & (\text{TT}) \\ 2d-2 & (\text{HT}) \end{cases} \qquad (1.1) $$

Technically, the rank vector is determined by a vector of error tolerances $\varepsilon_i$, $i = 1: \bar{d}$, whose sum yields the target accuracy $\varepsilon = \sum_{i=1}^{\bar{d}} \varepsilon_i$. In practice, the truncation error is uniformly distributed by setting $\varepsilon_i = \varepsilon/\bar{d}$ [29, 32, 15].

## 2. Iterative refinement for low-rank approximation.
In this section, we propose a mixed precision algorithm to compute LRA at high accuracy while performing most of the computations in low precision.

We consider the problem of computing low-rank factors $F$ of a matrix or tensor $X$ satisfying a relative accuracy $\varepsilon$, that is, $\|X - F\| \le \varepsilon\|X\|$. In finite precision arithmetic, we must carefully select the precision in order to achieve this accuracy. In a uniform precision context (where we perform all computations in the same precision), we must use a precision whose unit roundoff $u$ is safely smaller than $\varepsilon$, that is, $u = \theta\varepsilon$, with $\theta \le 1$ some parameter. Indeed, we may then expect the computed factors $\hat{F}$ to satisfy

$$ \|X - \hat{F}\| \le (\varepsilon + cu)\|X\| = (1 + c\theta)\varepsilon\|X\| $$

for some constant $c$, where the term $c\theta$ accounts for the effect of rounding errors and can be made as small as needed by decreasing $\theta$ (decreasing $u$, that is, increasing the precision).

For applications requiring relatively high accuracy (small values of $\varepsilon$), this uniform precision approach therefore cannot make use of lower precisions. This motivates us to propose a mixed precision method, outlined below, which uses low precision to compute a first approximate set of factors $F_0$, and then refines them into more accurate factors $F_1$ as follows.

1. Compute low-rank factors $F_0$ of $X$ at low accuracy $\varepsilon_\ell$ and in low precision $u_\ell = \theta\varepsilon_\ell$.

2. Compute the error $E = X - F_0$ in high precision $u$.
computing the residual $r = b - Ax_0$, and step 3 corresponds to solving another linear system $Ad = r$ for a correction term $d$, which is used to update $x_1 = x_0 + d$ (our step 4). We therefore baptize this approach “iterative refinement” for LRA.

One of the crucial insights that makes this method effective is that the error $E$ has a low numerical rank whenever $X$ has one too, since it is the sum of $X$ and $F_0$, which is low-rank by construction. To be specific, we will prove in subsection 2.2 that $\text{RANK}(E, \varepsilon_\ell)$ is at most $2\text{RANK}(X, \varepsilon_\ell^2)$ and, for the same reason, the rank of the refined factors $F_1$ is bounded by $3\text{RANK}(X, \varepsilon_\ell^2)$. The factors $F_1$ can then be cheaply recompressed into factors of optimal rank $\text{RANK}(X, \varepsilon_\ell^2)$ in high precision $u$ by exploiting their (suboptimal) low-rank structure.

If accuracy higher than $\varepsilon_\ell^2$ is needed, we can apply the method again on $F_1$ to obtain an accuracy of $\varepsilon_\ell^3$, and so on; we obtain Algorithm 2.1, which repeats this process until the desired accuracy is achieved. Algorithm 2.1 is described within a general iterative refinement framework that uses an arbitrary LRA algorithm: $\text{LRA}(X, \varepsilon)$ returns low-rank factors $F$ of $X$ at accuracy $\varepsilon$. We also require a $\text{DECOMPRESS}$ kernel which transforms low-rank factors $F$ back to a full-size object, and a $\text{RECOMPRESS}$ kernel which takes low-rank factors $F$ and $\varepsilon$ as input, and computes their optimal LRA at accuracy $\varepsilon$.

<caption>**Algorithm 2.1** Iterative refinement for LRA.</caption>

**Input:** $X$, the matrix or tensor of interest; 
$\varepsilon$, the desired relative accuracy for the approximation; 
$n_{\text{it}}$, the maximum number of iterations; 
LRA, a low-rank approximation algorithm; 
DECOMPRESS, RECOMPRESS: decompression and recompression algorithms.

**Output:** $F$, the low-rank factors $X$.

1: Compute $F = \text{LRA}(X, \varepsilon_\ell)$ in precision $u_\ell$.
2: **for** $i = 1$ **to** $n_{\text{it}}$ **do**
3: &nbsp;&nbsp;&nbsp;&nbsp;Compute $E = X - \text{DECOMPRESS}(F)$ in precision $u$.
4: &nbsp;&nbsp;&nbsp;&nbsp;Compute $\alpha = \|E\|$ in precision $u$.
5: &nbsp;&nbsp;&nbsp;&nbsp;If $\alpha \le \varepsilon\|X\|$, exit.
6: &nbsp;&nbsp;&nbsp;&nbsp;Scale $E \leftarrow \alpha^{-1}E$ in precision $u$.
7: &nbsp;&nbsp;&nbsp;&nbsp;Compute $F_E = \text{LRA}(E, \varepsilon_\ell)$ in precision $u_\ell$.
8: &nbsp;&nbsp;&nbsp;&nbsp;Scale back $F_E \leftarrow \alpha F_E$ in precision $u$.
9: &nbsp;&nbsp;&nbsp;&nbsp;Compute $F = \text{RECOMPRESS}(F + F_E, \varepsilon)$ in precision $u$.
10: **end for**

Algorithm 2.1 incorporates two additional steps to make the method effective. First, the error $E$ is scaled by the inverse of its norm (Line 6) before computing its low-rank factors $F_E$, and then scaled back after (Line 8). This is done to prevent the elements of $E$ from underflowing when converted to the lower precision, in the case where the arithmetic uses not only a reduced number of bits in the significand but also in the exponent, such as is the case for the IEEE half precision fp16 arithmetic. Second, the factors $F$ are recompressed (Line 9) to their optimal low-rank representation to avoid the rank growing out of control during the iterations. Indeed, as mentioned above, the rank of $F$ is bounded by $3\text{RANK}(X, \varepsilon_\ell^2)$ after one iteration. Therefore, in absence of recompression, the rank would grow as $3^k \text{RANK}(X, \varepsilon_\ell^2)$ after $k$ iterations, and would quickly make the method unaffordable. Fortunately, this recompression can be performed inexpensively for most LRA algorithms of interest.
3. Compute low-rank factors $F_E$ of $E$ at low accuracy $\varepsilon_\ell$ and in low precision $u_\ell = \theta\varepsilon_\ell$.

4. Update the low-rank factors of $X$ to $F_1 = F_0 + F_E$ in high precision $u$.

This approach is based on the observation that, while the first factors $F_0$ will achieve a low accuracy $\varepsilon_\ell$ relative to $\|X\|_2$,

$$ \|X - F_0\| \le \varepsilon_\ell \|X\|, $$

the factors $F_E$ of the error $E$ will achieve a low accuracy $\varepsilon_\ell$ relative to $\|E\|_2$,

$$ \|E - F_E\| \le \varepsilon_\ell \|E\|. $$

Neglecting for now the effect of rounding errors, we have $\|E\| = \|X - F_0\| \le \varepsilon_\ell \|X\|$, and therefore the combined factors $F_0 + F_E$ will satisfy an accuracy of about $\varepsilon_\ell^2$:

$$ \|X - F_1\| = \|X - F_0 - F_E\| = \|E - F_E\| \le \varepsilon_\ell \|E\| \le \varepsilon_\ell \|X - F_0\| \le \varepsilon_\ell^2 \|X\|. $$

This idea is therefore based on computing additive updates $F_0 + F_E$ in mixed precision. This is reminiscent of the iterative refinement for linear systems $Ax = b$, which is also based on additive updates of the solution $x$. Indeed, for $Ax = b$, step 1 corresponds to computing an initial solution $x_0$, step 2 corresponds to computing the residual $r = b - Ax_0$, step 3 corresponds to solving another linear system $Ad = r$ for a correction term $d$, and step 4 updates $x_1 = x_0 + d$. The term “iterative refinement” is most commonly used for the solution of linear systems, but is also found in least squares, singular value, or eigenvalue problems. Therefore, we believe that the method proposed in this article should be called “iterative refinement for LRA”.

One of the crucial insights that makes this method effective is that the error $E$ has a low numerical rank whenever $X$ has one too, since it is the sum of $X$ and $F_0$, which is low-rank by construction. To be specific, we will prove in subsection 2.2 that $\mathrm{RANK}(E, \varepsilon_\ell)$ is at most $2 \mathrm{RANK}(X, \varepsilon_\ell^2)$ and, for the same reason, the rank of the refined factors $F_1$ is bounded by $3 \mathrm{RANK}(X, \varepsilon_\ell^2)$. The factors $F_1$ can then be cheaply recompressed into factors of optimal rank $\mathrm{RANK}(X, \varepsilon_\ell^2)$ in high precision $u$ by exploiting their (suboptimal) low-rank structure.

If accuracy higher than $\varepsilon_\ell^2$ is needed, we can apply the method again on $F_1$ to obtain an accuracy of $\varepsilon_\ell^3$, and so on; we obtain Algorithm 2.1, which repeats this process until the desired accuracy is achieved.

Algorithm 2.1 is described within a general iterative refinement framework that uses an arbitrary LRA algorithm: LRA($X, \varepsilon$) returns low-rank factors $F$ of $X$ at accuracy $\varepsilon$. We also require a DECOMPRESS kernel which transforms low-rank factors $F$ back to a full-size object, and a RECOMPRESS kernel which takes low-rank factors $F$ and $\varepsilon$ as input, and computes their optimal LRA at accuracy $\varepsilon$ (that is, the LRA of the smallest rank that achieves an accuracy of at least $\varepsilon$). We consider these three abstract kernels in order for Algorithm 2.1 to be as general as possible. In particular, this will allow us to apply it to different LRA methods, both for matrices and tensors. However, to be concrete, we will also discuss specific examples in sections 3 and 4. Namely, we will consider five different LRA kernels: QRCP and randomized SVD for matrices, and HOSVD, TTSVD, and HTSVD for tensors. For matrices, the DECOMPRESS kernel is simply a matrix–matrix product; for tensors, it is similarly a contraction of the low-rank tensor along all inner dimensions. Finally, the RECOMPRESS kernel depends on the LRA kernel: for each case, we will show how to adapt the LRA kernel (which takes as input a full object) to obtain the RECOMPRESS kernel (which takes as input a low-rank object); see sections 3 and 4 for more details.
Algorithm 2.1 has four parameters that control its accuracy: $\varepsilon$, $\varepsilon_\ell$, $u$, and $u_\ell$.

* $\varepsilon$ indicates the target accuracy for the final factors, and is prescribed by the user as input to Algorithm 2.1.

* $u$ is the unit roundoff of the high precision, which should be taken to be the lowest possible such that $u$ is still safely smaller than $\varepsilon$: $u = \theta\varepsilon$, $\theta \le 1$.

* $u_\ell$ is the unit roundoff of the low precision, which is used to perform the most expensive parts of the computation, the calls to LRA. Its choice depends on both the available arithmetics on the target hardware, and the target accuracy $\varepsilon$. Indeed, lowering the precision makes each iteration faster but requires more of them.

* Finally, $\varepsilon_\ell$ is the tolerance used by the LRA kernel; since LRA is performed in precision $u_\ell$, $\varepsilon_\ell$ should be set such that $u_\ell$ is safely smaller than $\varepsilon_\ell$, that is, $\varepsilon_\ell = u_\ell/\theta_\ell$, $\theta_\ell \le 1$. We note that this is necessary because using an $\varepsilon_\ell$ too close to $u_\ell$ prevents most LRA algorithms to reliably detect the correct rank, due to the noise introduced by rounding errors.

In the rest of this section, we first perform an error analysis to determine the attainable accuracy and convergence rate of Algorithm 2.1. We then perform a complexity analysis to determine under which conditions the algorithm can be expected to be faster than a standard uniform precision LRA performed entirely in high precision $u$.

## 2.1. Error analysis.
In order to carry out the error analysis of Algorithm 2.1, we will use the standard model of floating-point arithmetic [16, sect. 2.2]. In addition, we also need to make the following three assumptions on the numerical behavior of the LRA, DECOMPRESS, and RECOMPRESS kernels. First, we assume that computing $F = \text{LRA}(X, \varepsilon_\ell)$ in precision $u_\ell$ yields computed factors $\hat{F}$ satisfying

$$ \|X - \hat{F}\| \le (\varepsilon_\ell + b_1 u_\ell) \|X\| = (1 + b_1 \theta_\ell) \varepsilon_\ell \|X\|. \quad (2.1) $$

Second, we assume that computing $F = \text{DECOMPRESS}(F)$ in precision $u$ yields computed factors $\hat{F}$ satisfying

$$ \|F - \hat{F}\| \le b_2 u \|F\| = b_2 \theta \varepsilon \|F\|. \quad (2.2) $$

Third, we assume that computing $F = \text{RECOMPRESS}(F, \varepsilon)$ in precision $u$ yields computed factors $\hat{F}$ satisfying

$$ \|F - \hat{F}\| \le (\varepsilon + b_3 u) \|F\| = (1 + b_3 \theta) \varepsilon \|F\|. \quad (2.3) $$

For simplicity, we ignore any rounding errors associated with the scaling by $\alpha = \|E\|$, which are negligible and can in fact be prevented by rounding $\alpha$ to the nearest power of two (in binary floating-point arithmetic).

We are now ready to prove the following result.

**THEOREM 2.1.** If Algorithm 2.1 is applied to $X$ with LRA, DECOMPRESS, and RECOMPRESS kernels satisfying (2.1)-(2.3), then after $k$ iterations, the computed factors $\hat{F}$ satisfy

$$ \|X - \hat{F}\| \le (\phi^{k+1} + \xi + O(\varepsilon_\ell \varepsilon)) \|X\|, \quad (2.4) $$

with $\phi = (1 + b_1 \theta_\ell) \varepsilon_\ell + (b_2 + 1) \theta_\ell \varepsilon$ and $\xi = (1 + (b_2 + b_3 + 2) \theta_\ell) \varepsilon$.

*Proof.* Defining $\hat{F}_i$ as the computed factors after $i$ iterations, our goal is to obtain a bound of the form

$$ \|X - \hat{F}_{i+1}\| \le \phi \|X - \hat{F}_i\| + \xi \|X\|, \quad \phi < 1, $$
<caption>**Algorithm 2.1** Iterative refinement for LRA.</caption>

**Input:** $X$, the matrix or tensor of interest;
$\varepsilon$, the desired relative accuracy for the approximation;
$n_{it}$, the maximum number of iterations;
LRA, a low-rank approximation algorithm;
DECOMPRESS, RECOMPRESS: decompression and recompression algorithms.
**Output:** $F$, the low-rank factors $X$.
1: Compute $F = \text{LRA}(X, \varepsilon)$ in precision $u_\ell$.
2: **for** $i = 1$ **to** $n_{it}$ **do**
3: &nbsp;&nbsp;&nbsp;&nbsp;Compute $E = X - \text{DECOMPRESS}(F)$ in precision $u$.
4: &nbsp;&nbsp;&nbsp;&nbsp;Compute $\alpha = \|E\|$ in precision $u$.
5: &nbsp;&nbsp;&nbsp;&nbsp;If $\alpha \le \varepsilon \|X\|$, exit.
6: &nbsp;&nbsp;&nbsp;&nbsp;Scale $E \leftarrow \alpha^{-1} E$ in precision $u$.
7: &nbsp;&nbsp;&nbsp;&nbsp;Compute $F_E = \text{LRA}(E, \varepsilon)$ in precision $u_\ell$.
8: &nbsp;&nbsp;&nbsp;&nbsp;Scale back $F_E \leftarrow \alpha F_E$ in precision $u$.
9: &nbsp;&nbsp;&nbsp;&nbsp;Compute $F = \text{RECOMPRESS}(F + F_E, \varepsilon)$ in precision $u$.
10: **end for**

Algorithm 2.1 incorporates two additional steps to make the method effective. First, the error $E$ is scaled by the inverse of its norm (Line 6) before computing its low-rank factors $F_E$, and then scaled back after (Line 8). This is done to prevent the elements of $E$ from underflowing when converted to the lower precision, in the case where the arithmetic uses not only a reduced number of bits in the significand but also in the exponent, such as is the case for the IEEE half precision fp16 arithmetic. We also note that the elements of $X$ itself are also susceptible to underflow or overflow; in this case, a similar scaling should be applied to $X$. Second, the factors $F$ are recompressed (Line 9) to their optimal low-rank representation to prevent the rank from growing out of control during the iterations. Indeed, as mentioned above, the rank of $F$ is bounded by $3 \text{RANK}(X, \varepsilon_\ell^2)$ after one iteration. Therefore, in the absence of recompression, the rank would grow as $3^k \text{RANK}(X, \varepsilon_\ell^2)$ after $k$ iterations, and quickly make the method unaffordable. This RECOMPRESS kernel can be implemented for an asymptotically negligible cost for many LRA algorithms, including all those considered in this article.

Algorithm 2.1 has four parameters that control its accuracy: $\varepsilon$, $\varepsilon_\ell$, $u$, and $u_\ell$.

* $\varepsilon$ indicates the target accuracy for the final factors, and is prescribed by the user as input to Algorithm 2.1.

* $u$ is the unit roundoff of the high precision, which should be taken to be the lowest possible such that $u$ is still safely smaller than $\varepsilon$: $u = \theta\varepsilon$, $\theta \le 1$.

* $u_\ell$ is the unit roundoff of the low precision, which is used to perform the most expensive parts of the computation, the calls to LRA. Its choice depends on both the available arithmetics on the target hardware, and the target accuracy $\varepsilon$. Indeed, lowering the precision makes each iteration faster but requires more of them.

* Finally, $\varepsilon_\ell$ is the tolerance used by the LRA kernel; since LRA is performed in precision $u_\ell$, $\varepsilon_\ell$ should be set such that $u_\ell$ is safely smaller than $\varepsilon_\ell$, that is, $\varepsilon_\ell = u_\ell/\theta_\ell$, $\theta_\ell \le 1$. We note that this is necessary because using an $\varepsilon_\ell$ too close to $u_\ell$ prevents most LRA algorithms from reliably detecting the correct rank, due to the noise introduced by rounding errors.
which will allow us to conclude that the error contracts by a factor $\phi$ at each iteration, until it converges to its maximum attainable accuracy $\xi$.

Given $\hat{F}_i$, the first step is to decompress it and compute $E$ at Line 3; by (2.2), the computed $\hat{E}$ satisfies

$$ \|X - \hat{F}_i - \hat{E}\| \le \theta\varepsilon((b_2 + 1)\|\hat{F}_i\| + \|X\|), \quad (2.5) $$

with an extra term $\theta\varepsilon(\|\hat{F}_i\| + \|X\|)$ on the right-hand side coming from the rounding error incurred by the subtraction. Using the triangle inequality $\|\hat{F}_i\| \le \|X - \hat{F}_i\| + \|X\|$, we can rearrange (2.5) as

$$ \|X - \hat{F}_i - \hat{E}\| \le \theta\varepsilon((b_2 + 1)\|X - \hat{F}_i\| + (b_2 + 2)\|X\|). \quad (2.6) $$

Then we compute LRA($\hat{E}$, $\varepsilon_\ell$) at Line 7; by (2.1) the computed $\hat{F}_E$ satisfies

$$ \|\hat{E} - \hat{F}_E\| \le (1 + b_1\theta_\ell)\varepsilon_\ell\|\hat{E}\|. \quad (2.7) $$

By using (2.6) and the triangle inequality, we rearrange (2.7) as

$$ \|\hat{E} - \hat{F}_E\| \le (1 + b_1\theta_\ell)\varepsilon_\ell(\|X - \hat{F}_i\| + \|X - \hat{F}_i - \hat{E}\|) \quad (2.8) $$

$$ = (1 + b_1\theta_\ell)\varepsilon_\ell\|X - \hat{F}_i\| + O(\varepsilon_\ell\varepsilon), \quad (2.9) $$

where we do not keep track explicitly of high order terms in $O(\varepsilon_\ell\varepsilon)$ for the sake of readability. Finally, we obtain the next iterate $F_{i+1}$ by recompressing $\hat{F}_i + \hat{F}_E$ at Line 9; by (2.3), the computed $\hat{F}_{i+1}$ satisfies

$$ \|\hat{F}_i + \hat{F}_E - \hat{F}_{i+1}\| \le (1 + b_3\theta)\varepsilon\|\hat{F}_i + \hat{F}_E\|. \quad (2.10) $$

By using (2.6), (2.9), and the triangle inequality, we rearrange (2.10) as

$$ \|\hat{F}_i + \hat{F}_E - \hat{F}_{i+1}\| \le (1 + b_3\theta)\varepsilon(\|X - \hat{F}_i - \hat{E}\| + \|\hat{E} - \hat{F}_E\| + \|X\|) \quad (2.11) $$

$$ = (1 + b_3\theta)\varepsilon\|X\| + O(\varepsilon_\ell\varepsilon). \quad (2.12) $$

Using the triangle inequality together with (2.6), (2.9), and (2.12), we obtain

$$ \|X - \hat{F}_{i+1}\| = \|X - \hat{F}_i - \hat{E} + \hat{E} - \hat{F}_E + \hat{F}_i + \hat{F}_E - \hat{F}_{i+1}\| \quad (2.13) $$

$$ \le \|X - \hat{F}_i - \hat{E}\| + \|\hat{E} - \hat{F}_E\| + \|\hat{F}_i + \hat{F}_E - \hat{F}_{i+1}\| \quad (2.14) $$

$$ \le \phi\|X - \hat{F}_i\| + \xi\|X\| + O(\varepsilon_\ell\varepsilon), \quad (2.15) $$

with

$$ \phi = (1 + b_1\theta_\ell)\varepsilon_\ell + (b_2 + 1)\theta\varepsilon \quad (2.16) $$

and

$$ \xi = (1 + (b_2 + b_3 + 2)\theta)\varepsilon. \quad (2.17) $$

Noting that the first factors $\hat{F}_0$ computed at Line 1 satisfy by (2.1)

$$ \|X - \hat{F}_0\| \le (1 + b_1\theta_\ell)\varepsilon_\ell\|X\| \le \phi\|X\| \quad (2.18) $$

concludes the proof. □
In the rest of this section, we first perform an error analysis to determine the attainable accuracy and convergence rate of Algorithm 2.1. We then perform a complexity analysis to determine under which conditions the algorithm can be expected to be faster than a standard uniform precision LRA performed entirely in high precision u.

**2.1. Error analysis.** In order to carry out the error analysis of Algorithm 2.1, we will use the standard model of floating-point arithmetic [18, sect. 2.2]. In addition, we also need to make the following three assumptions on the numerical behavior of the LRA, DECOMPRESS, and RECOMPRESS kernels. First, we assume that computing $F = \text{LRA}(X, \varepsilon_\ell)$ in precision $u_\ell$ yields computed factors $\hat{F}$ satisfying

$$ \|X - \hat{F}\| \le (\varepsilon_\ell + b_1 u_\ell) \|X\| = (1 + b_1 \theta_\ell) \varepsilon_\ell \|X\|. \quad (2.1) $$

Second, we assume that computing $F = \text{DECOMPRESS}(F)$ in precision $u$ yields computed factors $\hat{F}$ satisfying

$$ \|F - \hat{F}\| \le b_2 u \|F\| = b_2 \theta \varepsilon \|F\|. \quad (2.2) $$

Third, we assume that computing $F = \text{RECOMPRESS}(F, \varepsilon)$ in precision $u$ yields computed factors $\hat{F}$ satisfying

$$ \|F - \hat{F}\| \le (\varepsilon + b_3 u) \|F\| = (1 + b_3 \theta) \varepsilon \|F\|. \quad (2.3) $$

In (2.1)-(2.3), the error bounds are parameterized by the constants $b_1, b_2$, and $b_3$, which will usually depend on the dimensions of $X$.

For simplicity, we ignore any rounding errors associated with the scaling by $\alpha = \|E\|$, which are negligible and can in fact be prevented by rounding $\alpha$ to the nearest power of two (in binary floating-point arithmetic).

We are now ready to prove the following result.

**THEOREM 2.1.** If Algorithm 2.1 is applied to $X$ with LRA, DECOMPRESS, and RECOMPRESS kernels satisfying (2.1)-(2.3), then after $k$ iterations, the computed factors $\hat{F}$ satisfy

$$ \|X - \hat{F}\| \le (\phi^{k+1} + \xi + O(\varepsilon_\ell \varepsilon)) \|X\|, \quad (2.4) $$

with $\phi = (1 + b_1 \theta_\ell) \varepsilon_\ell + (b_2 + 1) \theta \varepsilon$ and $\xi = (1 + (b_2 + b_3 + 2) \theta) \varepsilon$.

*Proof.* Defining $\hat{F}_i$ as the computed factors after $i$ iterations, our goal is to obtain a bound of the form

$$ \|X - \hat{F}_{i+1}\| \le \phi \|X - \hat{F}_i\| + \xi \|X\|, \quad \phi < 1, $$

which will allow us to conclude that the error contracts by a factor $\phi$ at each iteration, until it converges to its maximum attainable accuracy $\xi$.

Given $\hat{F}_i$, the first step is to decompress it and compute $E$ at Line 3; by (2.2), the computed $\hat{E}$ satisfies

$$ \|X - \hat{F}_i - \hat{E}\| \le \theta \varepsilon ((b_2 + 1) \|\hat{F}_i\| + \|X\|), \quad (2.5) $$

with an extra term $\theta \varepsilon (\|\hat{F}_i\| + \|X\|)$ on the right-hand side coming from the rounding error incurred by the subtraction. Using the triangle inequality $\|\hat{F}_i\| \le \|X - \hat{F}_i\| + \|X\|$, we can rearrange (2.5) as

$$ \|X - \hat{F}_i - \hat{E}\| \le \theta \varepsilon ((b_2 + 1) \|X - \hat{F}_i\| + (b_2 + 2) \|X\|). \quad (2.6) $$
**Theorem 2.1** states that the approximation error $\|X - \hat{F}\|$ contracts by a factor $\phi = O(\varepsilon_\ell)$ at each iteration, until it reaches its maximum attainable accuracy $\xi = O(\varepsilon)$. Thus, after $k$ iterations, the error is of order $\varepsilon_\ell^{k+1} + \varepsilon$, which means that we can actually estimate in advance approximately how many iterations are needed to achieve the desired accuracy $\varepsilon$ (up to constants):

$$n_{it} = \log(\varepsilon)/\log(\varepsilon_\ell) - 1. \qquad (2.19)$$

It is worth noting that, unlike iterative refinement for linear systems, there is no dependence on the condition number of $X$ and thus the only condition for Algorithm 2.1 to converge is that $\phi < 1$, which should always be the case as long as $\theta_\ell$ is sufficiently small. Therefore, the algorithm is extremely general, can be applied to any matrix or tensor of low numerical rank, and can make use of potentially very low precisions.

From a numerical perspective, Algorithm 2.1 is therefore quite appealing. We next discuss under which conditions it is also attractive from a computational perspective.

**2.2. Complexity analysis.** The cost of Algorithm 2.1 will mainly depend on two factors: the relative speed for computing in different precisions on the target hardware, and the numerical ranks of the objects encountered during the iterations $(X, E, \text{and } F)$.

We begin by bounding these ranks solely as a function of the numerical ranks of $X$ at given accuracies. In order to do so, we will need to bound the numerical rank of $A+B$ as a function of that of $A$ and $B$ by using the following lemma. In what follows, recall that when $X$ is a tensor, $\text{RANK}(X, \varepsilon)$ is a vector (see subsection 1.1) and all the (in)equalities involving this quantity should be interpreted componentwise.

LEMMA 2.2.

$$\text{RANK}(A+B, \varepsilon) \le \text{RANK}(A, \frac{\varepsilon\|A+B\|}{2\|A\|}) + \text{RANK}(B, \frac{\varepsilon\|A+B\|}{2\|B\|}). \qquad (2.20)$$

*Proof.* For the purpose of this proof only, we introduce an alternative definition of numerical rank that measures accuracy in an absolute sense, rather than a relative one. Let $\text{RANKABS}(X, \varepsilon)$ be the smallest integer $r_\varepsilon$ such that there exists $F$ of rank $r_\varepsilon$ satisfying $\|X - F\| \le \varepsilon$. The $\text{RANKABS}$ operator satisfies

$$\text{RANKABS}(X, \varepsilon\|X\|) = \text{RANK}(X, \varepsilon)$$

and

$$\text{RANKABS}(A+B, \varepsilon) \le \text{RANKABS}(A, \varepsilon/2) + \text{RANKABS}(B, \varepsilon/2).$$

Therefore, we have that

$$\begin{align*}
\text{RANK}(A+B, \varepsilon) &= \text{RANKABS}(A+B, \varepsilon\|A+B\|) \\
&\le \text{RANKABS}(A, \varepsilon\|A+B\|/2) + \text{RANKABS}(B, \varepsilon\|A+B\|/2) \\
&= \text{RANK}(A, \varepsilon\|A+B\|/2\|A\|) + \text{RANK}(B, \varepsilon\|A+B\|/2\|B\|). \quad \square
\end{align*}$$

Let us consider the $i$th iteration of Algorithm 2.1. Using (2.20), we can bound the rank of $E = X - F$ as

$$\text{RANK}(E, \varepsilon_\ell) \le \text{RANK}(X, \varepsilon_\ell\|E\|/2\|X\|) + \text{RANK}(F). \qquad (2.21)$$
Then we compute LRA($\hat{E}$, $\varepsilon_\ell$) at Line 7; by (2.1) the computed $\hat{F}_E$ satisfies

$$ \|\hat{E} - \hat{F}_E\| \le (1 + b_1 \theta_\ell) \varepsilon_\ell \|\hat{E}\|. \quad (2.7) $$

By using (2.6) and the triangle inequality, we rearrange (2.7) as

$$ \|\hat{E} - \hat{F}_E\| \le (1 + b_1 \theta_\ell) \varepsilon_\ell (\|X - \hat{F}_i\| + \|X - \hat{F}_i - \hat{E}\|) \quad (2.8) $$

$$ = (1 + b_1 \theta_\ell) \varepsilon_\ell \|X - \hat{F}_i\| + O(\varepsilon_\ell \varepsilon), \quad (2.9) $$

where we do not keep track explicitly of high order terms in $O(\varepsilon_\ell \varepsilon)$ for the sake of readability. Finally, we obtain the next iterate $F_{i+1}$ by recompressing $\hat{F}_i + \hat{F}_E$ at Line 9; by (2.3), the computed $\hat{F}_{i+1}$ satisfies

$$ \|\hat{F}_i + \hat{F}_E - \hat{F}_{i+1}\| \le (1 + b_3 \theta) \varepsilon \|\hat{F}_i + \hat{F}_E\|. \quad (2.10) $$

By using (2.6), (2.9), and the triangle inequality, we rearrange (2.10) as

$$ \|\hat{F}_i + \hat{F}_E - \hat{F}_{i+1}\| \le (1 + b_3 \theta) \varepsilon (\|X - \hat{F}_i - \hat{E}\| + \|\hat{E} - \hat{F}_E\| + \|X\|) \quad (2.11) $$

$$ = (1 + b_3 \theta) \varepsilon \|X\| + O(\varepsilon_\ell \varepsilon). \quad (2.12) $$

Using the triangle inequality together with (2.6), (2.9), and (2.12), we obtain

$$ \|X - \hat{F}_{i+1}\| = \|X - \hat{F}_i - \hat{E} + \hat{E} - \hat{F}_E + \hat{F}_i + \hat{F}_E - \hat{F}_{i+1}\| \quad (2.13) $$

$$ \le \|X - \hat{F}_i - \hat{E}\| + \|\hat{E} - \hat{F}_E\| + \|\hat{F}_i + \hat{F}_E - \hat{F}_{i+1}\| \quad (2.14) $$

$$ \le \phi \|X - \hat{F}_i\| + \xi \|X\| + O(\varepsilon_\ell \varepsilon), \quad (2.15) $$

with

$$ \phi = (1 + b_1 \theta_\ell) \varepsilon_\ell + (b_2 + 1) \theta \varepsilon \quad (2.16) $$

and

$$ \xi = (1 + (b_2 + b_3 + 2)\theta)\varepsilon. \quad (2.17) $$

Noting that the first factors $\hat{F}_0$ computed at Line 1 satisfy by (2.1)

$$ \|X - \hat{F}_0\| \le (1 + b_1 \theta_\ell) \varepsilon_\ell \|X\| \le \phi \|X\| \quad (2.18) $$

concludes the proof. □

Theorem 2.1 states that the approximation error $\|X - \hat{F}\|$ contracts by a factor $\phi = O(\varepsilon_\ell)$ at each iteration, until it reaches its maximum attainable accuracy $\xi = O(\varepsilon)$. Thus, after $k$ iterations, the error is of order $\varepsilon_\ell^{k+1} + \varepsilon$, which means that we can actually estimate in advance approximately how many iterations are needed to achieve the desired accuracy $\varepsilon$ (up to constants):

$$ n_{it} = \log(\varepsilon)/\log(\varepsilon_\ell) - 1. \quad (2.19) $$

It is worth noting that, unlike iterative refinement for linear systems, there is no dependence on the condition number of $X$ and thus the only condition for Algorithm 2.1 to converge is that $\phi < 1$, which should always be the case as long as $\theta_\ell$ is sufficiently small. This is explained by the fact that the speed of convergence of Algorithm 2.1 depends on the backward error $\|A - XY^T\|$, rather than the forward one. Therefore, Algorithm 2.1 is extremely general: it can be applied to any matrix or tensor of low numerical rank and it can make use of potentially very low precisions.

From a numerical perspective, Algorithm 2.1 is therefore quite appealing. We next discuss under which conditions it is also attractive from a computational perspective.
Since $F$ is the low-rank factorization of $X$ after $i$ iterations, it achieves an accuracy of $\varepsilon_{\ell}^i$ and so

$$ \mathrm{RANK}(F) = \mathrm{RANK}(X, \varepsilon_{\ell}^i) \qquad (2.22) $$

(we assume here that $\varepsilon_{\ell}^i \le \varepsilon$ as otherwise the method would be stopped). After $i$ iterations, $\|E\| \le \varepsilon_{\ell}^i \|X\|$ and so overall (2.21) becomes

$$ \mathrm{RANK}(E, \varepsilon_{\ell}) \le \mathrm{RANK}(X, \varepsilon_{\ell}^{i+1}/2) + \mathrm{RANK}(X, \varepsilon_{\ell}^{i}). \qquad (2.23) $$

Since $\mathrm{RANK}(X, \varepsilon') \le \mathrm{RANK}(X, \varepsilon)$ for $\varepsilon \le \varepsilon'$, assuming that $\varepsilon \le \varepsilon_{\ell}^{i+1}/2$, we can obtain a simpler but weaker version of (2.23),

$$ \mathrm{RANK}(E, \varepsilon_{\ell}) \le 2 \mathrm{RANK}(X, \varepsilon). \qquad (2.24) $$

Thus, the rank of $E$ at any iteration is at most twice as large as the numerical rank of $X$ at accuracy $\varepsilon$.

With a similar argument, the rank of the refined factors $F + F_E$ is bounded by

$$ \mathrm{RANK}(F + F_E) \le \mathrm{RANK}(F) + \mathrm{RANK}(F_E) \qquad (2.25) $$

$$ = \mathrm{RANK}(F) + \mathrm{RANK}(E, \varepsilon_{\ell}) \qquad (2.26) $$

$$ \le \mathrm{RANK}(X, \varepsilon_{\ell}^{i+1}/2) + 2 \mathrm{RANK}(X, \varepsilon_{\ell}^{i}), \qquad (2.27) $$

where (2.27) follows from (2.22) and (2.23). Again, a simpler but weaker bound is

$$ \mathrm{RANK}(F + F_E) \le 3 \mathrm{RANK}(X, \varepsilon), \qquad (2.28) $$

showing that at any iteration the rank of the factors before recompression is at most three times the numerical rank of $X$ at accuracy $\varepsilon$.

Now that we have bounded the ranks of the objects appearing in Algorithm 2.1, we are ready to analyze its cost. We denote as $p$ and $s$ the product and the sum of the dimensions of $X$, respectively. Thus, if $X \in \mathbb{R}^{m \times n}$ is a matrix, $p = mn$ and $s = m + n$; if $X \in \mathbb{R}^{n_1 \times \cdots \times n_d}$ is a tensor of order $d$, $p = \prod n_i$ and $s = \sum n_i$. For readability, we only report the dominant term for the cost of each assumption/line, and assume a large scale setting, so that $p \gg s$.

We make the following assumptions on the flops cost of the different kernels:

$$ \text{Flops LRA}(X, \varepsilon) = c_1 p \mathrm{RANK}(X, \varepsilon)_1 + o(p), \qquad (2.29a) $$

$$ \text{Flops DECOMPRESS}(F) = c_2 p \mathrm{RANK}(F)_1 + o(p), \qquad (2.29b) $$

$$ \text{Flops RECOMPRESS}(F, \varepsilon) = o(p). \qquad (2.29c) $$

For tensors, the cost of LRA and DECOMPRESS depends on the order in which the dimensions are treated (see section 4). Here we assume they are treated in the natural order, without loss of generality since the dimensions can be arbitrarily reordered. In this case, the dominant cost of LRA and DECOMPRESS is proportional to the first coefficient of the rank vectors, denoted as $\mathrm{RANK}(X, \varepsilon)_1$. Note that for matrices, the rank is a scalar so that $\mathrm{RANK}(X, \varepsilon)_1 = \mathrm{RANK}(X, \varepsilon)$.

Let us now measure the flops cost of the $i$th iteration of Algorithm 2.1.

* Line 3: $c_2 p \mathrm{RANK}(X, \varepsilon_{\ell}^i)_1$ flops.

* Line 4: $2p$ flops (for the Frobenius norm).

* Line 6: $p$ flops.

* Line 7: $c_1 p (\mathrm{RANK}(X, \varepsilon_{\ell}^{i+1}/2)_1 + \mathrm{RANK}(X, \varepsilon_{\ell}^i)_1)$ flops.
**2.2. Complexity analysis.** The cost of Algorithm 2.1 will mainly depend on two factors: the relative speed for computing in different precisions on the target hardware, and the numerical ranks of the objects encountered during the iterations $(X, E, \text{and } F)$.

We begin by bounding these ranks solely as a function of the numerical ranks of $X$ at given accuracies. In order to do so, we will need to bound the numerical rank of $A+B$ as a function of that of $A$ and $B$ by using the following lemma. In what follows, recall that when $X$ is a tensor, $\text{RANK}(X, \varepsilon)$ is a vector (see subsection 1.1) and all the (in)equalities involving this quantity should be interpreted componentwise.

LEMMA 2.2.

$$ \text{RANK}(A+B, \varepsilon) \le \text{RANK}(A, \frac{\varepsilon\|A+B\|}{2\|A\|}) + \text{RANK}(B, \frac{\varepsilon\|A+B\|}{2\|B\|}). \quad (2.20) $$

*Proof.* For the purpose of this proof only, we introduce an alternative definition of numerical rank that measures accuracy in an absolute sense, rather than a relative one. Let $\text{RANKABS}(X, \varepsilon)$ be the smallest integer $r_\varepsilon$ such that there exists $F$ of rank $r_\varepsilon$ satisfying $\|X - F\| \le \varepsilon$. The $\text{RANKABS}$ operator satisfies

$$ \text{RANKABS}(X, \varepsilon\|X\|) = \text{RANK}(X, \varepsilon) $$

and

$$ \text{RANKABS}(A+B, \varepsilon) \le \text{RANKABS}(A, \varepsilon/2) + \text{RANKABS}(B, \varepsilon/2) $$

by the triangle inequality. Therefore, we have that

$$ \begin{aligned} \text{RANK}(A+B, \varepsilon) &= \text{RANKABS}(A+B, \varepsilon\|A+B\|) \\ &\le \text{RANKABS}(A, \varepsilon\|A+B\|/2) + \text{RANKABS}(B, \varepsilon\|A+B\|/2) \\ &= \text{RANK}(A, \varepsilon\|A+B\|/2\|A\|) + \text{RANK}(B, \varepsilon\|A+B\|/2\|B\|). \quad \square \end{aligned} $$

Let us consider the $i$th iteration of Algorithm 2.1. Using (2.20), we can bound the rank of $E = X - F$ as

$$ \text{RANK}(E, \varepsilon_i) \le \text{RANK}(X, \varepsilon_i\|E\|/2\|X\|) + \text{RANK}(F). \quad (2.21) $$

Since $F$ is the low-rank factorization of $X$ after $i$ iterations, it achieves an accuracy of $\varepsilon_i^i$ and so

$$ \text{RANK}(F) = \text{RANK}(X, \varepsilon_i^i) \quad (2.22) $$

(we assume here that $\varepsilon_i^i \le \varepsilon$ as otherwise the method would be stopped). After $i$ iterations, $\|E\| \le \varepsilon_i^i \|X\|$ and so overall (2.21) becomes

$$ \text{RANK}(E, \varepsilon_i) \le \text{RANK}(X, \varepsilon_i^{i+1}/2) + \text{RANK}(X, \varepsilon_i^i). \quad (2.23) $$

Since $\text{RANK}(X, \varepsilon') \le \text{RANK}(X, \varepsilon)$ for $\varepsilon \le \varepsilon'$, assuming that $\varepsilon \le \varepsilon_i^{i+1}/2$, we can obtain a simpler but weaker version of (2.23),

$$ \text{RANK}(E, \varepsilon_i) \le 2 \text{RANK}(X, \varepsilon). \quad (2.24) $$

Thus, the rank of $E$ at any iteration is at most twice as large as the numerical rank of $X$ at accuracy $\varepsilon$.
• Line 8: $o(p)$ flops.

• Line 9: $o(p)$ flops.

The dominant steps are the calls to LRA (Line 7), which is performed in low precision, and to DECOMPRESS (Line 3), which is performed in high precision. The computation of $\alpha$ and the scaling by $\alpha$ (Lines 4 and 6) could also be significant if the ranks are very small.

Summing the costs of these dominant steps across all $n_{it}$ iterations, plus the cost of the initial LRA (Line 1), gives a total flops cost of

$$ \text{Flops IR} = c_1 p \text{RANK}(X, \varepsilon_\ell)_1 + p \sum_{i=1}^{n_{it}} \left( c_1 \text{RANK}(X, \varepsilon_\ell^{i+1}/2)_1 + c_1 \text{RANK}(X, \varepsilon_\ell^i)_1 + c_2 \text{RANK}(X, \varepsilon_\ell^i)_1 + 3 \right). \quad (2.30) $$

Since some of these flops are performed in low precision and some in high precision, we must account for the different speeds of different arithmetics. To do so, we ponderate the low precision flops by a weight $\omega_\ell < 1$ which indicates the relative performance ratio between the low and the high precision. We obtain

$$ \text{Time IR} \propto \omega_\ell c_1 p \text{RANK}(X, \varepsilon_\ell)_1 + p \sum_{i=1}^{n_{it}} \left( \omega_\ell c_1 \text{RANK}(X, \varepsilon_\ell^{i+1}/2)_1 + \omega_\ell c_1 \text{RANK}(X, \varepsilon_\ell^i)_1 + c_2 \text{RANK}(X, \varepsilon_\ell^i)_1 + 3 \right), \quad (2.31) $$

where the “Time” formula should be taken as a rough estimator of the time performance of the algorithm, although in practice the actual execution time depends on a number of other complex factors such as the arithmetic intensity of the operations, the data locality, and the parallelism.

This cost is to be compared with the cost of simply computing LRA($X, \varepsilon$) in the high precision $u$, given by

$$ \text{Time Ref.} \propto c_1 p \text{RANK}(X, \varepsilon)_1. \quad (2.32) $$

This complexity analysis reveals two situations where Algorithm 2.1 can outperform the uniform precision approach.

Numerical ranks $\text{RANK}(X, \varepsilon_\ell^i)$ rapidly decreasing as $\varepsilon_\ell$ increases. The first situation is when the numerical ranks of $X$ at accuracy lower than the final target $\varepsilon$ are much smaller than $\text{RANK}(X, \varepsilon)$. This can certainly be the case in some applications. For example, in the case of matrices, the numerical rank of $X$ at any given accuracy is determined by its singular values. If the singular values decay rapidly, $\text{RANK}(X, \varepsilon_\ell)$ will in general be significantly smaller than $\text{RANK}(X, \varepsilon)$. The most extreme example is a matrix with one large singular value and all the remaining $\text{RANK}(X, \varepsilon) - 1$ singular values just above $\varepsilon\|X\|$. In this extreme case, all the iterations except the last will only need to compute on rank-1 matrices, so the overall cost of the method will be dominated by that of the last iteration, which is roughly $\omega_\ell c_1 p \text{RANK}(X, \varepsilon)$, and so is always smaller than (2.32). In a more realistic setting where the singular values decay exponentially, we may still expect (2.31) to be smaller than (2.32) even for reasonably traditional values of $\omega_\ell$.

This situation also extends to tensors, although we do not have such a simple characterization as one based on singular values. Essentially, if the underlying matrices used in the tensor representation exhibit rapidly decaying singular values, then the rank vectors $\text{RANK}(X, \varepsilon_\ell)$ will take much smaller values than $\text{RANK}(X, \varepsilon)$ when $\varepsilon_\ell \gg \varepsilon$.
With a similar argument, the rank of the refined factors $F + F_E$ is bounded by

$$ \mathrm{RANK}(F + F_E) \leq \mathrm{RANK}(F) + \mathrm{RANK}(F_E) \quad (2.25) $$

$$ = \mathrm{RANK}(F) + \mathrm{RANK}(E, \varepsilon_{\ell}) \quad (2.26) $$

$$ \leq \mathrm{RANK}(X, \varepsilon_{\ell}^{i+1}/2) + 2 \mathrm{RANK}(X, \varepsilon_{\ell}^{i}), \quad (2.27) $$

where (2.27) follows from (2.22) and (2.23). Again, a simpler but weaker bound is

$$ \mathrm{RANK}(F + F_E) \leq 3 \mathrm{RANK}(X, \varepsilon), \quad (2.28) $$

showing that at any iteration the rank of the factors before recompression is at most three times the numerical rank of $X$ at accuracy $\varepsilon$.

Now that we have bounded the ranks of the objects appearing in Algorithm 2.1, we are ready to analyze its cost. We denote as $p$ and $s$ the product and the sum of the dimensions of $X$, respectively. Thus, if $X \in \mathbb{R}^{m \times n}$ is a matrix, $p = mn$ and $s = m + n$; if $X \in \mathbb{R}^{n_1 \times \cdots \times n_d}$ is a tensor of order $d$, $p = \prod n_i$ and $s = \sum n_i$. For readability, we only report the dominant term for the cost of each assumption/line, and assume a large scale setting, so that $p \gg s$.

We make the following assumptions on the flops (floating-point operations) cost of the different kernels:

$$ \mathrm{Flops\ LRA}(X, \varepsilon) = c_1 p \mathrm{RANK}(X, \varepsilon)_1 + o(p), \quad (2.29a) $$

$$ \mathrm{Flops\ DECOMPRESS}(F) = c_2 p \mathrm{RANK}(F)_1 + o(p), \quad (2.29b) $$

$$ \mathrm{Flops\ RECOMPRESS}(F, \varepsilon) = o(p). \quad (2.29c) $$

For tensors, the cost of LRA and DECOMPRESS depends on the order in which the dimensions are treated (see section 4). Here we assume they are treated in the natural order, without loss of generality since the dimensions can be arbitrarily reordered. In this case, the dominant cost of LRA and DECOMPRESS is proportional to the first coefficient of the rank vectors, denoted as $\mathrm{RANK}(X, \varepsilon)_1$. Note that for matrices, the rank is a scalar so that $\mathrm{RANK}(X, \varepsilon)_1 = \mathrm{RANK}(X, \varepsilon)$.

Let us now measure the flops cost of the $i$th iteration of Algorithm 2.1.

* Line 3: $c_2 p \mathrm{RANK}(X, \varepsilon_{\ell}^i)_1$ flops.

* Line 4: $2p$ flops (for the Frobenius norm).

* Line 6: $p$ flops.

* Line 7: $c_1 p (\mathrm{RANK}(X, \varepsilon_{\ell}^{i+1}/2)_1 + \mathrm{RANK}(X, \varepsilon_{\ell}^i)_1)$ flops.

* Line 8: $o(p)$ flops.

* Line 9: $o(p)$ flops.

The dominant steps are the calls to LRA (Line 7), which is performed in low precision, and to DECOMPRESS (Line 3), which is performed in high precision. The computation of $\alpha$ and the scaling by $\alpha$ (Lines 4 and 6) could also be significant if the ranks are very small.

Summing the costs of these dominant steps across all $n_{it}$ iterations, plus the cost of the initial LRA (Line 1), gives a total flops cost of

$$ \mathrm{Flops\ IR} = c_1 p \mathrm{RANK}(X, \varepsilon_{\ell})_1 + p \sum_{i=1}^{n_{it}} \left( c_1 \mathrm{RANK}(X, \varepsilon_{\ell}^{i+1}/2)_1 + c_2 \mathrm{RANK}(X, \varepsilon_{\ell}^i)_1 + 3 \right). \quad (2.30) $$

Since some of these flops are performed in low precision and some in high precision, we must account for the different speeds of different arithmetics. To do so, we ponderate
Very fast low precision arithmetic (or very slow high precision arithmetic). The second situation is when the low precision arithmetic is much faster than the high precision one, that is, when $\omega_{\ell} \ll 1$. This is becoming increasingly common for low precisions on modern hardware, especially accelerators. For example, the fp16 and bfloat16 arithmetics can be up to 16 times faster than fp32 arithmetic on recent NVIDIA GPUs. Similarly, fp16 arithmetic can be up to 8 times faster than fp32 on the AMD Instinct MI250X GPUs. In this situation, the time cost (2.31) can indeed be lower than (2.32), even in the worst case scenario where the ranks of all objects throughout the iterations attain their upper bound of $2 \text{RANK}(X, \varepsilon)$ for E. Indeed, in this case, (2.31) becomes

$$p \text{RANK}(X, \varepsilon)_1 (\omega_{\ell} c_1 (2n_{\text{it}} + 1) + n_{\text{it}} c_2) + 3p. \quad (2.33)$$

Neglecting the $3p$ term, (2.33) is smaller than (2.32) if

$$\omega_{\ell}(2n_{\text{it}} + 1) + n_{\text{it}}c_2/c_1 \le 1. \quad (2.34)$$

For most LRA algorithms of interest, the cost of compressing (LRA kernel) a full object is higher than the cost of decompressing it (DECOMPRESS kernel), that is, $c_1 > c_2$. Therefore (2.34) can certainly be satisfied for sufficiently small $\omega_{\ell}$.

Moreover, we now describe a case where the DECOMPRESS kernel can also be performed in low precision. Indeed, while Algorithm 2.1 requires the DECOMPRESS kernel Line 3 to be performed in high precision, this kernel has the particularity of taking as input the factors F, which are stored in low precision. Therefore, what we really need is to compute in high precision with low precision numbers; this happens to be an easier task than the more general problem of computing in high precision with high precision numbers. In fact, an increasing range of modern hardware provides the capability to perform this task inexpensively. For example, the NVIDIA GPUs are equipped with so-called tensor core units that can carry out matrix multiplication with half precision (16-bit) matrices at the accuracy of single precision (fp32) arithmetic but at the much higher speed of half precision arithmetic. Similar instructions are available on several other architectures, which have been analyzed by Blanchard et al. [4] under a common framework called block FMA.

In our context, assuming the DECOMPRESS kernel can make use of these block FMA units changes the cost analysis quite significantly since all the dominant operations are then performed at the speed of the low precision. Even under the same worst case scenario where the ranks of all objects throughout the iterations attain their upper bound, (2.34) now becomes

$$\omega_{\ell}(2n_{\text{it}} + 1 + n_{\text{it}}c_2/c_1) \le 1. \quad (2.35)$$

Thus, if block FMA hardware can be exploited, even the ratio $c_2/c_1$ need not be necessarily small for Algorithm 2.1 to be effective.

Clearly, the two situations discussed above are not exclusive, so in practice it is very possible that we have both smaller ranks and a fast low precision. To further assess under which condition Algorithm 2.1 can outperform the standard approach, we now specialize it to specific matrix or tensor algorithms in section 3 and section 4.

**3. Application to matrix low-lank approximation.** In this section we specialize Algorithm 2.1 to the case where $X \in \mathbb{R}^{m \times n}$ is a matrix.

For any unitarily invariant norm, the optimal LRA algorithm for matrices is to compute their singular value decomposition (SVD) and truncate it to the desired
the low precision flops by a weight $\omega_{\ell} < 1$ which indicates the relative performance ratio between the low and the high precision. We obtain

$$ \text{Time IR} \propto \omega_{\ell} c_1 p \text{RANK}(X, \varepsilon_{\ell})_1 + p \sum_{i=1}^{n_{\text{it}}} \left( \omega_{\ell} c_1 \text{RANK}(X, \varepsilon_{\ell}^{i+1}/2)_1 + \omega_{\ell} c_1 \text{RANK}(X, \varepsilon_{\ell}^{i})_1 + c_2 \text{RANK}(X, \varepsilon_{\ell}^{i})_1 + 3 \right), \quad (2.31) $$

where the “Time” formula should be taken as a rough estimator of the time performance of the algorithm, although in practice the actual execution time depends on a number of other complex factors such as the arithmetic intensity of the operations, the data locality, and the parallelism.

This cost is to be compared with the cost of simply computing LRA($X, \varepsilon$) in the high precision $u$, given by

$$ \text{Time Ref.} \propto c_1 p \text{RANK}(X, \varepsilon)_1. \quad (2.32) $$

This complexity analysis reveals two situations where Algorithm 2.1 can outperform the uniform precision approach.

Numerical ranks RANK($X, \varepsilon_{\ell}^i$) rapidly decreasing as $\varepsilon_{\ell}$ increases. The first situation is when the numerical ranks of $X$ at accuracy lower than the final target $\varepsilon$ are much smaller than RANK($X, \varepsilon$). This can certainly be the case in some applications. For example, in the case of matrices, the numerical rank of $X$ at any given accuracy is determined by its singular values: $r = \text{RANK}(X, \varepsilon_{\ell})$ is the number of singular values that need to be kept for the truncated SVD to approximate $X$ with accuracy at least $\varepsilon_{\ell}$. Thus, if the singular values decay rapidly, RANK($X, \varepsilon_{\ell}$) will in general be significantly smaller than RANK($X, \varepsilon$). The most extreme example is a matrix with one large singular value and all the remaining RANK($X, \varepsilon$) - 1 singular values just above $\varepsilon \|X\|$. In this extreme case, all the iterations except the last will only need to compute on rank-1 matrices, so the overall cost of the method will be dominated by that of the last iteration, which is roughly $\omega_{\ell} c_1 p \text{RANK}(X, \varepsilon)$, and so is always smaller than (2.32). In a more realistic setting where the singular values decay exponentially, we may still expect (2.31) to be smaller than (2.32) even for reasonably traditional values of $\omega_{\ell}$.

This situation also extends to tensors, although we do not have such a simple characterization as one based on singular values. Essentially, if the underlying matrices used in the tensor representation exhibit rapidly decaying singular values, then the rank vectors RANK($X, \varepsilon_{\ell}$) will take much smaller values than RANK($X, \varepsilon$) when $\varepsilon_{\ell} \gg \varepsilon$.

Very fast low precision arithmetic (or very slow high precision arithmetic). The second situation is when the low precision arithmetic is much faster than the high precision one, that is, when $\omega_{\ell} \ll 1$. This is becoming increasingly common for low precisions on modern hardware, especially accelerators. For example, the fp16 and bfloat16 arithmetics can be up to 16 times faster than fp32 arithmetic on recent NVIDIA GPUs. Similarly, fp16 arithmetic can be up to 8 times faster than fp32 on the AMD Instinct MI250X GPUs.

Moreover, such hardware provides another crucial feature which allows the DECOMPRESS kernel to also be performed in low precision. Indeed, while Algorithm 2.1 requires the DECOMPRESS kernel Line 3 to be performed in high precision, this kernel has the particularity of taking as input the factors $F$, which are stored in low precision. Therefore, what we really need is to compute in high precision with low
accuracy, a result known as Eckart-Young theorem [10]. Specifically, given the SVD $X = U\Sigma V^T$, the optimal rank-k approximation of X is

$$ \arg \min_{\substack{M \in \mathbb{R}^{m \times n} \\ \mathrm{RANK}(M)=k}} \|X - M\| = U_k \Sigma_k V_k^T, $$

where $U_k \Sigma_k V_k^T$ is the truncated SVD of X, formed by the first k singular vectors and values only.

While the truncated SVD provides the best accuracy, it is expensive due to the need to compute the full SVD before truncating it: it requires $O(mn^2)$ flops. For this reason, low-rank matrix approximations are often computed using slightly less accurate but cheaper alternatives. In this paper, we focus on two widely popular choices: the truncated QR decomposition with column pivoting (QRCP, subsection 3.1), and the randomized SVD (subsection 3.2). For each algorithm, we discuss their use as LRA kernel in our mixed precision iterative refinement approach (Algorithm 2.1): in particular, we check that the algorithm satisfies the assumptions (2.1)-(2.3) of the error analysis and the assumptions (2.29a)-(2.29c) of the complexity analysis. We also explain how to perform the RECOMPRESS kernel based on the specific LRA choice.

**3.1. Truncated QRCP decomposition.** Our first choice of LRA algorithm is a truncated QRCP decomposition. QRCP decomposes the original matrix $X \in \mathbb{R}^{m \times n}$ as

$$ XP = QR $$

where $Q \in \mathbb{R}^{m \times n}$ is a matrix with orthonormal columns, $R \in \mathbb{R}^{n \times n}$ is an upper triangular matrix, and $P \in \mathbb{R}^{n \times n}$ is a permutation matrix. Thanks to pivoting, the QRCP decomposition can be used as a rank-revealing algorithm because the norm of the trailing submatrix at each step k of the QR factorization, $\|XP - Q_k R_k\|$, is monotonically decreasing. Therefore, we can stop the QR factorization as soon as this norm becomes smaller than the target tolerance $\epsilon_\ell$. We obtain

$$ X \approx Q_k V_k^T, \quad Q_k \in \mathbb{R}^{m \times k}, \quad V_k = P R_k^T \in \mathbb{R}^{n \times k}, \qquad (3.1) $$

where $k = \mathrm{RANK}(X, \epsilon_\ell)$. To avoid the need to apply the permutation P each time we need to apply X, we form and store $V_k = P R_k^T$.

We can implement a suitable RECOMPRESS algorithm using this truncated QRCP decomposition. Several versions are possible; we describe in Algorithm 3.1 the one that we will use in this article. Given a (non-optimal) truncated QRCP decomposition $Q_{in}V_{in}^T$, Algorithm 3.1 recompresses it into an optimal LRA as follows. First, a thin QR factorization $QR$ of $V_{in}$ is computed. Then, the product $Q_{in}R^T$ is formed and a truncated QRCP decomposition $Q_{out}\bar{V}^T$ is computed at the desired target accuracy $\epsilon$. This yields the recompressed left factor $Q_{out}$; the recompressed right factor $V_{out}$ is obtained by forming $\bar{Q}\bar{V}$.

We now check that the assumptions (2.1)-(2.3) of the error analysis are satisfied for truncated QRCP and determine the value of the constants in the error bounds. For (2.1), we can analyze the truncated QRCP by separating the truncation and rounding errors. As explained above, we can control the size of the truncation error by stopping the QRCP decomposition once the approximation error falls below the desired threshold $\epsilon_\ell$, so that in exact arithmetic we obtain

$$ Q_k R_k = XP + E, \quad \|E\| \le \epsilon_\ell \|X\|. $$
precision numbers; this happens to be an easier task than the more general problem of computing in high precision with high precision numbers. In fact, an increasing range of modern hardware provides the capability to perform this task inexpensively. For example, the NVIDIA GPUs are equipped with so-called tensor core units that can carry out matrix multiplication with half precision (16-bit) matrices at the accuracy of single precision (fp32) arithmetic but at the much higher speed of half precision arithmetic. Similar instructions are available on several other architectures, which have been analyzed by Blanchard et al. [6] under a common framework called block FMA. In our context, a block FMA unit therefore allows all the dominant operations to be performed at the speed of the low precision, which is much higher than that of the high precision. In this situation, the time cost of Algorithm 2.1 can be much lower than (2.32), even in the worst case scenario where the ranks of all objects throughout the iterations attain their upper bound of 2 RANK(X, ε) for E. Indeed, in this case, we have

$$ \text{Time IR} \propto \omega_{\ell} p \text{RANK}(X, \varepsilon)_1 (c_1(2n_{\text{it}} + 1) + n_{\text{it}}c_2) + 3p. \quad (2.33) $$

Neglecting the 3p term, (2.33) is smaller than (2.32) if

$$ \omega_{\ell}(2n_{\text{it}} + 1 + n_{\text{it}}c_2/c_1) \le 1. \quad (2.34) $$

For many LRA algorithms, including all those considered in this article, the cost of compressing (LRA kernel) a full object is higher than the cost of decompressing it (DECOMPRESS kernel), that is, $c_1 > c_2$. Therefore (2.34) is certainly satisfied if $\omega_{\ell}(3n_{\text{it}} + 1) \le 1$.

To illustrate this cost analysis, let us take the NVIDIA GPU tensor cores as an example, for which fp16 arithmetic ($u_{\ell} \approx 3 \times 10^{-4}$) is up to 16 times faster ($\omega_{\ell} = 1/16$) than fp32 arithmetic ($u \approx 6 \times 10^{-8}$). Since $u_{\ell}^2 \approx u$, a single refinement step ($n_{\text{it}} = 1$) suffices to recover fp32 accuracy. Thus, the term $\omega_{\ell}(3n_{\text{it}} + 1)$ is equal to $4/16 = 1/4$, which suggests that a speedup of up to a factor 4× could be expected.

Clearly, the two situations discussed above are not exclusive, so in practice it is very possible that we have both smaller ranks and a fast low precision. To further assess under which condition Algorithm 2.1 can outperform the standard approach, we now specialize it to specific matrix or tensor algorithms in section 3 and section 4.

**3. Application to matrix low-lank approximation.** In this section we specialize Algorithm 2.1 to the case where $X \in \mathbb{R}^{m \times n}$ is a matrix.

For any unitarily invariant norm, the optimal LRA algorithm for matrices is to compute their singular value decomposition (SVD) and truncate it to the desired accuracy, a result known as Eckart-Young theorem [12]. Specifically, given the SVD $X = U \Sigma V^T$, the optimal rank-k approximation of X is

$$ \arg \min_{M \in \mathbb{R}^{m \times n} \atop \text{RANK}(M)=k} \|X - M\| = U_k \Sigma_k V_k^T, $$

where $U_k \Sigma_k V_k^T$ is the truncated SVD of X, formed by the first k singular vectors and values only.

While the truncated SVD provides the best accuracy, it is expensive due to the need to compute the full SVD before truncating it: it requires $O(mn^2)$ flops. For this reason, low-rank matrix approximations are often computed using slightly less accurate but cheaper alternatives. In this paper, we focus on two widely popular choices: the truncated QR decomposition with column pivoting (QRCP, subsection 3.1), and the randomized SVD (subsection 3.2). For each algorithm, we discuss their use as
**Algorithm 3.1 RECOMPRESS algorithm using truncated QRCP decomposition.**

**Input:** a truncated QRCP decomposition $Q_{in}V_{in}^T$ of the form (3.1);
$\varepsilon$, the target accuracy.

**Output:** a recompressed QRCP decomposition $Q_{out}V_{out}^T$.

1. Compute the thin QR factorization $\bar{Q}R = V_{in}$.

2. Compute the truncated QRCP decomposition $Q_{out}\bar{V}^T = Q_{in}R^T$.

3. $V_{out} \leftarrow \bar{Q}\bar{V}$

To account for the rounding errors, we use standard analysis of QR decomposition [16, p. 361], which shows that the computed QR factors satisfy the backward error bound

$$ \hat{Q}_k \hat{R}_k = XP + \Delta X + E, \quad \|\Delta X\| \le \sqrt{k} \tilde{\gamma}_{mk} \|X\|, \qquad (3.2) $$

where $\tilde{\gamma}_{mk} = cmku_\ell/(1-cmku_\ell) \approx cmku_\ell$, for a modest constant c independent of the dimensions of the problem. Note that this bound is valid even with column pivoting, since computing the QRCP decomposition $XP = QR$ of X is equivalent to computing the unpivoted QR decomposition of $XP$. Assumption (2.1) is thus satisfied with $b_1 \approx cmk\sqrt{k}$. The DECOMPRESS kernel is a standard matrix multiplication between the low-rank factors Q and V hence assumption (2.2) is satisfied with $b_2 = k$ by [16, p. 71]. Finally, to bound the error introduced by the RECOMPRESS kernel we must analyze Algorithm 3.1. Since the algorithm simply consists of standard QR, QRCP, and matrix multiplication operations, this analysis is straightforward and relies on standard error bounds from the literature. Although we omit it here for the sake of conciseness, we have performed this analysis and found that (2.3) is satisfied with $b_3 \approx cn(k^{3/2} + \tilde{k}^{3/2}) + k$, where k is the rank of $Q_{in}V_{in}^T$ before recompression and $\tilde{k}$ is the rank of $Q_{out}V_{out}^T$ after recompression.

Finally, we discuss the cost of Algorithm 2.1 when using truncated QRCP decomposition. The truncated QRCP decomposition (3.1) can be computed in $4mnk + o(mnk)$ flops [12], so that assumption (2.29a) is satisfied with $c_1 = 4$. The DECOMPRESS kernel is a matrix multiplication which requires $2mnk$ flops, so that assumption (2.29b) is satisfied with $c_2 = 2$. For the RECOMPRESS kernel, Algorithm 3.1 requires $(6m + 2n)k^2 + o(mk^2 + nk^2)$ flops, so that assumption (2.29c) is satisfied since $(m+n)k^2 = o(mn)$.

**3.2. Randomized SVD decomposition.** Our second choice of LRA algorithm is the randomized SVD decomposition. Several variants have been proposed; in this article we use the one described in Algorithm 3.2, which was proposed by Martinsson and Voronin [22].

The algorithm consists of two phases. The first phase (Lines 3 to 11) iteratively builds a QB decomposition of the matrix such that $\|X - QB\| \le \varepsilon_\ell\|X\|$, where $Q \in \mathbb{R}^{m \times k}$ has orthonormal columns and where the dimension k is hopefully close to $\text{RANK}(X, \varepsilon_\ell)$. To do so, a random Gaussian matrix $\Omega \in \mathbb{R}^{n \times b}$ is drawn for a given block size b and used to sample the matrix $Y = X\Omega$ (Line 5). Then Y is added to the basis Q while preserving the orthonormality: this can for example be accomplished using the Gram-Schmidt algorithm (Line 6). This process is repeated until the columns of Q are a sufficiently good approximation of the column space of X, that is, until $\|X - QQ^T X\|$ is small. In order to efficiently compute this quantity, the matrix $B = Q^T X$ is formed block by block (Line 7). Then, in the second phase (Lines 12 and 13), a truncated SVD of X can easily be obtained by computing the
LRA kernel in our mixed precision iterative refinement approach (Algorithm 2.1): in particular, we check that the algorithm satisfies the assumptions (2.1)–(2.3) of the error analysis and the assumptions (2.29a)–(2.29c) of the complexity analysis. We also explain how to perform the RECOMPRESS kernel based on the specific LRA choice.

**3.1. Truncated QRCP decomposition.** Our first choice of LRA algorithm is a truncated QRCP decomposition. QRCP decomposes the original matrix $X \in \mathbb{R}^{m \times n}$ as

$$XP = QR$$

where $Q \in \mathbb{R}^{m \times n}$ is a matrix with orthonormal columns, $R \in \mathbb{R}^{n \times n}$ is an upper triangular matrix, and $P \in \mathbb{R}^{n \times n}$ is a permutation matrix. Thanks to pivoting, the QRCP decomposition can be used as a rank-revealing algorithm because the norm of the trailing submatrix at each step $k$ of the QR factorization, $\|XP - Q_k R_k\|$, is monotonically decreasing. Therefore, we can stop the QR factorization as soon as this norm becomes smaller than the target tolerance $\varepsilon_\ell$. We obtain

$$X \approx Q_k V_k^T, \quad Q_k \in \mathbb{R}^{m \times k}, \quad V_k = P R_k^T \in \mathbb{R}^{n \times k}, \quad (3.1)$$

where $k = \text{RANK}(X, \varepsilon_\ell)$. To avoid the need to apply the permutation $P$ whenever we need to apply $X$, we form and store $V_k = P R_k^T$.

We can implement a suitable RECOMPRESS algorithm using this truncated QRCP decomposition. Several versions are possible; we describe in Algorithm 3.1 the one that we will use in this article. Given a (non-optimal) truncated QRCP decomposition $Q_{in}V_{in}^T$, Algorithm 3.1 recompresses it into an optimal LRA as follows. First, a thin QR factorization $\overline{QR}$ of $V_{in}$ is computed. Then, the product $Q_{in}R^T$ is formed and a truncated QRCP decomposition $Q_{out}\overline{V}^T$ is computed at the desired target accuracy $\varepsilon$. This yields the recompressed left factor $Q_{out}$; the recompressed right factor $V_{out}$ is obtained by forming $\overline{QV}$.

**Algorithm 3.1 RECOMPRESS algorithm using truncated QRCP decomposition.**

**Input:** a truncated QRCP decomposition $Q_{in}V_{in}^T$ of the form (3.1); $\varepsilon$, the target accuracy.

**Output:** a recompressed QRCP decomposition $Q_{out}V_{out}^T$.

1. Compute the thin QR factorization $\overline{QR} = V_{in}$.

2. Compute the truncated QRCP decomposition $Q_{out}\overline{V}^T = Q_{in}R^T$ with low-rank truncation threshold $\varepsilon$.

3. $V_{out} \leftarrow \overline{QV}$

We now check that the assumptions (2.1)–(2.3) of the error analysis are satisfied for truncated QRCP and determine the value of the constants in the error bounds. For (2.1), we can analyze the truncated QRCP by separating the truncation and rounding errors. As explained above, we can control the size of the truncation error by stopping the QRCP decomposition once the approximation error falls below the desired threshold $\varepsilon_\ell$, so that in exact arithmetic we obtain

$$Q_k R_k = XP + E, \quad \|E\| \le \varepsilon_\ell \|X\|.$$

To account for the rounding errors, we use standard analysis of QR decomposition [18, p. 361], which shows that the computed QR factors satisfy the backward error bound

$$\hat{Q}_k \hat{R}_k = XP + \Delta X + E, \quad \|\Delta X\| \le \sqrt{k} \tilde{\gamma}_{mk} \|X\|, \quad (3.2)$$
**Algorithm 3.2 Randomized SVD decomposition.**

**Input:** $X \in \mathbb{R}^{m \times n}$, a target accuracy $\varepsilon_\ell$, and a block size $b$.

**Output:** a truncated SVD $U\Sigma V^T$ decomposition of $X$.

1. Initialize $Q$ and $B$ to empty matrices.

2. $n_X = \|X\|$

3. **repeat**

4. Draw a random Gaussian matrix $\Omega \in \mathbb{R}^{n \times b}$.

5. $Y = X\Omega$

6. $Q_b = qr(Y - Q(Q^T Y))$

7. $B_b = Q_b^T X$

8. $Q \leftarrow [Q \ Q_b]$

9. $B \leftarrow \begin{bmatrix} B \\ B_b \end{bmatrix}$

10. $X \leftarrow \bar{X} - Q_b B_b$

11. **until** $\|X\| \le \varepsilon_\ell n_X$

12. Compute the truncated SVD decomposition $B \approx \bar{U}\Sigma V^T$ at accuracy $\varepsilon_\ell$.

13. $U = Q\bar{U}$

truncated SVD of the lower dimensional matrix $B \in \mathbb{R}^{k \times n}$: if $B = \bar{U}\Sigma V^T$, then
$X = U\Sigma V^T$ with $U = Q\bar{U}$.

499

500

501

502

503

504

505

506

507

We chose to use this specific randomized SVD variant because it presents several
advantages. It is blocked, which allows for an efficient implementation on modern
hardware. Moreover it provides a cheap yet reliable error estimation. Thanks to
blocking and error estimation, Algorithm 3.2 can adaptively reveal the numerical
rank of the matrix at the requested accuracy $\varepsilon_\ell$; no a priori knowledge on the rank is
thus necessary. Finally, we have experimentally observed Algorithm 3.2 to be more
accurate than other randomized SVD variants that we have tested.

508

509

510

511

512

513

514

515

516

In order to perform the RECOMPRESS operation using randomized SVD, Algorithm 3.2 needs to be slightly adapted. The algorithm takes as input a (non-optimal) low-rank matrix $X = U\Sigma V^T$ and seeks to recompress it, without forming the full $X$. Lines 5 and 7 of Algorithm 3.2 are matrix products and can thus exploit the low-rank structure of $X$. Line 10 is a subtraction between two low-rank matrices, which requires no operations (the low-rank factors can simply be concatenated). The only difficulty lies on Line 11, which requires to compute $\|X\|$ to control the error and stop the algorithm. In order to compute $\|X\|$ without forming $X$, we orthonormalize one of the two low-rank factors and compute the norm of the other one.

517

518

519

520

521

522

523

524

We now check if the assumptions (2.1)–(2.3) of the error analysis are satisfied for
the randomized SVD. To do so, we rely on the error analysis of Connolly, Higham,
and Pranesh [7], which determines error bounds for Algorithm 3.2 in floating-point
arithmetic. By [7, Corollary 2.4], assumption (2.1) is satisfied with $b_1 \approx \sqrt{n}mk$.
Assumption (2.2) is satisfied with $b_2 = k + 1$ since decompressing $U\Sigma V^T$ can be
achieved with a matrix–matrix product and a scaling. Finally, the analysis of [7] does
not directly apply to the RECOMPRESS version of the algorithm discussed above, but
we expect its numerical behavior to be similar to the original algorithm.

525

526

527

528

Finally, we discuss the cost of Algorithm 3.2. If $X$ is a full $m \times n$ matrix, Algorithm 3.2 costs $6mnk + o(mnk)$ flops [22, Eqn. (25)], so that assumption (2.29a) is satisfied with $c_1 = 6$. It is easy to extend [22, Eqn. (25)] to the case where $X = U\Sigma V^T$ is represented under low-rank form to check that the RECOMPRESS vari-
where $\tilde{\gamma}_{mk} = cmku_{\ell}/(1-cmku_{\ell}) \approx cmku_{\ell}$, for a modest constant $c$ independent of the dimensions of the problem. Note that this bound is valid even with column pivoting, since computing the QRCP decomposition $XP = QR$ of $X$ is equivalent to computing the unpivoted QR decomposition of $XP$. Assumption (2.1) is thus satisfied with $b_1 \approx cmk\sqrt{k}$. The DECOMPRESS kernel is a standard matrix multiplication between the low-rank factors $Q$ and $V$ hence assumption (2.2) is satisfied with $b_2 = k$ by [18, p. 71]. Finally, to bound the error introduced by the RECOMPRESS kernel we must analyze Algorithm 3.1. Since the algorithm simply consists of standard QR, QRCP, and matrix multiplication operations, this analysis is straightforward and relies on standard error bounds from the literature. Although we omit it here for the sake of conciseness, we have performed this analysis and found that (2.3) is satisfied with $b_3 \approx cn(k^{3/2} + \tilde{k}^{3/2}) + k$, where $k$ is the rank of $Q_{in}V_{in}^T$ before recompression and $\tilde{k}$ is the rank of $Q_{out}V_{out}^T$ after recompression.

Finally, we discuss the cost of Algorithm 2.1 when using truncated QRCP decomposition. The truncated QRCP decomposition (3.1) can be computed in $4mnk + o(mnk)$ flops [14], so that assumption (2.29a) is satisfied with $c_1 = 4$. The DECOMPRESS kernel is a matrix multiplication which requires $2mnk$ flops, so that assumption (2.29b) is satisfied with $c_2 = 2$. For the RECOMPRESS kernel, Algorithm 3.1 requires $(6m+2n)k^2 + o(mk^2+nk^2)$ flops, so that assumption (2.29c) is satisfied since $(m+n)k^2 = o(mn)$.

**3.2. Randomized SVD decomposition.** Our second choice of LRA algorithm is the randomized SVD decomposition. Several variants have been proposed; in this article we use the one described in Algorithm 3.2, which was proposed by Martinsson and Voronin [24].

<caption>**Algorithm 3.2 Randomized SVD decomposition.**</caption>

**Input:** $X \in \mathbb{R}^{m \times n}$, a target accuracy $\varepsilon_\ell$, and a block size $b$.

**Output:** a truncated SVD $U\Sigma V^T$ decomposition of $X$.

1: Initialize $Q$ and $B$ to empty matrices.

2: $n_X = \|X\|$

3: **repeat**

4: Draw a random Gaussian matrix $\Omega \in \mathbb{R}^{n \times b}$.

5: $Y = X\Omega$

6: $Q_b = \mathbf{qr}(Y - Q(Q^T Y))$

7: $B_b = Q_b^T X$

8: $Q \leftarrow [Q Q_b]$

9: $B \leftarrow \begin{bmatrix} B \\ B_b \end{bmatrix}$

10: $X \leftarrow X - Q_b B_b$

11: **until** $\|X\| \le \varepsilon_\ell n_X$

12: Compute the truncated SVD decomposition $B \approx \bar{U}\Sigma V^T$ at accuracy $\varepsilon_\ell$.

13: $U = Q\bar{U}$

The algorithm consists of two phases. The first phase (Lines 3 to 11) iteratively builds a QB decomposition of the matrix such that $\|X - QB\| \le \varepsilon_\ell \|X\|$, where $Q \in \mathbb{R}^{m \times k}$ has orthonormal columns and where the dimension $k$ is hopefully close to $\text{RANK}(X, \varepsilon_\ell)$. To do so, a random Gaussian matrix $\Omega \in \mathbb{R}^{n \times b}$ is drawn for a given block size $b$ and used to sample the matrix $Y = X\Omega$ (Line 5). Then $Y$ is added to the basis $Q$ while preserving the orthonormality: this can for example be
ant of Algorithm 3.2 described above has a cost in $O((m+n)k^2) = o(mn)$ flops, so that assumption (2.29c) is indeed satisfied. Finally, assumption (2.29b) is satisfied with $c_2 = 2$ since the DECOMPRESS kernel is simply a matrix–matrix product.

### 4. Application to tensor low-rank approximation.
In this section, we explore the application of our method to three different tensor decompositions, namely Tucker [29], TT [27], and HT [13]. All three decompositions provide direct methods for the LRA and RECOMPRESS kernels that can guarantee a prescribed accuracy $\epsilon$, as well as fast DECOMPRESS routines, rendering them amenable to Algorithm 2.1. We first give a brief description of these decompositions as well as their corresponding LRA, RECOMPRESS, and DECOMPRESS kernels, then provide an error and complexity analysis with respect to the assumptions (2.1)–(2.3) and (2.29a)–(2.29c).

#### 4.1. Tucker, TT, and HT decompositions.
The Tucker decomposition represents a tensor $X \in \mathbb{R}^{n_1 \times \cdots \times n_d}$ with a core tensor $G \in \mathbb{R}^{r_1 \times \cdots \times r_d}$ and with factor matrices $U = \{U^{(i)}\}$ with orthonormal columns, where $U^{(i)} \in \mathbb{R}^{n_i \times r_i}$ for $i \in \{1, \dots, d\}$. Each element of a Tucker tensor is given by

$$X_{i_1, \dots, i_d} = \sum_{\alpha_1, \dots, \alpha_d}^{r_1, \dots, r_d} G_{\alpha_1, \dots, \alpha_d} U_{i_1, \alpha_1}^{(1)} \cdots U_{i_d, \alpha_d}^{(d)} \quad (4.1)$$

The TT format consists of a sequence of 3rd-order tensors $G^{(i)}$ linked by the first and third dimensions that correspond to ranks. Each element of a TT tensor is given by

$$X_{i_1, \dots, i_d} = \sum_{\alpha_0, \dots, \alpha_d}^{r_0, \dots, r_d} G_{\alpha_0, i_1, \alpha_1}^{(1)} \cdots G_{\alpha_{d-1}, i_d, \alpha_d}^{(d)} \quad (4.2)$$

The HT decomposition represents $X$ as a binary tree of tensors whose leaves are matrices $U^{(i)} \in \mathbb{R}^{n_i \times r_i}$ and the rest are 3rd order tensors. Each element of an HT tensor is given by

$$X_{i_1, \dots, i_d} = \sum_{\alpha_0, \dots, \alpha_{2d-2}}^{r_0, \dots, r_{2d-2}} G_{\alpha_0, \alpha_1, \alpha_2}^{(1)} \cdots G_{\alpha_{d-2}, \alpha_{2d-3}, \alpha_{2d-2}}^{(d-1)} U_{i_1, \alpha_{d-1}}^{(1)} \cdots U_{i_d, \alpha_{2d-2}}^{(d)} \quad (4.3)$$

We illustrate these decompositions in Figure 4.1 using a tensor network diagram [5] where each circle represents a tensor and each edge corresponds to a dimension. For each decomposition, the corresponding LRA method yields a network of core tensors (or just “cores”), whose contraction/multiplication along the inner dimensions (of size $r_i$, $i = 1: \bar{d}$) yields the full tensor. The outer dimensions (of size $n_i$, $i = 1: d$) correspond to dimensions of the full tensor $X$.

The HOSVD [29], TTSVD [27], and HTSVD [13] algorithms can be used as LRA kernels to compute the Tucker, TT, and HT decompositions of a given full tensor $X$, respectively. All three are direct methods that can achieve any prescribed accuracy $\epsilon$ in exact arithmetic. To the best of our knowledge, the effect of rounding errors on these methods in finite precision arithmetic has not been analyzed in the literature. While these methods are therefore not known to be stable, empirical experiments suggest that they still reliably achieve an accuracy of order $\epsilon$ when run in a finite precision with a unit roundoff $u$ sufficiently smaller than $\epsilon$. Formally proving this result is an open problem that deserves a dedicated study.

We provide a high level pseudo-code that encapsulates these methods in Algorithm 4.1. In essence, all three algorithms perform $\bar{d}$ (defined in (1.1)) successive SVDs
accomplished using the Gram–Schmidt algorithm (Line 6). This process is repeated
until the columns of $Q$ are a sufficiently good approximation of the column space of
$X$, that is, until $\|X - QQ^T X\|$ is small. In order to efficiently compute this quantity,
the matrix $B = Q^T X$ is formed block by block (Line 7). Then, in the second phase
(Lines 12 and 13), a truncated SVD of $X$ can easily be obtained by computing the
truncated SVD of the lower dimensional matrix $B \in \mathbb{R}^{k \times n}$: if $B = \bar{U}\Sigma V^T$, then
$X = U\Sigma V^T$ with $U = \bar{Q}\bar{U}$.

We chose to use this specific randomized SVD variant because it presents several
advantages. It is blocked, which allows for an efficient implementation on modern
hardware. Moreover it provides a cheap yet reliable error estimation. Thanks to
blocking and error estimation, Algorithm 3.2 can adaptively reveal the numerical
rank of the matrix at the requested accuracy $\epsilon_\ell$; no a priori knowledge on the rank is
thus necessary. Finally, we have experimentally observed Algorithm 3.2 to be more
accurate than other randomized SVD variants that we have tested.

In order to perform the RECOMPRESS operation, one could simply use Algorithm 3.1; however, if one wishes to use randomized SVD not only for the LRA but also the RECOMPRESS operation, then Algorithm 3.2 can still be used with some modifications. These modifications are described in Appendix A and result in Algorithm A.1.

We now check if the assumptions (2.1)–(2.3) of the error analysis are satisfied for
the randomized SVD. To do so, we rely on the error analysis of Connolly, Higham,
and Pranesh [9], which determines error bounds for Algorithm 3.2 in floating-point
arithmetic. Note that their analysis assumes that the deterministic SVD on line 12
is computed with a numerically stable algorithm. By [9, Corollary 2.4], assumption
(2.1) is satisfied with $b_1 \approx \sqrt{n}mk$. Assumption (2.2) is satisfied with $b_2 = k+1$ since
decompressing $U\Sigma V^T$ can be achieved with a matrix–matrix product and a scaling.
Finally, the analysis of [9] does not directly apply to the RECOMPRESS version of the
algorithm discussed above, but we expect its numerical behavior to be similar to the
original algorithm.

Next, we discuss the cost of Algorithm 3.2. If $X$ is a full $m \times n$ matrix, Algorithm 3.2 costs $6mnk + o(mnk)$ flops [24, Eqn. (25)], so that assumption (2.29a) is satisfied with $c_1 = 6$. It is easy to extend [24, Eqn. (25)] to the case where $X = U\Sigma V^T$ is represented under low-rank form to check that the RECOMPRESS variant of Algorithm 3.2 described above has a cost in $O((m+n)k^2) = o(mn)$ flops, so that assumption (2.29c) is indeed satisfied. Finally, assumption (2.29b) is satisfied with $c_2 = 2$ since the DECOMPRESS kernel is simply a matrix–matrix product.

**4. Application to tensor low-rank approximation.** In this section, we explore the application of our method to three different tensor decompositions, namely Tucker [32], TT [29], and HT [15]. All three decompositions provide direct methods for the LRA and RECOMPRESS kernels that can guarantee a prescribed accuracy $\epsilon$, as well as fast DECOMPRESS routines, rendering them amenable to Algorithm 2.1. We first give a brief description of these decompositions as well as their corresponding LRA, RECOMPRESS, and DECOMPRESS kernels, then provide an error and complexity analysis with respect to the assumptions (2.1)–(2.3) and (2.29a)–(2.29c).

**4.1. Tucker, TT, and HT decompositions.** The Tucker decomposition represents a tensor $X \in \mathbb{R}^{n_1 \times \cdots \times n_d}$ with a core tensor $G \in \mathbb{R}^{r_1 \times \cdots \times r_d}$ and with factor matrices $U = \{U^{(i)}\}$ with orthonormal columns, where $U^{(i)} \in \mathbb{R}^{n_i \times r_i}$ for $i \in \{1, \dots, d\}$.
<figure><img src="image_4.png" /><figcaption>Fig. 4.1: A *d*-order tensor (left) and its representation in Tucker, TT, and HT decompositions (right) using tensor network diagram.</figcaption></figure>

on matricizations of $X$ to isolate each core of the decomposition and determine its corresponding rank, as indicated on Line 4. Then, the left orthogonal component $U$ of the SVD is reshaped and added to the decomposition $F_X$, and the non-orthogonal component is reshaped into a tensor to be further factored in the subsequent iterations. At the end, the decomposed tensor $F_X$ consisting of $\bar{d} + 1$ tensors is returned.

As the LRA kernel isolates the cores through a series of $\bar{d}$ SVDs on $X$, the DECOMPRESS kernel performs successive tensor contractions on neighboring cores of $F_X$ through the inner dimensions that link them to obtain the full tensor $X$. In Algorithm 4.2, we provide a generic DECOMPRESS kernel that performs these contractions
<figure><img src="image_3.png" /><figcaption>Fig. 4.1: A d-order tensor (left) and its representation in Tucker, TT, and HT decompositions (right) using tensor network diagram.</figcaption></figure>

Each element of a Tucker tensor is given by

$$X_{i_1, \dots, i_d} = \sum_{\alpha_1, \dots, \alpha_d}^{r_1, \dots, r_d} G_{\alpha_1, \dots, \alpha_d} U_{i_1, \alpha_1}^{(1)} \dots U_{i_d, \alpha_d}^{(d)} \quad (4.1)$$

The TT format consists of a sequence of 3rd-order tensors $G^{(i)}$ linked by the first and third dimensions that correspond to ranks. Each element of a TT tensor is given
**Algorithm 4.1 LRA algorithm for a tensor.**

**Input:** a tensor $X$ and a target accuracy $\epsilon$.

**Output:** $F_X$: a low-rank tensor decomposition of $X$.

1: $\bar{\epsilon} \leftarrow \epsilon/\bar{d}$

2: $F_X \leftarrow []$

3: **for** $i = 1$ **to** $\bar{d}$ **do**

4: &nbsp;&nbsp;&nbsp;&nbsp;[$U, Y$] $\leftarrow$ SVD(MATRICIZE($X$), $\bar{\epsilon}$)

5: &nbsp;&nbsp;&nbsp;&nbsp;$F_X \leftarrow [F_X, \RESHAPE(U)]$

6: &nbsp;&nbsp;&nbsp;&nbsp;$X \leftarrow \RESHAPE(Y)$

7: **end for**

8: &nbsp;&nbsp;&nbsp;&nbsp;$F_X \leftarrow [F_X, \RESHAPE(X)]$

**Algorithm 4.2 DECOMPRESS algorithm for a tensor decomposition.**

**Input:** a low-rank tensor decomposition $F_X$.

**Output:** the full tensor $X$ corresponding to $F_X$.

1: $X \leftarrow F_X[1]$

2: **for** $i = \bar{d}$ **to** 1 in reverse order **do**

3: &nbsp;&nbsp;&nbsp;&nbsp;$X \leftarrow X \times_{\{r_i\}} F_X[i+1]$

4: **end for**

579 in the reverse order of SVDs in Algorithm 4.1 for simplicity. This is done on Line 3 where $\times_{\{r_i\}}$ represents the contraction/multiplication across the inner dimension of size $r_i$. In practice, these contractions could be performed in an arbitrary order among the $\bar{d}$ inner dimensions as tensor contraction is associative.

583 Finally, the RECOMPRESS kernel of the TT decomposition consists of an orthogonalization and a rank pruning phase [27]. The goal of the orthogonalization phase is to “sweep” all non-orthogonal components in $F_X$ into a single (the first or the last) core. To do this, a QR factorization is performed on a matricization of the core at the opposite end. The orthogonal component $Q$ is reshaped to replace the treated core, and the non-orthogonal component $R$ is passed over to the next core through a tensor-matrix contraction. This process is repeated $\bar{d}$ times until all cores except one become orthogonal, which finalizes the orthogonalization phase. Next, the rank pruning phase performs an SVD on a matricization of the non-orthogonal core to reduce the size of an inner dimension. The orthogonal part $U$ obtained from the SVD is kept as the pruned core whereas the non-orthogonal part $\Sigma V^T$ is integrated into the next core to be truncated through tensor-matrix contraction, thereby still keeping all cores except one orthogonal. The process is similarly repeated on all $\bar{d}$ inner dimensions in the reverse order of orthogonalization, which finalizes the rank pruning phase. The standard RECOMPRESS kernel in the HT format [13] similarly sweeps all non-orthogonal components into a single core (the root of the tree) with a sequence of QR factorizations followed by matrix multiplies through a bottom-up traversal of the tree. However, in the rank pruning phase, it processes inner cores level-by-level starting from the root, and performs the rank truncation on the matricization of the tensor sub-network involving all tensors in the path from the root down to the inner core that is being truncated. Though this variant remains efficient (thanks to the re-use of partial tensor contractions from the previous level through memorization) and provides more coarse-grain parallelism (with respect to TTSVD for instance, as each node in a tree level can be truncated independently), it requires to form a Gram
by

$$X_{i_1, \dots, i_d} = \sum_{\alpha_0, \dots, \alpha_d}^{r_0, \dots, r_d} G_{\alpha_0, i_1, \alpha_1}^{(1)} \dots G_{\alpha_{d-1}, i_d, \alpha_d}^{(d)} \quad (4.2)$$

The HT decomposition represents $X$ as a binary tree of tensors whose leaves are matrices $U^{(i)} \in \mathbb{R}^{n_i \times r_i}$ and the rest are 3rd order tensors. Each element of an HT tensor is given by

$$X_{i_1, \dots, i_d} = \sum_{\alpha_0, \dots, \alpha_{2d-2}}^{r_0, \dots, r_{2d-2}} G_{\alpha_0, \alpha_1, \alpha_2}^{(1)} \dots G_{\alpha_{d-2}, \alpha_{2d-3}, \alpha_{2d-2}}^{(d-1)} U_{i_1, \alpha_{d-1}}^{(1)} \dots U_{i_d, \alpha_{2d-2}}^{(d)} \quad (4.3)$$

We illustrate these decompositions in Figure 4.1 using a tensor network diagram [7] where each circle represents a tensor and each edge corresponds to a dimension. For each decomposition, the corresponding LRA method yields a network of core tensors (or just “cores”), whose contraction/multiplication along the inner dimensions (of size $r_i$, $i = 1: \bar{d}$) yields the full tensor. The outer dimensions (of size $n_i$, $i = 1: d$) correspond to dimensions of the full tensor $X$.

**Algorithm 4.1 LRA algorithm for a tensor.**

**Input:** a tensor $X$ and a target accuracy $\varepsilon$.

**Output:** $F_X$: a low-rank tensor decomposition of $X$.

1: $\bar{\varepsilon} \leftarrow \varepsilon/\bar{d}$

2: $F_X \leftarrow []$

3: **for** $i = 1$ **to** $\bar{d}$ **do**

4: &nbsp;&nbsp;&nbsp;&nbsp;[$U, Y$] $\leftarrow$ SVD(MATRICIZE($X$), $\bar{\varepsilon}$)

5: &nbsp;&nbsp;&nbsp;&nbsp;$F_X \leftarrow [F_X, \text{RESHAPE}(U)]$

6: &nbsp;&nbsp;&nbsp;&nbsp;$X \leftarrow \text{RESHAPE}(Y)$

7: **end for**

8: $F_X \leftarrow [F_X, \text{RESHAPE}(X)]$

The HOSVD [32], TTSVD [29], and HTSVD [15] algorithms can be used as LRA kernels to compute the Tucker, TT, and HT decompositions of a given full tensor $X$, respectively. All three are direct methods that can achieve any prescribed accuracy $\varepsilon$ in exact arithmetic. To the best of our knowledge, the effect of rounding errors on these methods in finite precision arithmetic has not been analyzed in the literature. While these methods are therefore not known to be stable, empirical experiments suggest that they still reliably achieve an accuracy of order $\varepsilon$ when run in a finite precision with a unit roundoff $u$ sufficiently smaller than $\varepsilon$. Formally proving this result is an open problem that deserves a dedicated study.

We provide a high level pseudo-code that encapsulates these methods in Algorithm 4.1. In essence, all three algorithms perform $\bar{d}$ (defined in (1.1)) successive SVDs on matricizations of $X$ to isolate each core of the decomposition and determine its corresponding rank, as indicated on Line 4. Then, the left orthogonal component $U$ of the SVD is reshaped and added to the decomposition $F_X$, and the non-orthogonal component is reshaped into a tensor to be further factored in the subsequent iterations. At the end, the decomposed tensor $F_X$ consisting of $\bar{d} + 1$ tensors is returned.

As the LRA kernel isolates the cores through a series of $\bar{d}$ SVDs on $X$, the DECOMPRESS kernel performs successive tensor contractions on neighboring cores of $F_X$ through the inner dimensions that link them to obtain the full tensor $X$. In Algorithm 4.2, we provide a generic DECOMPRESS kernel that performs these contractions
matrix to remain efficient when performing rank truncation (otherwise the size of the matricized tensor explodes, yet its Gramian stays compact). We observed that this method greatly suffers from numerical instability due to this inherent limitation, particularly when soliciting a high precision in RECOMPRESS, as it squares the condition number of the matricized tensor. We instead devised an alternative RECOMPRESS kernel, reminiscent of that of the TT decomposition, for HT and Tucker decompositions. This variant similarly involves an orthogonalization and rank pruning through $\tilde{d}$ QRs and SVDs followed by matrix multiplies, albeit with no Gram matrix forming, and proved to be much more robust numerically yet still efficient. We defer the detailed presentation of these algorithms and their numerical and cost analysis to a future work. Due to these intricacies, we do not provide a high level RECOMPRESS kernel encapsulating all three decompositions; however, for the purposes of this paper, we only use the fact that all three decompositions provide a RECOMPRESS kernel performing $\tilde{d}$ QR and SVD computations followed by matrix multiplies on its cores, which suffices for our analysis in the next section.

**4.2. Error and complexity analysis.** For the error analysis, as previously mentioned the stability of tensor computations in finite precision arithmetic is an open problem that deserves a dedicated study. However, we may mention that the RESHAPE and MATRICIZE operations do not introduce any error as they merely involve data shuffling, nor do they have any impact on the Frobenius norm of the successive matricizations. Moreover, the LRA kernel is based on a chain of SVDs, the DECOMPRESS kernel is based on a chain of matrix multiplies, and the RECOMPRESS kernel on a chain of QRs, SVDs, and matrix multiplies. All of these basic linear algebra building blocks have stable implementations and so we can reasonably expect the numerical behavior of the overall tensor computation to satisfy assumptions (2.1)–(2.3).

For the complexity analysis of the LRA step, we assume that the first SVD reduces the size of the tensor significantly, that is, to $o(p)$, which renders the cost of subsequent SVDs negligible. In this case, we can use the same constant $c_1 = 6$ as in the case of randomized SVD in subsection 3.2. With the same assumption, the cost of DECOMPRESS in Algorithm 4.2 will also be dominated by the last contraction, which is essentially a matrix multiplication on permuted tensors across the first inner dimension. Thus, we can similarly use the constant $c_2 = 2$ as in subsection 3.2, since the cost of the previous contractions is in $o(p)$ in this case as well. Even without this assumption, we can find constants $c_1 = 6t$ and $c_2 = 2t$ with the same $t > 1$, since intermediate tensor sizes in each SVD and contraction steps in Algorithm 4.1 and Algorithm 4.2 stay the same. This keeps the cost ratio of these two steps constant across all $\tilde{d}$ dimensions yielding the same $t$ for $c_1$ and $c_2$; we skip further details for brevity. Finally, the cost of RECOMPRESS remains in $o(p)$ as it involves QRs and SVDs on the matricizations of 2nd or 3rd order cores in the decomposition, whose size is in $o(p)$ by the low-rank assumption. Therefore, we satisfy the assumptions (2.29a)–(2.29c) of the complexity analysis in the tensor case.

# 5. Numerical experiments.

**5.1. Experimental setting.** We now test our iterative refinement approach experimentally. We developed a MATLAB code that implements Algorithm 2.1 and can use as LRA kernel any of the matrix and tensor LRA algorithms discussed in the previous two sections: QRCP, randomized SVD, HOSVD, TTSVD, and HTSVD.

For the matrix algorithms (QRCP and randomized SVD), we use our own imple-
**Algorithm 4.2 DECOMPRESS algorithm for a tensor decomposition.**

**Input:** a low-rank tensor decomposition $F_X$.

**Output:** the full tensor $X$ corresponding to $F_X$.

1: $X \leftarrow F_X[1]$

2: **for** $i = \bar{d}$ **to** 1 in reverse order **do**

3: &nbsp;&nbsp;&nbsp;&nbsp;$X \leftarrow X \times_{\{r_i\}} F_X[i+1]$

4: **end for**

in the reverse order of SVDs in Algorithm 4.1 for simplicity. This is done on Line 3 where $\times_{\{r_i\}}$ represents the contraction/multiplication across the inner dimension of size $r_i$. In practice, these contractions could be performed in an arbitrary order among the $\bar{d}$ inner dimensions as tensor contraction is associative.

Finally, for the TT format we use the standard RECOMPRESS kernel (called “rounding” in the tensor literature) which relies on a sequence of QRs and SVDs on the matricizations of the cores (we refer to [29] for a detailed description). For HT and Tucker decompositions, we employ an adaptation of this TT RECOMPRESS kernel to these formats, instead of the standard RECOMPRESS kernel proposed in [15] mainly to avoid the Gram matrix formation therein.

**4.2. Discussion of the error and complexity.** For the error analysis, as previously mentioned, tensor computations are experimentally observed to behave stably in finite precision arithmetic, but a formal proof of this fact is not known. The LRA kernel is based on a chain of SVDs, the DECOMPRESS kernel is based on a chain of matrix multiplies, and the RECOMPRESS kernel on a chain of QRs, SVDs, and matrix multiplies. Despite the fact that all of these basic linear algebra building blocks have stable implementations, it is not sufficient to directly conclude on the stability of the overall tensor computation, because the composition of stable algorithms is not necessarily stable [5]. Therefore, the stability of tensor algorithms is an open problem that deserves a dedicated study.

For the complexity analysis of the LRA step, we assume that the first SVD reduces the size of the tensor significantly, that is, to $o(p)$, which renders the cost of subsequent SVDs negligible. In this case, we can use the same constant $c_1 = 6$ as in the case of randomized SVD in subsection 3.2. With the same assumption, the cost of DECOMPRESS in Algorithm 4.2 will also be dominated by the last contraction, which is essentially a matrix multiplication on permuted tensors across the first inner dimension. Thus, we can similarly use the constant $c_2 = 2$ as in subsection 3.2, since the cost of the previous contractions is in $o(p)$ in this case as well. Even without this assumption, we can find constants $c_1 = 6t$ and $c_2 = 2t$ with the same $t > 1$, since intermediate tensor sizes in each SVD and contraction steps in Algorithm 4.1 and Algorithm 4.2 stay the same. This keeps the cost ratio of these two steps constant across all $\bar{d}$ dimensions yielding the same $t$ for $c_1$ and $c_2$; we skip further details for brevity. Finally, the cost of RECOMPRESS remains in $o(p)$ as it involves QRs and SVDs on the matricizations of 2nd or 3rd order cores in the decomposition, whose size is in $o(p)$ by the low-rank assumption. Therefore, we satisfy the assumptions (2.29a)-(2.29c) of the complexity analysis in the tensor case.

# 5. Numerical experiments.

**5.1. Experimental setting.** We now test our iterative refinement approach experimentally. We developed a MATLAB code that implements Algorithm 2.1 and
mentation. For the tensor algorithms (HOSVD, TTSVD, and HTSVD), we rely on the implementations from the libraries described in [3], [26], and [21], respectively, with some adjustments. We use MATLAB version R2019a throughout the experiments.

In the experiments, the high precision $u$ is set to double precision (fp64 arithmetic, with unit roundoff $u = 2^{-53} \approx 1 \times 10^{-16}$) and the use of various low precisions $u_\ell$ is simulated with the chop library of Higham and Pranesh [18]: single precision (fp32 arithmetic, with unit roundoff $u_\ell = 2^{-24} \approx 6 \times 10^{-8}$) and half precision (fp16 and bfloat16 arithmetics, with unit roundoff $u_\ell = 2^{-11} \approx 5 \times 10^{-4}$ and $u_\ell = 2^{-8} \approx 4 \times 10^{-3}$, respectively).

For the low-rank truncation thresholds, we set $\varepsilon = 10^{-13}$ and $\varepsilon_\ell = \theta_\ell u_\ell$, where $\theta_\ell \le 1$ is a scaling factor necessary to control the effect of rounding errors on the ability of the LRA to detect the correct numerical rank. We analyze in detail the role of this parameter $\theta_\ell$ in subsection 5.3; based on the conclusions of this analysis, we have set $\theta_\ell = 2^{-1}$ for all LRA kernels except randomized SVD and QRCP, for which we have used $\theta_\ell = 2^{-2}$ and $\theta_\ell = 2^{-3}$, respectively.

To test the algorithms, we use randomly generated matrices and tensors with various singular value distributions. More precisely, we compare three types of distributions for the singular values $\sigma_i$:

$$ \sigma_i = \max(\hat{\sigma}_i, 10^{-16}), \quad \hat{\sigma}_i = \begin{cases} 1/i & (\text{linear}) \\ i^{-10} & (\text{power}) \\ e^{-i} & (\text{exponential}) \end{cases} \quad (5.1) $$

These three distributions are illustrated in Figure 5.1. We generate the matrices as $Q_1\Sigma Q_2$ where $Q_1, Q_2$ are random orthogonal matrices and $\Sigma$ is a diagonal matrix with the specified singular values as coefficients. For the tensor experiments involving Tucker and TT decompositions, we generate the low-rank tensor in the Tucker format with $Q^{(1)}, \dots, Q^{(d)}$ random orthonormal matrices and a $d$th-order core tensor $G$ whose coefficient $G_{i_1, \dots, i_d}$ is given by $\sigma_{\max(i_1, \dots, i_d)}$ in (5.1). For the experiments involving the HT format, we generate the low-rank tensor in the HT format whose leaf nodes are matrices with orthogonal columns, and each element $G_{i,j,k}^{(t)}$ of its internal cores is set to $\sigma_{\max(i,j,k)}$.

**5.2. Experimental results.** In the first experiment, we analyze the behavior of Algorithm 2.1 for the matrix case, which is provided in Figure 5.2. We consider three types of matrices based on their singular value distribution (linear, power, and exponential) and two LRA kernels (QRCP and randomized SVD). In each case, we plot the relative error

$$ \eta_i = \frac{\|X - F_i\|}{\|X\|} \quad (5.2) $$

where $F_i$ is the computed low-rank factor of $X$ after $i$ refinement steps ($F_0$ is thus a standard LRA of $X$ in precision $u_\ell$). The number next to each marker indicates the rank of $F_i$ after recompression.

Figure 5.2 shows that for all three matrices and for both QRCP and randomized SVD, Algorithm 2.1 behaves as expected: the error $\eta_i$ is roughly equal to $\varepsilon_\ell^{i+1} = (\theta u_\ell)^{i+1}$ after $i$ refinement steps. Thus, using single precision as the low precision $u_\ell$, we can achieve an accuracy close to double precision with only one refinement step. Naturally, since for QRCP $\theta = 2^{-3}$ is relatively small, $(u_\ell/\theta)^2$ is noticeably larger than the accuracy achieved by a standard double precision LRA; the refined factors are thus not completely as accurate as if computed directly in double precision. This
can use as LRA kernel any of the matrix and tensor LRA algorithms discussed in the previous two sections: QRCP, randomized SVD, HOSVD, TTSVD, and HTSVD.

For the matrix algorithms (QRCP and randomized SVD), we use our own implementation. For the tensor algorithms (HOSVD, TTSVD, and HTSVD), we rely on the implementations from the libraries described in [4], [28], and [23], respectively, with some adjustments. We use MATLAB version R2019a throughout the experiments.

In the experiments, the high precision $u$ is set to double precision (fp64 arithmetic, with unit roundoff $u = 2^{-53} \approx 1 \times 10^{-16}$) and the use of various low precisions $u_\ell$ is simulated with the chop library of Higham and Pranesh [20]: single precision (fp32 arithmetic, with unit roundoff $u_\ell = 2^{-24} \approx 6 \times 10^{-8}$) and half precision (fp16 and bfloat16 arithmetics, with unit roundoff $u_\ell = 2^{-11} \approx 5 \times 10^{-4}$ and $u_\ell = 2^{-8} \approx 4 \times 10^{-3}$, respectively).

For the low-rank truncation thresholds, we set $\varepsilon = 10^{-13}$ and $\varepsilon_\ell = \theta_\ell u_\ell$, where $\theta_\ell \le 1$ is a scaling factor necessary to control the effect of rounding errors on the ability of the LRA to detect the correct numerical rank. We analyze in detail the role of this parameter $\theta_\ell$ in subsection 5.3; based on the conclusions of this analysis, we have set $\theta_\ell = 2^{-1}$ for all LRA kernels except randomized SVD and QRCP, for which we have used $\theta_\ell = 2^{-2}$ and $\theta_\ell = 2^{-3}$, respectively.

To test the algorithms, we use randomly generated matrices and tensors with various singular value distributions. More precisely, we compare three types of distributions for the singular values $\sigma_i$:

$$ \sigma_i = \max(\hat{\sigma}_i, 10^{-16}), \quad \hat{\sigma}_i = \begin{cases} 1/i & (\text{linear}) \\ i^{-10} & (\text{power}) \\ e^{-i} & (\text{exponential}) \end{cases}. \qquad (5.1) $$

These three distributions are illustrated in Figure 5.1. All random matrices are square and of size $100 \times 100$. All random tensors are of fourth order ($d = 4$) and of size $100 \times 100 \times 100 \times 100$. We generate the matrices as $Q_1 \Sigma Q_2$ where $Q_1, Q_2$ are random orthogonal matrices and $\Sigma$ is a diagonal matrix with the specified singular values as coefficients. For the tensor experiments involving Tucker and TT decompositions, we generate the low-rank tensor in the Tucker format with $Q^{(1)}, \dots, Q^{(d)}$ random orthonormal matrices and a $d$th-order core tensor $G$ whose coefficient $G_{i_1 \dots i_d}$ is given by $\sigma_{\max(i_1, \dots, i_d)}$ in (5.1). For the experiments involving the HT format, we generate the low-rank tensor in the HT format whose leaf nodes are matrices with orthogonal columns, and each element $G_{i,j,k}^{(t)}$ of its internal cores is set to $\sigma_{\max(i,j,k)}$.

In addition to these random data sets, we also report results for data sets coming from two real-life applications in section 5.4.

## 5.2. Experimental results.
In the first experiment, we analyze the behavior of Algorithm 2.1 for the matrix case, which is provided in Figure 5.2. We consider three types of matrices based on their singular value distribution (linear, power, and exponential) and two LRA kernels (QRCP and randomized SVD). In each case, we plot the relative error

$$ \eta_i = \frac{\|X - F_i\|}{\|X\|} \qquad (5.2) $$

where $F_i$ is the computed low-rank factor of $X$ after $i$ refinement steps ($F_0$ is thus a standard LRA of $X$ in precision $u_\ell$). The number next to each marker indicates the rank of $F_i$ after recompression.
<figure><img src="image_5.png" /><figcaption>Fig. 5.1: Three types of singular value distributions used in the experiments.</figcaption></figure>

gap can be filled by performing a second refinement step, although this would likely be unnecessary in most practical scenarios.

Similar results are obtained using half precision as the low precision format (fp16 or bfloat16 arithmetics). An error close to the single precision accuracy can be achieved in just one or a few steps. Moreover, we can even reach double precision accuracy if needed, which illustrates an attractive property of Algorithm 2.1: it can reach an accuracy essentially as high as desired while performing most of its operations in a precision as low as desired.

Finally, we discuss the rank behavior of $F_i$ across IR steps. We see that it is roughly equal to RANK($X$, $\eta_i$), the numerical rank of $X$ at accuracy $\eta_i$. Thus, for matrices with rapidly decaying singular values such as the exponential case in Figure 5.1, these ranks tend to be much smaller at the early steps in low precision. Reflecting on the cost of the algorithm, we can expect the use of low precision plus refinement to be particularly cost efficient for such matrices.

Figure 5.3 shows similar plots for the tensor decompositions (TTSVD, HOSVD, and HTSVD). The results for tensors follow the same trend, and lead to the same conclusion as for matrices. The main difference is that the rank of $F_i$ (text labels) is now a vector instead.

Overall, our experiments confirm that Algorithm 2.1 is able to rapidly converge to a high accuracy, while using low precision for the LRA kernel. This observation is valid for all types of matrices and tensors in our test set, and all five LRA kernels that we tested, which shows that Algorithm 2.1 is robust and can work in a wide variety of settings.

**5.3. Role of $\theta_\ell$.** There is a tradeoff in choosing the scaling factor $\theta_\ell$: the larger it is, the faster the convergence (since the error is reduced by a factor roughly equal to $\theta_\ell u_\ell$ at each refinement step), but if it is too large, the rank will no longer be correctly detected and this will lead to a significant rank growth.

We illustrate this in Table 5.1 and Table 5.2, where we compare the convergence of Algorithm 2.1 for different values of $\theta_\ell$. Table 5.1 is for the TTSVD kernel and Table 5.2 is for the QRCP kernel; fp16 arithmetic is used as the low precision $u_\ell$ in both cases. The tables show that the method converges faster as $\theta_\ell$ increases, as expected:
<figure><img src="image_4.png" /><figcaption>Fig. 5.1: Three types of singular value distributions used in the experiments.</figcaption></figure>

Figure 5.2 shows that for all three matrices and for both QRCP and randomized SVD, Algorithm 2.1 behaves as expected: the error $η_i$ is roughly equal to $ε_i^{i+1} = (θu_i)^{i+1}$ after $i$ refinement steps. Thus, using single precision as the low precision $u_i$, we can achieve an accuracy close to double precision with only one refinement step. Naturally, since for QRCP $θ = 2^{-3}$ is relatively small, $(u_i/θ)^2$ is noticeably larger than the accuracy achieved by a standard double precision LRA; the refined factors are thus not completely as accurate as if computed directly in double precision. This gap can be filled by performing a second refinement step, although this would likely be unnecessary in most practical scenarios.

Similar results are obtained using half precision as the low precision format (fp16 or bfloat16 arithmetic). An error close to the single precision accuracy can be achieved in just one or a few steps. Moreover, we can even reach double precision accuracy if needed, which illustrates an attractive property of Algorithm 2.1: it can reach an accuracy which is only limited by the high precision, but not by the low precision, despite this low precision being used to perform most of its operations.

Finally, we discuss the rank behavior of $F_i$ across IR steps. We see that it is roughly equal to $\text{RANK}(X, η_i)$, the numerical rank of $X$ at accuracy $η_i$. Thus, for matrices with rapidly decaying singular values such as the exponential case in Figure 5.1, these ranks tend to be much smaller at the early steps in low precision. Reflecting on the cost of the algorithm, we can expect the use of low precision plus refinement to be particularly cost efficient for such matrices.

Figure 5.3 shows similar plots for the tensor decompositions (TTSVD, HOSVD, and HTSVD). The results for tensors follow the same trend, and lead to the same conclusion as for matrices. The main difference is that the rank of $F_i$ (text labels) is now a vector instead.

Overall, our experiments confirm that Algorithm 2.1 is able to rapidly converge to a high accuracy, while using low precision for the LRA kernel. This observation is valid for all types of matrices and tensors in our test set, and all five LRA kernels that we tested, which shows that Algorithm 2.1 is robust and can work in a wide variety of settings.
<figure><img src="image_6.png" /><figcaption>Fig. 5.2: Convergence of Algorithm 2.1 for three types of matrices (with different singular value distributions, see Figure 5.1) and for two different LRA kernels (QRCP or randomized SVD).</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 5.2: Convergence of Algorithm 2.1 for three types of matrices (with different singular value distributions, see Figure 5.1) and for two different LRA kernels (QRCP or randomized SVD). The number next to each marker indicates the rank of $F_i$ after recompression.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>(a) TTSVD (exponential).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(b) HOSVD (exponential).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(c) TTSVD (power).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(d) HOSVD (power).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(e) TTSVD (linear).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(f) HOSVD (linear).</figcaption></figure>

Fig. 5.3: Convergence of Algorithm 2.1 for three types of tensors (depending on the singular value distribution, see Figure 5.1) and with TTSVD or HOSVD as LRA kernel.
<figure><img src="image_6.png" /><figcaption>(a) TTSVD (exponential).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(b) HOSVD (exponential).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(c) TTSVD (power).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(d) HOSVD (power).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(e) TTSVD (linear).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(f) HOSVD (linear).</figcaption></figure>

Fig. 5.3: Convergence of Algorithm 2.1 for three types of tensors (depending on the singular value distribution, see Figure 5.1) and with TTSVD or HOSVD as LRA kernel. The numbers next to each marker indicate the rank vector of $F_i$ after recompression.
<figure><img src="image_13.png" /><figcaption>Fig. 5.4: Convergence of Algorithm 2.1 for three types of tensors (depending on the singular value distribution, see Figure 5.1) and with HTSVD as LRA kernel.</figcaption></figure>

the error $\eta_i$ is smaller for larger values of $\theta_i$. For TTSVD (Table 5.1), this faster convergence is achieved without compromising the correct rank detection, which remains contained throughout the iterations and for all values of $\theta_i \le 1$. Thus, in this case, a large value of $\theta_i$ is recommended. The situation is different for QRCP (Table 5.2), for which a too large value of $\theta_i$ (here $\theta_i \ge 2^{-2}$) prevents the rank to be correctly detected, thus leads to a rank explosion.

We therefore conclude that the optimal choice of $\theta_i$ depends on the LRA algorithm, and more specifically, on its sensitivity to rounding errors when detecting the numerical rank. Empirically, we have observed QRCP to be the most sensitive of the LRA kernels, and to a lesser extent the randomized SVD kernel; the other kernels behaved well even for large $\theta_i$. For this reason, in our experiments we have set $\theta_i = 2^{-1}$ for all kernels except QRCP and randomized SVD, for which we have set $\theta_i = 2^{-3}$ and $\theta_i = 2^{-2}$, respectively. This setting allows to keep the ranks contained except for a few sporadic cases when using half precision.

**6. Conclusion.** We have presented a new mixed precision iterative refinement algorithm for computing low-rank matrix and tensor approximations. The algorithm first computes a low-rank approximation in low precision, and then computes another low-rank approximation of the error term, also in low precision, to refine the accuracy
<figure><img src="image_12.png" /><figcaption>Fig. 5.4: Convergence of Algorithm 2.1 for three types of tensors (depending on the singular value distribution, see Figure 5.1) and with HTSVD as LRA kernel. The numbers next to each marker indicate the rank vector of $F_i$ after recompression.</figcaption></figure>

**5.3. Role of $\theta_\ell$.** There is a tradeoff in choosing the scaling factor $\theta_\ell$: the larger it is, the faster the convergence (since the error is reduced by a factor roughly equal to $\theta_\ell u_\ell$ at each refinement step), but if it is too large, the rank will no longer be correctly detected and this will lead to a significant rank growth.

We illustrate this in Table 5.1 and Table 5.2, where we compare the convergence of Algorithm 2.1 for different values of $\theta_\ell$. Table 5.1 is for the TTSVD kernel and Table 5.2 is for the QRCP kernel; fp16 arithmetic is used as the low precision $u_\ell$ in both cases. The tables show that the method converges faster as $\theta_\ell$ increases, as expected: the error $\eta_i$ is smaller for larger values of $\theta_\ell$. For TTSVD (Table 5.1), this faster convergence is achieved without compromising the correct rank detection, which remains contained throughout the iterations and for all values of $\theta_\ell \le 1$. Thus, in this case, a large value of $\theta_\ell$ is recommended. The situation is different for QRCP (Table 5.2), for which a too large value of $\theta_\ell$ (here $\theta_\ell \ge 2^{-2}$) prevents the rank from being correctly detected, thus leading to a rank explosion.

We therefore conclude that the optimal choice of $\theta_\ell$ depends on the LRA algorithm, and more specifically, on its sensitivity to rounding errors when detecting the
<table><caption>Table 5.1: Relative error $\eta_i$ and RANK($F_i$) at different steps $i$ and for different values of $\theta_\ell$, for TTSVD (using fp16 as $u_\ell$ and exponential distribution of singular values).</caption><thead><tr><th rowspan="2">i \ &theta;<sub>l</sub></th><th colspan="5">Relative error &eta;<sub>i</sub></th><th colspan="5">RANK(F<sub>i</sub>)</th></tr><tr><th>2<sup>0</sup></th><th>2<sup>-1</sup></th><th>2<sup>-2</sup></th><th>2<sup>-3</sup></th><th>2<sup>-4</sup></th><th>2<sup>0</sup></th><th>2<sup>-1</sup></th><th>2<sup>-2</sup></th><th>2<sup>-3</sup></th><th>2<sup>-4</sup></th></tr></thead><tbody><tr><td>0</td><td>6e-04</td><td>8e-04</td><td>9e-04</td><td>2e-03</td><td>5e-03</td><td>23,39,16</td><td>07,07,06</td><td>06,06,06</td><td>05,05,05</td><td>04,04,04</td></tr><tr><td>1</td><td>3e-07</td><td>1e-06</td><td>4e-06</td><td>1e-05</td><td>9e-05</td><td>14,14,15</td><td>12,12,12</td><td>11,11,11</td><td>10,10,10</td><td>08,08,08</td></tr><tr><td>2</td><td>2e-10</td><td>1e-09</td><td>8e-09</td><td>7e-08</td><td>6e-07</td><td>21,21,21</td><td>19,19,19</td><td>17,17,17</td><td>15,15,15</td><td>13,13,13</td></tr><tr><td>3</td><td>9e-14</td><td>8e-13</td><td>2e-11</td><td>4e-10</td><td>4e-09</td><td>29,29,29</td><td>26,26,26</td><td>23,23,23</td><td>20,20,20</td><td>18,18,18</td></tr><tr><td>4</td><td></td><td>6e-14</td><td>6e-14</td><td>1e-12</td><td>4e-11</td><td></td><td>28,28,28</td><td>28,28,28</td><td>26,26,26</td><td>22,22,22</td></tr><tr><td>5</td><td></td><td></td><td></td><td>6e-14</td><td>3e-13</td><td></td><td></td><td></td><td>28,28,28</td><td>27,27,27</td></tr><tr><td>6</td><td></td><td></td><td></td><td></td><td>6e-14</td><td></td><td></td><td></td><td></td><td>28,28,28</td></tr></tbody></table>

<table><caption>Table 5.2: Relative error $\eta_i$ and RANK($F_i$) at different steps $i$ and for different values of $\theta$, for QRCP decomposition (using fp16 as $u_\ell$ and exponential distribution of singular values).</caption><thead><tr><th rowspan="2">i \ &theta;<sub>l</sub></th><th colspan="5">Relative error &eta;<sub>i</sub></th><th colspan="5">RANK(F<sub>i</sub>)</th></tr><tr><th>2<sup>0</sup></th><th>2<sup>-1</sup></th><th>2<sup>-2</sup></th><th>2<sup>-3</sup></th><th>2<sup>-4</sup></th><th>2<sup>0</sup></th><th>2<sup>-1</sup></th><th>2<sup>-2</sup></th><th>2<sup>-3</sup></th><th>2<sup>-4</sup></th></tr></thead><tbody><tr><td>0</td><td>2e-03</td><td>2e-03</td><td>2e-03</td><td>2e-03</td><td>3e-03</td><td>13</td><td>9</td><td>9</td><td>7</td><td>6</td></tr><tr><td>1</td><td>2e-05</td><td>1e-05</td><td>1e-05</td><td>5e-05</td><td>1e-04</td><td>25</td><td>18</td><td>13</td><td>10</td><td>9</td></tr><tr><td>2</td><td>5e-08</td><td>8e-08</td><td>7e-08</td><td>1e-07</td><td>1e-06</td><td>72</td><td>43</td><td>20</td><td>16</td><td>14</td></tr><tr><td>3</td><td>2e-10</td><td>2e-10</td><td>4e-10</td><td>8e-10</td><td>8e-09</td><td>95</td><td>86</td><td>36</td><td>21</td><td>19</td></tr><tr><td>4</td><td>2e-12</td><td>2e-12</td><td>3e-12</td><td>2e-12</td><td>4e-11</td><td>27</td><td>27</td><td>29</td><td>27</td><td>24</td></tr><tr><td>5</td><td>7e-13</td><td>7e-13</td><td>7e-13</td><td>7e-13</td><td>1e-12</td><td>28</td><td>28</td><td>28</td><td>28</td><td>27</td></tr><tr><td>6</td><td></td><td></td><td></td><td></td><td>7e-13</td><td></td><td></td><td></td><td></td><td>28</td></tr></tbody></table>

of the approximation. The process can be repeated to further refine the accuracy, and we ensure the rank of the approximation remains bounded by using inexpensive recompression operations. We have performed the error analysis of this algorithm, which proves that the low precision determines the convergence speed, whereas the attainable accuracy only depends on the high precision. Therefore, any desired level of accuracy can be attained, even though most of the operations are performed in low precision. This makes the algorithm computationally attractive, and we have performed a complexity analysis to determine specific conditions under which we can expect it to be cheaper than simply computing a low-rank approximation directly in high precision. We have applied our algorithm to various matrix and tensor low-rank approximations algorithms, and performed MATLAB experiments that confirm its robustness and convergence in a wide range of settings.

This paper lays the theoretical and algorithmic foundations for a mixed precision iterative refinement framework for matrix and tensor low-rank approximations. In future work, we will develop high performance implementations of these algorithms to tackle the low-rank approximation of large scale matrices and tensors on modern parallel architectures, especially those equipped with GPU accelerators, which provide extremely fast low precision arithmetic.

**7. Acknowledgements.** This work was performed thanks to generous grants from ANR (JCJC SELESTE, ANR-20-CE46-008-01) and Paris Ile-de-France Region (DIM RFSI RC-TENSOR No 2021-05).
<table><caption>Table 5.1: Relative error $\eta_i$ and RANK($F_i$) at different steps $i$ and for different values of $\theta_\ell$, for TTSVD (using fp16 as $u_\ell$ and exponential distribution of singular values).</caption><thead><tr><th rowspan="2">i|θ<sub>ℓ</sub></th><th colspan="5">Relative error η<sub>i</sub></th><th colspan="5">RANK(F<sub>i</sub>)</th></tr><tr><th>2<sup>0</sup></th><th>2<sup>-1</sup></th><th>2<sup>-2</sup></th><th>2<sup>-3</sup></th><th>2<sup>-4</sup></th><th>2<sup>0</sup></th><th>2<sup>-1</sup></th><th>2<sup>-2</sup></th><th>2<sup>-3</sup></th><th>2<sup>-4</sup></th></tr></thead><tbody><tr><td>0</td><td>6e-04</td><td>8e-04</td><td>9e-04</td><td>2e-03</td><td>5e-03</td><td>23,39,16</td><td>7, 7, 6</td><td>6, 6, 6</td><td>5, 5, 5</td><td>4, 4, 4</td></tr><tr><td>1</td><td>3e-07</td><td>1e-06</td><td>4e-06</td><td>1e-05</td><td>9e-05</td><td>14,14,15</td><td>12,12,12</td><td>11,11,11</td><td>10,10,10</td><td>8, 8, 8</td></tr><tr><td>2</td><td>2e-10</td><td>1e-09</td><td>8e-09</td><td>7e-08</td><td>6e-07</td><td>21,21,21</td><td>19,19,19</td><td>17,17,17</td><td>15,15,15</td><td>13,13,13</td></tr><tr><td>3</td><td>9e-14</td><td>8e-13</td><td>2e-11</td><td>4e-10</td><td>4e-09</td><td>29,29,29</td><td>26,26,26</td><td>23,23,23</td><td>20,20,20</td><td>18,18,18</td></tr><tr><td>4</td><td></td><td>6e-14</td><td>6e-14</td><td>1e-12</td><td>4e-11</td><td></td><td>28,28,28</td><td>28,28,28</td><td>26,26,26</td><td>22,22,22</td></tr><tr><td>5</td><td></td><td></td><td></td><td>6e-14</td><td>3e-13</td><td></td><td></td><td></td><td>28,28,28</td><td>27,27,27</td></tr><tr><td>6</td><td></td><td></td><td></td><td></td><td>6e-14</td><td></td><td></td><td></td><td></td><td>28,28,28</td></tr></tbody></table>

<table><caption>Table 5.2: Relative error $\eta_i$ and RANK($F_i$) at different steps $i$ and for different values of $\theta$, for QRCP decomposition (using fp16 as $u_\ell$ and exponential distribution of singular values).</caption><thead><tr><th rowspan="2">i|θ<sub>ℓ</sub></th><th colspan="5">Relative error η<sub>i</sub></th><th colspan="5">RANK(F<sub>i</sub>)</th></tr><tr><th>2<sup>0</sup></th><th>2<sup>-1</sup></th><th>2<sup>-2</sup></th><th>2<sup>-3</sup></th><th>2<sup>-4</sup></th><th>2<sup>0</sup></th><th>2<sup>-1</sup></th><th>2<sup>-2</sup></th><th>2<sup>-3</sup></th><th>2<sup>-4</sup></th></tr></thead><tbody><tr><td>0</td><td>2e-03</td><td>2e-03</td><td>2e-03</td><td>2e-03</td><td>3e-03</td><td>13</td><td>9</td><td>9</td><td>7</td><td>6</td></tr><tr><td>1</td><td>2e-05</td><td>1e-05</td><td>1e-05</td><td>5e-05</td><td>1e-04</td><td>25</td><td>18</td><td>13</td><td>10</td><td>9</td></tr><tr><td>2</td><td>5e-08</td><td>8e-08</td><td>7e-08</td><td>1e-07</td><td>1e-06</td><td>72</td><td>43</td><td>20</td><td>16</td><td>14</td></tr><tr><td>3</td><td>2e-10</td><td>2e-10</td><td>4e-10</td><td>8e-10</td><td>8e-09</td><td>95</td><td>86</td><td>36</td><td>21</td><td>19</td></tr><tr><td>4</td><td>2e-12</td><td>2e-12</td><td>3e-12</td><td>2e-12</td><td>4e-11</td><td>27</td><td>27</td><td>29</td><td>27</td><td>24</td></tr><tr><td>5</td><td>7e-13</td><td>7e-13</td><td>7e-13</td><td>7e-13</td><td>1e-12</td><td>28</td><td>28</td><td>28</td><td>28</td><td>27</td></tr><tr><td>6</td><td></td><td></td><td></td><td></td><td>7e-13</td><td></td><td></td><td></td><td></td><td>28</td></tr></tbody></table>

numerical rank. Empirically, we have observed QRCP to be the most sensitive of the LRA kernels, and to a lesser extent the randomized SVD kernel; the other kernels behaved well even for large $\theta_\ell$. For this reason, in our experiments we have set $\theta_\ell = 2^{-1}$ for all kernels except QRCP and randomized SVD, for which we have set $\theta_\ell = 2^{-3}$ and $\theta_\ell = 2^{-2}$, respectively. This setting allows to keep the ranks contained except for a few sporadic cases when using half precision. It must be noted that the sensitivity to rounding errors is also matrix (or tensor) dependent; however, we do not tune $\theta_\ell$ individually for each matrix/tensor, because this would not be representative of a realistic scenario where it is unknown in advance how sensitive the data at hand actually is.

**5.4. Results on real-life data.** Finally, we also experiment our method on data sets (one matrix and one tensor) coming from real-life applications. The matrix corresponds to a subdomain in the discretization of the Poisson equation. It is of size $253 \times 252$ and its singular values decay relatively rapidly (not shown). The tensor comes from a quantum chemistry application and corresponds to an eigenfunction of the Hamiltonian operator for computing the vibrational spectrum of the H2CO molecule [30]. Here, the tensor is 6th-order ($d = 6$), of size $17 \times 17 \times 13 \times 13 \times 9 \times 9$, and obtained from an eigenvalue computation in which both the operator and vectors are compressed in the TT format throughout (and we obtain the full tensor with DECOMPRESS in high precision in the end), and the precision is set to $10^{-9}$ to prevent rank explosion.

We report the results for these two problems in Figure 5.5-Figure 5.6 with all
[1] E. AGULLO, O. COULAUD, L. GIRAUD, M. IANNACITO, G. MARAIT, AND N. SCHENKELS, *The backward stable variants of GMRES in variable accuracy*, (2022), https://inria.hal.science/ hal-03776837.

[2] P. AMESTOY, O. BOITEAU, A. BUTTARI, M. GEREST, F. JÉZÉQUEL, J.-Y. L'EXCELLENT, AND T. MARY, *Mixed Precision Low Rank Approximations and their Application to Block Low Rank LU Factorization*, IMA J. Numer. Anal., (2022), https://doi.org/10.1093/imanum/ drac037.

[3] B. W. BADER AND T. G. KOLDA, *Algorithm 862: MATLAB tensor classes for fast algorithm prototyping*, ACM Trans. Math. Software, 32 (2006), pp. 635–653, https://doi.org/10.1145/ 1186785.1186794.

[4] P. BLANCHARD, N. J. HIGHAM, F. LOPEZ, T. MARY, AND S. PRANESH, *Mixed precision block fused multiply-add: Error analysis and application to GPU tensor cores*, SIAM J. Sci. Comput., 42 (2020), pp. C124–C141, https://doi.org/10.1137/19M1289546.

[5] A. CICHOCKI, N. LEE, I. OSELEDETS, A.-H. PHAN, Q. ZHAO, AND D. P. MANDIC, *Tensor networks for dimensionality reduction and large-scale optimization: Part 1 low-rank tensor decompositions*, Foundations and Trends® in Machine Learning, 9 (2016), pp. 249–429, https://doi.org/http://dx.doi.org/10.1561/2200000059.

[6] P. COMON, *Tensors: a brief introduction*, IEEE Signal Processing Magazine, 31 (2014), pp. 44–53, https://doi.org/10.1109/MSP.2014.2298533.

[7] M. P. CONNOLLY, N. J. HIGHAM, AND S. PRANESH, *Randomized low rank matrix approximation: Rounding error analysis and a mixed precision algorithm*, MIMS EPrint 2022.5, Manchester Institute for Mathematical Sciences, The University of Manchester, UK, July 2022, http: //eprints.maths.manchester.ac.uk/2884/. Revised March 2023.

[8] O. COULAUD, L. GIRAUD, AND M. IANNACITO, *A robust gmres algorithm in tensor train format*, arXiv preprint arXiv:2210.14533, (2022), https://doi.org/10.48550/arXiv.2210.14533.

[9] L. DE LATHAUWER, B. DE MOOR, AND J. VANDEWALLE, *A multilinear singular value decomposition*, SIAM J. Matrix Anal. Appl., 21 (2000), pp. 1253–1278, https://doi.org/10.1137/ S0895479896305696.

[10] C. ECKART AND G. YOUNG, *The approximation of one matrix by another of lower rank*, Psychometrika, 1 (1936), pp. 211–218, https://doi.org/10.1007/BF02288367.

[11] M. FASI, N. J. HIGHAM, F. LOPEZ, T. MARY, AND M. MIKAITIS, *Matrix multiplication in multiword arithmetic: Error analysis and application to GPU tensor cores*, SIAM J. Sci. Comput., 45 (2023), pp. C1–C19, https://doi.org/10.1137/21m1465032.

[12] G. H. GOLUB AND C. F. VAN LOAN, *Matrix computations*, JHU press, 2013, https://doi.org/ 10.1137/1028073.

[13] L. GRASEDYCK, *Hierarchical singular value decomposition of tensors*, SIAM J. Matrix Anal. Appl., 31 (2010), pp. 2029–2054, https://doi.org/10.1137/090764189.

[14] L. GRASEDYCK, D. KRESSNER, AND C. TOBLER, *A literature survey of low-rank tensor approximation techniques*, GAMM-Mitteilungen, 36 (2013), pp. 53–78, https://doi.org/10.1002/ gamm.201310004.

[15] W. HACKBUSCH AND S. KÜHN, *A new scheme for the tensor representation*, Journal of Fourier analysis and applications, 15 (2009), pp. 706–722, https://doi.org/10.1007/ s00041-009-9094-9.

[16] N. J. HIGHAM, *Accuracy and Stability of Numerical Algorithms*, Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, second ed., 2002, https://doi.org/10.1137/ 1.9780898718027.

[17] N. J. HIGHAM AND T. MARY, *Mixed precision algorithms in numerical linear algebra*, Acta Numerica, 31 (2022), pp. 347–414, https://doi.org/10.1017/s0962492922000022.

[18] N. J. HIGHAM AND S. PRANESH, *Simulating low precision floating-point arithmetic*, SIAM J. Sci. Comput., 41 (2019), pp. C585–C602, https://doi.org/10.1137/19M1251308.

[19] N. KISHORE KUMAR AND J. SCHNEIDER, *Literature survey on low rank approximation of matrices*, Linear and Multilinear Algebra, 65 (2017), pp. 2212–2244, https://doi.org/10.1080/ 03081087.2016.1267104.

[20] T. G. KOLDA AND B. W. BADER, *Tensor decompositions and applications*, SIAM Rev., 51 (2009), pp. 455–500, https://doi.org/10.1137/07070111X.

[21] D. KRESSNER AND C. TOBLER, *htucker—a MATLAB toolbox for tensors in hierarchical tucker format*, Mathicse, EPF Lausanne, (2012).

[22] P.-G. MARTINSSON AND S. VORONIN, *A randomized blocked algorithm for efficiently computing rank-revealing factorizations of matrices*, SIAM J. Sci. Comput., 38 (2016), pp. S485–S507, https://doi.org/10.1137/15M1026080.
the corresponding LRA kernels, using the same legend and methodology as for the random data experiments. Overall, the results lead to the same trends and conclusions as the previous experiments; ranks progressively increase with the IR steps, and our method manages to achieve high precision levels in all cases. In the TTSVD case, ranks stop increasing around the compute precision of the application ($10^{-9}$), whereas in HOSVD and HTSVD cases we still observe a significant rank increase beyond this threshold. This is mostly due to the fact that the underlying tensor decomposition of this tensor is TT by construction, and this topology mismatch results in higher ranks.

<figure><img src="image_13.png" /><figcaption>Fig. 5.5: Convergence of Algorithm 2.1 for a real-life matrix (Poisson) of size 253 × 252, for different LRA kernels. The numbers next to each marker indicate the rank of $F_i$ after recompression.</figcaption></figure>

**5.5. Estimation of the time cost and role of $\omega_\ell$.** Based on the results obtained in the previous experiments, we can also obtain a rough estimation of the time cost of each method, using the formulas (2.31) (for classical hardware) or (2.34) (for hardware with block FMA units). This estimation depends on the number of iterations needed to converge to the desired accuracy, the ranks at each iteration, and the relative speed of the low precision arithmetic $\omega_\ell$.

We illustrate how this estimation can be obtained by using the case of randomized SVD with the Poisson matrix (Figure 5.5b) as an example. We consider two possible settings:

* Setting 1: classical (CPU) hardware, where the relative speed of each arithmetic is proportional to the number of bits. For this setting we target an accuracy of $\epsilon = 10^{-12}$ and we use fp64 as the high precision and fp32 as the low precision (hence $\omega_\ell = 0.5$).

* Setting 2: GPU hardware with tensor core units (a type of block FMA). We target an accuracy of $\epsilon = 10^{-5}$ and we use fp32 as the high precision and fp16 as the low precision. On most tensor core architectures (A100 and H100 in particular), fp16 arithmetic is 16 times faster than fp32, hence we take $\omega_\ell = \frac{1}{16} = 0.0625$.

In setting 1, our IR method (red curve in Figure 5.5b) converges in 2 iterations to an accuracy of $\epsilon = 10^{-12}$. We use formula (2.31) with $\text{RANK}(X, 10^{-6}) = 14$,
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848

[23] R. OOI, T. IWASHITA, T. FUKAYA, A. IDA, AND R. YOKOTA, *Effect of mixed precision computing on H-matrix vector multiplication in BEM analysis*, in Proceedings of the International Conference on High Performance Computing in Asia-Pacific Region, 2020, pp. 92–101, https://doi.org/10.1145/3368474.3368479.

[24] H. OOTOMO AND R. YOKOTA, *Recovering single precision accuracy from tensor cores while surpassing the FP32 theoretical peak performance*, Int. J. High Perform. Comput. Appl., 36 (2022), pp. 475–491, https://doi.org/10.1177/10943420221090256.

[25] H. OOTOMO AND R. YOKOTA, *Mixed-precision random projection for RandNLA on tensor cores*, 2023, https://arxiv.org/abs/2304.04612.

[26] I. OSELEDETS, *TT-Toolbox v2.2.2*, 2023, https://github.com/oseledets/TT-Toolbox.

[27] I. V. OSELEDETS, *Tensor-train decomposition*, SIAM J. Sci. Comput., 33 (2011), pp. 2295–2317, https://doi.org/10.1137/090752286.

[28] L. R. TUCKER, *Some mathematical notes on three-mode factor analysis*, Psychometrika, 31 (1966), pp. 279–311, https://doi.org/10.1007/BF02289464.

[29] N. VANNIEUWENHOVEN, R. VANDEBRIL, AND K. MEERBERGEN, *A new truncation strategy for the higher-order singular value decomposition*, SIAM J. Sci. Comput., 34 (2012), pp. A1027–A1052, https://doi.org/10.1137/110836067.

[30] Z. YANG, J. SHAN, AND Z. ZHANG, *Hardware-efficient mixed-precision CP tensor decomposition*, arXiv preprint arXiv:2209.04003, (2022), https://doi.org/10.48550/arXiv.2209.04003.
<figure><img src="image_14.png" /><figcaption>Fig. 5.6: Convergence of Algorithm 2.1 for a real-life tensor (H2CO) of size 17 × 17 × 13 × 13 × 9 × 9, for different LRA kernels. The numbers next to each marker indicate the rank of $F_i$ after recompression.</figcaption></figure>

$\text{RANK}(X, 10^{-12}) = 42$, $c_1 = 6$, and $c_2 = 2$ to obtain an estimated cost of

$$p \times (0.5 \times 6 \times (14 + 14 + 42) + 2 \times 14 + 3) = 241p.$$

This is to be compared with simply computing the LRA directly in fp64 arithmetic, for which formula (2.32) gives an estimated time cost of $6 \times 42 \times p = 252p$. Hence, we can expect a time reduction of about 4%. This example illustrates that even with classical hardware where the low precision is “only” twice as fast as the high precision, our IR method can still achieve moderate speedups in the case where the ranks decrease rapidly at low accuracies, such as for this Poisson matrix. This is however not the case for other matrices, including the random synthetic ones used previously.

In setting 2, our IR method (yellow curve in Figure 5.5b) converges in 2 iterations to an accuracy of $\epsilon = 10^{-5}$. We use formula (2.34) with $\text{RANK}(X, 10^{-3}) = 5$, $\text{RANK}(X, 10^{-5}) = 11$, $c_1 = 6$, and $c_2 = 2$ to obtain an estimated cost of

$$p \times (0.0625 \times (6 \times (5 + 5 + 11) + 2 \times 5) + 3) = 11.5p.$$

This is to be compared with simply computing the LRA directly in fp32 arithmetic,
for which formula (2.32) gives an estimated time cost of $6 \times 11 \times p = 66p$. Hence we can expect a significant time reduction by a factor of almost $6 \times$. This example illustrates that our IR method is particularly attractive on hardware with very fast low precision arithmetic.

The above examples are only meant to illustrate the trends that one might expect on different types of hardware and for different target accuracies and low precisions. This analysis can be similarly applied to the other matrices, LRA algorithms, and to tensors.

**6. Conclusion.** We have presented a new mixed precision iterative refinement algorithm for computing low-rank matrix and tensor approximations. The algorithm first computes a low-rank approximation in low precision, and then computes another low-rank approximation of the error term, also in low precision, to refine the accuracy of the approximation. The process can be repeated to further refine the accuracy, and we ensure the rank of the approximation remains bounded by using inexpensive recompression operations. We have performed the error analysis of this algorithm, which proves that the low precision determines the convergence speed, whereas the attainable accuracy only depends on the high precision. Therefore, high accuracy can be attained, even though most of the operations are performed in low precision. This makes the algorithm computationally attractive, and we have performed a complexity analysis to determine specific conditions under which we can expect it to be cheaper than simply computing a low-rank approximation directly in high precision. We have applied our algorithm to various matrix and tensor low-rank approximations algorithms, and performed MATLAB experiments that confirm its robustness and convergence in a wide range of settings.

This paper lays the theoretical and algorithmic foundations for a mixed precision iterative refinement framework for matrix and tensor low-rank approximations. In future work, we will develop high performance implementations of these algorithms to tackle the low-rank approximation of large scale matrices and tensors on modern parallel architectures, especially those equipped with GPU accelerators, which provide extremely fast low precision arithmetic. We have presented preliminary performance results in [3], which show that Algorithm 2.1 with randomized SVD, using fp16 as the low precision and fp32 as the high precision, can achieve up to a $2.2 \times$ speedup with respect to randomized SVD in fp32, while achieving an accuracy that is equivalent or better.

Our error analysis is based on a general model that assumes the use of numerically stable implementations for the basic kernels (in particular, low-rank approximation and recompression). For tensors, such computations are experimentally observed to behave stably in practice, but formally proving their stability remains an open problem which provides an interesting perspective for future work.

**Acknowledgements.** This work was performed thanks to generous grants from ANR (JCJC SELESTE, ANR-20-CE46-0008-01) and Paris Ile-de-France Region (DIM RFSI RC-TENSOR No 2021-05). We thank Olivier Coulaud for providing the H2CO tensor data.
REFERENCES

[1] E. AGULLO, O. COULAUD, L. GIRAUD, M. IANNACITO, G. MARAIT, AND N. SCHENKELS, *The backward stable variants of GMRES in variable accuracy*, Research Report RR-9483, Inria, Sept. 2022, https://inria.hal.science/hal-03776837.

[2] P. R. AMESTOY, O. BOITEAU, A. BUTTARI, M. GEREST, F. JÉZÉQUEL, J.-Y. L'EXCELLENT, AND T. MARY, *Mixed precision low rank approximations and their application to block low rank LU factorization*, IMA J. Numer. Anal., (2022), https://doi.org/10.1093/imanum/drac037.

[3] M. BABOULIN, S. DONFACK, O. KAYA, T. MARY, AND M. ROBEYNS, *Mixed precision randomized low-rank approximation with GPU tensor cores*, in Euro-Par 2024: Parallel Processing, Cham, 2024, Springer Nature Switzerland, pp. 31–44, https://doi.org/10.1007/978-3-031-69583-4_3.

[4] B. W. BADER AND T. G. KOLDA, *Algorithm 862: MATLAB tensor classes for fast algorithm prototyping*, ACM Trans. Math. Software, 32 (2006), pp. 635–653, https://doi.org/10.1145/1186785.

[5] C. BELTRÁN, V. NOFERINI, AND N. VANNIEUWENHOVEN, *When can forward stable algorithms be composed stably?*, IMA J. Numer. Anal., 44 (2024), pp. 886–919, https://doi.org/10.48550/arXiv.2109.10610.

[6] P. BLANCHARD, N. J. HIGHAM, F. LOPEZ, T. MARY, AND S. PRANESH, *Mixed precision block fused multiply-add: Error analysis and application to GPU tensor cores*, SIAM J. Sci. Comput., 42 (2020), pp. C124–C141, https://doi.org/10.1137/19M1289546.

[7] A. CICHOCKI, N. LEE, I. OSELEDETS, A.-H. PHAN, Q. ZHAO, AND D. P. MANDIC, *Tensor networks for dimensionality reduction and large-scale optimization: Part 1 low-rank tensor decompositions*, Foundations and Trends® in Machine Learning, 9 (2016), pp. 249–429, https://doi.org/http://dx.doi.org/10.1561/2200000059.

[8] P. COMON, *Tensors: a brief introduction*, IEEE Signal Processing Magazine, 31 (2014), pp. 44–53, https://doi.org/10.1109/MSP.2014.2298533.

[9] M. P. CONNOLLY, N. J. HIGHAM, AND S. PRANESH, *Randomized low rank matrix approximation: Rounding error analysis and a mixed precision algorithm*, MIMS EPrint 2022.5, Manchester Institute for Mathematical Sciences, The University of Manchester, UK, July 2022, http://eprints.maths.manchester.ac.uk/2884/. Revised March 2023.

[10] O. COULAUD, L. GIRAUD, AND M. IANNACITO, *A robust GMRES algorithm in tensor train format*, arXiv preprint arXiv:2210.14533, (2022).

[11] L. DE LATHAUWER, B. DE MOOR, AND J. VANDEWALLE, *A multilinear singular value decomposition*, SIAM J. Matrix Anal. Appl., 21 (2000), pp. 1253–1278, https://doi.org/10.1137/S0895479896305696.

[12] C. ECKART AND G. YOUNG, *The approximation of one matrix by another of lower rank*, Psychometrika, 1 (1936), pp. 211–218, https://doi.org/10.1007/BF02288367.

[13] M. FASI, N. J. HIGHAM, F. LOPEZ, T. MARY, AND M. MIKAITIS, *Matrix multiplication in multiword arithmetic: Error analysis and application to GPU tensor cores*, SIAM J. Sci. Comput., 45 (2023), pp. C1–C19, https://doi.org/10.1137/21m1465032.

[14] G. H. GOLUB AND C. F. VAN LOAN, *Matrix computations*, JHU press, 2013, https://doi.org/10.1137/1.9781421407944.

[15] L. GRASEDYCK, *Hierarchical singular value decomposition of tensors*, SIAM J. Matrix Anal. Appl., 31 (2010), pp. 2029–2054, https://doi.org/10.1137/090764189.

[16] L. GRASEDYCK, D. KRESSNER, AND C. TOBLER, *A literature survey of low-rank tensor approximation techniques*, GAMM-Mitteilungen, 36 (2013), pp. 53–78, https://doi.org/10.1002/gamm.201310003.

[17] W. HACKBUSCH AND S. KÜHN, *A new scheme for the tensor representation*, Journal of Fourier analysis and applications, 15 (2009), pp. 706–722, https://doi.org/10.1007/s00041-009-9094-9.

[18] N. J. HIGHAM, *Accuracy and Stability of Numerical Algorithms*, Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, second ed., 2002, https://doi.org/10.1137/1.9780898718027.

[19] N. J. HIGHAM AND T. MARY, *Mixed precision algorithms in numerical linear algebra*, Acta Numerica, 31 (2022), pp. 347–414, https://doi.org/10.1017/s0962492922000022.

[20] N. J. HIGHAM AND S. PRANESH, *Simulating low precision floating-point arithmetic*, SIAM J. Sci. Comput., 41 (2019), pp. C585–C602, https://doi.org/10.1137/19M1251308.

[21] N. KISHORE KUMAR AND J. SCHNEIDER, *Literature survey on low rank approximation of matrices*, Linear Multilinear Algebra, 65 (2017), pp. 2212–2244, https://doi.org/10.1080/03081087.2016.1267104.

[22] T. G. KOLDA AND B. W. BADER, *Tensor decompositions and applications*, SIAM Rev., 51
(2009), pp. 455–500, https://doi.org/10.1137/07070111X.

[23] D. KRESSNER AND C. TOBLER, *Algorithm 941: htucker—a MATLAB toolbox for tensors in hierarchical tucker format*, ACM Trans. Math. Softw., (2014), https://doi.org/10.1145/2538688.

[24] P.-G. MARTINSSON AND S. VORONIN, *A randomized blocked algorithm for efficiently computing rank-revealing factorizations of matrices*, SIAM J. Sci. Comput., 38 (2016), pp. S485–S507, https://doi.org/10.1137/15M1026080.

[25] R. OOI, T. IWASHITA, T. FUKAYA, A. IDA, AND R. YOKOTA, *Effect of mixed precision computing on H-matrix vector multiplication in BEM analysis*, in Proceedings of the International Conference on High Performance Computing in Asia-Pacific Region, HPCAsia '20, New York, NY, USA, 2020, Association for Computing Machinery, p. 92–101, https://doi.org/10.1145/3368474.3368479.

[26] H. OOTOMO AND R. YOKOTA, *Recovering single precision accuracy from Tensor Cores while surpassing the fp32 theoretical peak performance*, Int. J. High Perform. Comput. Appl., 36 (2022), pp. 475–491, https://doi.org/10.1177/10943420221090256.

[27] H. OOTOMO AND R. YOKOTA, *Mixed-precision random projection for RandNLA on tensor cores*, 2023, https://doi.org/10.1145/3592979.3593413.

[28] I. OSELEDETS, *TT-Toolbox v2.2.2*, 2023, https://github.com/oseledets/TT-Toolbox.

[29] I. V. OSELEDETS, *Tensor-train decomposition*, SIAM J. Sci. Comput., 33 (2011), pp. 2295–2317, https://doi.org/10.1137/090752286.

[30] M. RAKHUBA AND I. OSELEDETS, *Calculating vibrational spectra of molecules using tensor train decomposition*, The Journal of Chemical Physics, 145 (2016), p. 124101, https://doi.org/10.1063/1.4962420.

[31] L. R. TUCKER, *Some mathematical notes on three-mode factor analysis*, Psychometrika, 31 (1966), pp. 279–311, https://doi.org/10.1007/BF02289464.

[32] N. VANNIEUWENHOVEN, R. VANDEBRIL, AND K. MEERBERGEN, *A new truncation strategy for the higher-order singular value decomposition*, SIAM J. Sci. Comput., 34 (2012), pp. A1027–A1052, https://doi.org/10.1137/110836067.

[33] Z. YANG, J. SHAN, AND Z. ZHANG, *Hardware-efficient mixed-precision CP tensor decomposition*, arXiv preprint arXiv:2209.04003, (2022), https://doi.org/10.48550/arXiv.2209.04003.
Appendix A. Randomized SVD-based Recompress algorithm.

To recompress a (non-optimal) low-rank matrix $X = U_{in}Σ_{in}V_{in}^T$ to an optimal representation $U_{out}Σ_{out}V_{out}^T$ using randomized SVD, Algorithm 3.2 must be adapted as outlined in Algorithm A.1. The matrix products on Lines 6 and 8 of Algorithm A.1 can be efficiently performed by exploiting the low-rank structure of $X$. We stop the main loop when the norm of the low-rank matrix $Q_E R_E$ becomes smaller than the desired accuracy $ε$. This $Q_E R_E$ matrix is initially equal to the input matrix $X$ and is incrementally deflated at each step by $Q_b B_b$. In order to measure its norm, we reorthonormalize $Q_E$ at each step and measure the norm of $R_E$.

**Algorithm A.1 RECOMPRESS** algorithm using randomized SVD decomposition.

**Input:** a truncated SVD $U_{in}Σ_{in}V_{in}^T ∈ R^{m×n}$, the target accuracy $ε$, and a block size $b$.

**Output:** a recompressed SVD $U_{out}Σ_{out}V_{out}^T$ decomposition.

1. Initialize $Q$ and $B$ to empty matrices

2. Initialize $Q_E = U_{in}$ and $R_E = Σ_{in}V_{in}^T$.

3. $n_X = ||Σ_{in}||$

4. **repeat**

5. Draw a random Gaussian matrix $Ω ∈ R^{n×b}$.

6. $Y = U_{in}Σ_{in}(V_{in}^TΩ)$

7. $Q_b = qr(Y - Q(Q^T Y))$

8. $B_b = (Q_b^T U_{in})Σ_{in}V_{in}^T$

9. $Q ← [Q Q_b]$

10. $B ← \begin{bmatrix} B \\ B_b \end{bmatrix}$

11. $Q_E ← [Q_E Q_b]$

12. $R_E ← \begin{bmatrix} R_E \\ -B_b \end{bmatrix}$

13. $Q_E \bar{R} = qr(Q_E)$

14. $R_E ← \bar{R} R_E$

15. **until** $||R_E|| ≤ εn_X$

16. Compute the truncated SVD decomposition $B ≈ UΣ_{out}V_{out}^T$ at accuracy $ε$.

17. $U_{out} = Q \bar{U}$