# CROQuant: Complex Rank-One Quantization Algorithm

Maël Chaumette, Rémi Gribonval, Elisa Riccietti

► To cite this version:

Maël Chaumette, Rémi Gribonval, Elisa Riccietti. CROQuant: Complex Rank-One Quantization Algorithm. GRETSI 2025 – XXXème Colloque Francophone de Traitement du Signal et des Images, Aug 2025, Strasbourg, France. hal-05140371

HAL Id: hal-05140371

https://hal.science/hal-05140371v1

Submitted on 2 Jul 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# CROQuant: Complex Rank-One Quantization Algorithm

Maël CHAUMETTE<fnref n="1" /> Rémi GRIBONVAL<fnref n="1" /> Elisa RICCIETTI<fnref n="2" />

**Résumé** – Cet article présente un algorithme de quantification pour les matrices de rang un à valeurs complexes qui exploite les invariances par remise à l'échelle du problème pour obtenir un meilleur résultat que la quantification au plus proche voisin. En s'appuyant sur cet algorithme, on propose aussi une approche pour la quantification des matrices creuses à structure papillon à valeurs complexes apparaissant par exemple dans la transformée de Fourier rapide. Comparé à la quantification au plus proche voisin élément par élément, on obtient ainsi une réduction de 30% du nombre de bits nécessaires pour une précision donnée sur les matrices papillons, tout en maintenant une complexité polynomiale en la dimension des matrices.

**Abstract** – This paper presents a quantization algorithm for complex-valued rank-one matrices that exploits rescaling-invariances of the problem to obtain better results than round-to-nearest strategy. This algorithm is also used as a building block for an heuristic strategy to quantize complex-valued butterfly-structured sparse matrices appearing for example in the fast Fourier transform. Compared to element-wise round-to-nearest quantization, the number of bits is reduced by 30% the number of bits for a given precision on butterfly matrices, while maintaining a polynomial time complexity in the dimension of the matrices.

## 1 Introduction

Quantization is a fundamental issue in computer science and machine learning. With the ever-increasing size of models and data, the need to reduce memory by lowering numerical precision while maintaining acceptable performance is becoming crucial. In particular, quantization reduces the memory required to store matrices and speeds up the computations. These advantages are widely needed for deep neural networks that are mainly composed of huge, dense matrices.

Rank-one matrices play a central role in many compression problems in machine learning and signal processing, e.g., the Singular Value Decomposition [9]. Recently, an optimal quantization approach for real-valued rank-one matrices exploiting rescaling-invariances to minimize the quantization error, was introduced in [8]. In this work, we extend these results to the case of complex-valued matrices, an essential generalization to treat several applications in signal processing, such as the Fast Fourier Transform (FFT). Indeed, the FFT involves a butterfly structure, corresponding to a product of structured sparse matrices often used to factorize dense matrices, such as Hadamard and discrete cosine transform matrices [10], to speed up matrix-vector products and reduce memory space [12]. Quantization of such matrices further amplifies these computational gains. Their particular structure allows quantization to be heuristically decomposed into a succession of rank-one quantization problems exploiting our approach.

In Section 2, we formalize the quantization problem for complex rank-one matrices and explain why this problem cannot be solved simply with the optimal quantization algorithm applied to the real and imaginary parts.

We propose in Section 3 an algorithm adapted to this com-

plex setting. We show that the properties of invariance by rescaling persist in the complex setting, and allow for more accurate quantization than a naive nearest-neighbor rounding approach. Our algorithm depends on a parameter $b_m$ that controls the accuracy and affects the computational time (the higher $b_m$ the higher the accuracy but the longer the computational time). However, we empirically show that small values of this parameter give good results. This algorithm is also much faster than the brute-force method of testing all possible combinations for each element of the rank-one matrix: the latter provides an algorithm with exponential time complexity in $t$, the number of significant bits, and in the dimensions $m, n$ of the matrix; whereas our algorithm has a time complexity of $O(nm \min(n, m)b_m^2 2^{2t})$, which is simply polynomial in the dimension of the rank-one matrix. The exponential dependence on $t$ is tractable in this context because we are interested in small values of $t$, typically $t = 4$ or 6 for applications in modern float formats typically used in machine learning.

In Section 4, we discuss the quantization of complex butterfly matrices that play a central role in the FFT. Our algorithm provides quantization that is also more efficient than the naive rounding approach. At a given accuracy, our algorithm will need 30% fewer bits than the naive rounding method.

**Notations.** $\mathbb{F}_t$ is a set of floating-point numbers with $t$-bit of significance, as in [8]. We denote vectors in lowercase boldface ($\mathbf{x}$) and matrices in uppercase boldface letters ($\mathbf{X}$). $\|\cdot\|$ is the Frobenius norm. Re, Im denote the real and imaginary parts of a complex number, which we quantize separately. We therefore define $\mathbb{C}\mathbb{F}_t := \mathbb{F}_t + i\mathbb{F}_t$.

## 2 Problem formulation and baselines

Given $\mathbf{x} \in \mathbb{C}^m$ and $\mathbf{y} \in \mathbb{C}^n$, we formulate the quantization problem as the following minimization problem:

$$ \hat{\mathbf{x}}, \hat{\mathbf{y}} \in \arg \min_{\hat{\mathbf{x}} \in \mathbb{C}^{m}, \hat{\mathbf{y}} \in \mathbb{C}^{n}} \underbrace{\|\mathbf{x}\mathbf{y}^{\mathrm{H}} - \hat{\mathbf{x}}\hat{\mathbf{y}}^{\mathrm{H}}\|^{2}}_{C_{\infty, \mathbf{y}}(\hat{\mathbf{x}}, \hat{\mathbf{y}})}. \quad (1) $$

<footnote n="1">¹Inria, CNRS, ENS de Lyon, Université Claude Bernard Lyon 1, LIP, UMR 5668, 69342, Lyon cedex 07, France</footnote>

<footnote n="2">²ENS de Lyon, CNRS, Inria, Université Claude Bernard Lyon 1, LIP, UMR 5668, 69342, Lyon cedex 07, France</footnote>

<footnote>This project was supported by the SHARP project of the PEPR-IA (ANR-23-PEIA-0008, granted by France 2030), and the project ANR MEPHISTO ANR-24-CE23-7039. The authors would like to thank the Centre Blaise Pascal's platform at ENS de Lyon (Lyon, France) for computing facilities. The platform uses the SIDUS solution [14] developed by Emmanuel Quemener.</footnote>
Like in [8], the naive approach, called round-to-nearest (RTN), involves mapping each element of $x$ and $y$ to its nearest neighbor in $\mathbb{C}F_t$. We define the function round($\cdot$) that maps the real and the imaginary parts of a complex number to their nearest neighbor in $\mathbb{F}_t$. For a vector, the round function is applied element-wise. This method does not take into account that the problem is invariant by rescaling: $\forall \lambda \in \mathbb{C}^*, xy^H = (\lambda x)(\frac{1}{\lambda}y)^H$. Another possible method is to apply separately the optimal algorithm of [8] to the real and imaginary parts of $x$ and $y$. However, this method is not optimal either, since it omits cross-terms in the product of the complex numbers. We therefore need to find another method to obtain optimal quantization of complex rank-one matrices.

### 3 Proposed quantization algorithm

In this section, after characterizing the optimal solution to the quantization problem (1), we build an algorithm providing an efficient approximate solution, discuss its complexity, and finally present its performance against the RTN method.

**Characterizing the optimal solution.** The following lemma is a key result to address (1).

**Lemma 1.** Given $x \in \mathbb{C}^m$, $y \in \mathbb{C}^n$ and $t \ge 1$, it holds

$$ \inf_{\hat{x} \in \mathbb{C}F_t^m, \hat{y} \in \mathbb{C}F_t^n} \|x\hat{y}^H - \hat{x}\hat{y}^H\|^2 = \inf_{\lambda \in \mathbb{C}^n} f(\lambda) \quad (2) $$

$$ \text{with } f(\lambda) := \max_{\hat{x} \in \text{round}(\lambda x)} \|x\hat{y}^H - \hat{x} \text{round}(\mu(\hat{x})y)^H\|^2 \quad (3) $$

$$ \text{where } \mu(\hat{x}) := \begin{cases} \frac{\langle x, \hat{x} \rangle}{\|\hat{x}\|^2} & \text{if } \hat{x} \neq 0 \\ 0 & \text{otherwise} \end{cases} \quad (4) $$

The lemma and its proof are the complex-valued analog of [8, Lemma 4.2], that was shown to allow reducing the 4mn-variable optimization problem (1) to a one-variable problem: to find an optimal $\lambda^*$ of $f$ defined in (3), assuming it exists.

**Using the piecewise constant nature of $f$.** Since the function $\lambda \mapsto \text{round}(\lambda x)$ is piecewise constant, $f$ is also piecewise constant, so finding a minimizer on $\mathbb{C}$ amounts to finding on which of the corresponding pieces $f$ attains the smallest value. The simplest way to individuate these pieces is to study the *breakpoints* of the function $\lambda \in \mathbb{C} \mapsto \text{round}(\lambda x)$ where $x \in \mathbb{C}^m$. These breakpoints correspond to $\lambda$ values for which there is at least one entry $x_j$ such that $\text{round}(\lambda x_j)$ corresponds to a tie in the choice of the nearest neighbor for the real or the imaginary part.

**Characterization of breakpoints through breaklines.** The following result describes these breakpoints.

**Lemma 2.** Consider $x \in \mathbb{C}^m$. For each entry $x_j := u + iv$, $j = 1, \dots, m$, the breakpoints of $\lambda \in \mathbb{C} \mapsto \text{round}(\lambda x_j)$ are straight lines in the complex plane with equations

$$ \begin{cases} u \operatorname{Im}(\lambda) = -v \operatorname{Re}(\lambda) + s(k + \frac{1}{2})2^{-b-t} \\ v \operatorname{Im}(\lambda) = u \operatorname{Re}(\lambda) + s(k + \frac{1}{2})2^{-b-t} \end{cases} \quad (5) \\ \forall k \in [2^{t-1}, 2^t - 1], \forall s \in \{-1, +1\}, \forall b \in \mathbb{Z} $$

The proof is an adaptation of the proof in [8, Lemma 5.1]. We call *breaklines* the lines defined by (5). There are infinitely many such breaklines approximating “accumulation lines” defined by taking $b \to +\infty$, with equations $u \operatorname{Im}(\lambda) + v \operatorname{Re}(\lambda) = 0$ (resp. $v \operatorname{Im}(\lambda) - u \operatorname{Re}(\lambda) = 0$).

<figure><img src="image_2.png" /><figcaption>Figure 1 – Function $\lambda \in \mathbb{C} \mapsto \frac{\|x\hat{y}^H - \hat{x}(\lambda)\hat{y}(\lambda)^H\|}{\|x\hat{y}^H\|}$ with $\hat{x}(\lambda) := \text{round}(\lambda x)$, $\hat{y}(\lambda) := \text{round}(\mu(\hat{x}(\lambda))y)$ and $\mu(\hat{x}(\lambda))$ defined in (4). Breaklines of the function $\lambda \in \Omega(x) \mapsto \text{round}(\lambda x)$ in black (corresponding to $b \le 5$ in Lemma 2), centroids in green, accumulation lines in orange and the centroids on the accumulation lines in magenta. The search space $\Omega(x)$ is delimited by the red lines. Here, $x, y \in \mathbb{C}^2$ were drawn with real and imaginary parts following a uniform distribution in $[0, 1]$ and $t = 3$.</figcaption></figure>

**Minimizing $f$: restricting to a compact domain.** Furthermore, it is easy to show that $f(2\lambda) = f(i\lambda) = f(\lambda)$ for every $\lambda \in \mathbb{C}$, so studying function $f$ can be limited to any compact domain $\Omega \subset \mathbb{C}$ generating a tiling of the complex plane via dilations and quarter turn rotation (i.e., with the sets $i^k 2^j \Omega$, $0 \le k \le 3$, $j \in \mathbb{Z}$). Figure 1 illustrates that such a domain can be chosen to be a trapezium with sides defined by breaklines of $f$. We denote it $\Omega(x)$. Figure 1 also shows the typical shape of $f$ on $\Omega(x)$ with its breaklines (corresponding to $b \le 5$), and accumulation lines displayed in orange.

**Construction of the algorithm.** As in the real-valued case, since $f$ is piecewise constant, minimizing it amounts to find on which of the *polygonal pieces* (delimited by its breaklines) it attains its smallest value. One also needs to consider the (constant) values of $f$ on finitely many *line segments* lying on the intersection of the accumulation lines with $\Omega(x)$. However, contrary from the real-valued case where it suffices to visit a finite number of breakpoints, here in the complex-valued case there are infinitely many pieces even in the bounded region $\Omega(x)$. It is therefore not possible to explicitly evaluate $f$ in all these pieces in order to find $\hat{x}^*$ and $\hat{y}^*$. For this reason we define a parameter $b_m$, and choose to consider just the breaklines corresponding to $b \le b_m$ in Lemma 2. This defines a finite number of pieces, which we can enumerate and iterate over using the python library shapely<fnref n="1" />. Each of the resulting pieces has a centroid, and the collection of such centroids, completed with the midpoints of the line segments described above (displayed in magenta on Figure 1) is denoted $C_{b_m}(x)$. Algorithm 1 describes the resulting near-minimization algorithm.

**Role of the parameter.** The larger the parameter $b_m$, the closer the resulting pair $(\hat{x}, \hat{y})$ is to being optimal, but this will have an impact on time complexity because this will increase the number of tested centroids. We will see below that $b_m = 1$ is enough to get better results than the RTN approach.

<footnote n="1">1. https://shapely.readthedocs.io/en/stable/</footnote>
**Algorithm 1: Complex rank-one quantization algorithm**

**Data:** $x \in \mathbb{C}^m$, $y \in \mathbb{C}^n$, $t \ge 1$, $b_m \in \mathbb{Z}$.

**Result:** $\hat{x}^* \in \mathbb{C}F_t^m$, $\hat{y}^* \in \overline{\mathbb{C}}F_t^n$, $\lambda^* \in \mathbb{C}$.

1. Initialize $\hat{x}^* \leftarrow 0, \hat{y}^* \leftarrow 0$;

2. **if** $x, y \ne 0$ **then**

3. Build set of centroids $C_{b_m}(x)$ (Lemma 2, $\lambda$ values on accumulation lines and call to shapely);

4. **for** $\lambda \in C_{b_m}(x)$ **do**

5. $\hat{x} \leftarrow \text{round}(\lambda x)$;

6. $\mu \leftarrow \frac{\langle x, \hat{x} \rangle}{\|\hat{x}\|^2}$;

7. $\hat{y} \leftarrow \text{round}(\mu y)$;

8. **if** $C_{x,y}(\hat{x}, \hat{y}) < C_{x,y}(\hat{x}^*, \hat{y}^*)$ **then**

9. $\hat{x}^* \leftarrow \hat{x}, \hat{y}^* \leftarrow \hat{y}, \lambda^* \leftarrow \lambda$;

<figure><img src="image_3.png" /><figcaption>Figure 2 – Scatterplot of $\rho_{b_m}$, $\rho_{rtn}$ for 100 pairs $(x, y) \in \mathbb{C}^n \times \mathbb{C}^n$ with uniformly distributed real and imaginary parts in $[0, 1]$ for different values of $n$ and $t$ with $b_m = 1$.</figcaption></figure>

**Complexity.** The complexity of our algorithm is determined in several stages. Firstly, the number of breaklines is bounded by $4mb_m2^{t-1}$. Then, according to [6], the maximum number of regions defined by $k$ straight lines is $\frac{k(k+1)}{2} + 1$. In our case, we can conclude that there are $O(m^2b_m^22^{2t})$ centroids. For each centroid, we need to calculate $C_{x,y}(\hat{x}, \hat{y})$, which has a cost of $O(m+n)$, by exploiting the same trick as in the real-valued case [8]. The total cost of the algorithm is therefore $O((m+n)m^2b_m^22^{2t})$. But, since $x$ and $y$ play symmetrical roles, we can swap their roles by looping on $C_{b_m}(y)$, which will result in a time complexity of $O((m+n)n^2b_m^22^{2t})$. Thus Algorithm 1 has a time complexity $O(nm \min(n, m)b_m^22^{2t})$.

**Experiments and results.** Experiments are run on an Intel(R) Xeon (R) Gold 5218 CPU @ 2.30GHz. Our implementation is available in open source [1]. We consider 100 pairs $(x, y) \in \mathbb{C}^{12} \times \mathbb{C}^{12}$, where the real and imaginary parts of each component follow a uniform distribution on $[0, 1]$. We set $t = 4$. Table 1 shows the average of the relative quantization error of our algorithm $\rho_{b_m} := \frac{\|xy^H - \hat{x}^*(\hat{y}^*)^H\|}{\|xy^H\|}$ and of the computation time for different values of $b_m$. The computation time increases rapidly, while the relative variation between the error at $b_m = 1$ and $b_m = 6$ is around 5%. This empirical proof encourages us to use $b_m = 1$ in the rest of our experiments.

<table><caption>Table 1 – Average of the quantization error and the of calculation time of Algorithm 1 over 100 pairs $(x, y) \in \mathbb{C}^{12} \times \mathbb{C}^{12}$ with uniformly distributed real and imaginary parts in $[0, 1]$ for different values of $b_m$ with $t = 4$.</caption><thead><tr><th></th><th>b<sub>m</sub> = -2</th><th>b<sub>m</sub> = 1</th><th>b<sub>m</sub> = 3</th><th>b<sub>m</sub> = 6</th></tr></thead><tbody><tr><td>Error (×10<sup>-2</sup>)</td><td>2.398</td><td>2.308</td><td>2.202</td><td>2.173</td></tr><tr><td>Time (s)</td><td>0.913</td><td>22.46</td><td>64.07</td><td>161.5</td></tr></tbody></table>

We compare the quantization error $\rho_{b_m}$ of our algorithm and the RTN strategy error $\rho_{rtn} := \frac{\|xy^H - \text{round}(x) \text{round}(y)^H\|}{\|xy^H\|}$ in Figure 2 with the same setup. As in the real-valued case [8], the performance of Algorithm 1 depends on $n$, unlike that of RTN, but it is more efficient especially for small $n$. Indeed, when $n = 4$, $\rho_{b_m}$ is approximately 5 times smaller than $\rho_{rtn}$. Furthermore, as $n$ increases, the point clouds become more concentrated.

We also highlight that the naive application of the optimal real-valued algorithm to the real and imaginary parts of $x$ and $y$ (results not displayed) yields a quantization error approximately 10 times larger than $\rho_{b_m}$.

# 4 Application to butterfly quantization

We now exploit our algorithm to quantize butterfly matrices similar to those appearing in the FFT. The FFT accelerates the discrete Fourier transform (DFT) using the Cooley-Tukey algorithm [2]. The main idea is to use the divide-and-conquer strategy, which recursively factorizes the DFT of size $n \times n$ into $\log_2(n)$ sparse matrices (for simplicity we consider here a dimension $n$ that is a power of two, but more flexible butterfly factorizations can be defined [11]), the butterfly matrices. If we denote $F \in \mathbb{C}^{n \times n}$ the DFT matrix, then its factorization associated to the Cooley-Tukey algorithm is

$$F = B_1 \cdots B_L$$

where $L := \log_2(n)$ and the $B_i$ are structured sparse matrices with the so-called Kronecker-sparse structure [11].

To quantize the factors $B_i$, in the spirit of what was done in the real-valued case [8], we aim to leverage Algorithm 1 to work on a product of complex-valued butterfly factors.

**Quantization algorithm for butterfly matrices.** First, let us take a detour into optimal quantization for a product of two (possibly sparse) matrices $X, Y \in \mathbb{C}^{n \times r}$: the quantization problem then reads

$$X^*, Y^* \in \arg \min_{\hat{X}, \hat{Y} \in \mathcal{CF}_t} \|XY^H - \hat{X}\hat{Y}^H\|^2 \quad (6)$$

where $\mathcal{CF}_t$ is the set of matrices with coefficient in $\mathbb{C}F_t$ and the same support as $X, Y$. The support of a matrix $A$ is the set of coordinates $(i, j)$ where $A_{i,j} \ne 0$. To the best of our knowledge, there is no method for finding the optimal quantization for this problem. However, if the rank-one matrices, $x_i y_i^*$, where $x_i$ and $y_i$ are the corresponding columns of $X$ and $Y$, $1 \le i \le r$, have disjoint supports, then the problem (6) can be decomposed into $r$ independent complex rank-one matrix quantization problems:

$$\|XY^H - \hat{X}\hat{Y}^H\|^2 = \sum_{i=1}^{r} \|x_i y_i^H - \hat{x}_i \hat{y}_i^H\|^2. \quad (7)$$

Each sub-problem can be addressed using Algorithm 1.

We can leverage the two-factor setting: the crux is that when considering certain so-called *chainable* Kronecker-sparse factors [11], for any subset of consecutive factors, the
<figure><img src="image_4.png" /><figcaption>Figure 3 – Quantization error on the butterfly decomposition of 10 random matrices as a function of $t$ ($n = 256$, $b_m = 1$).</figcaption></figure>

product between $X = B_{l_0} \cdots B_{l_1} \in \mathbb{C}^{n \times n}$ and $Y^H = B_{l_1+1} \cdots B_{l_2} \in \mathbb{C}^{n \times n}$, with $1 \le l_0 \le l_1 \le l_2 \le L$, can be written as a sum of $n$ rank-one matrices with disjoint support [10]. As suggested for the real-valued case [8], this can be used to heuristically decompose the product of $L$ butterfly matrices into several products of two matrices to apply the optimal quantization of the problem (6). The quantization of each product will therefore be optimal, but the global quantization will not. As in [8] we consider two heuristics:

- pairwise: writing $(B_1B_2)(B_3B_4)\cdots(B_{L-1}B_L)$, quantization is applied to each pair of consecutive factors (if $L$ is odd, the last matrix is quantized using the RTN strategy).

- Left-to-Right (LTR): writing $B_1(B_2(\cdots(B_{L-1}B_L)))$, quantization is performed from left to right (details in [8]).

We also implement stochastic rounding [4] and fixed-point rounding [5], which work like RTN, i.e., quantization is applied element-wise without any rescaling.

**Evaluation on random butterfly matrices.** To compare these methods, we generate 10 random collections of Kronecker-sparse factors $B_1, \dots, B_L$ and quantize them with the methods described before to obtain $\hat{B}_1 \cdots \hat{B}_L$. Figure 3 shows the evolution of the quantization error $\rho := \frac{\|\hat{B}_1 \cdots \hat{B}_L - \hat{B}_1 \cdots \hat{B}_L\|}{\|\hat{B}_1 \cdots \hat{B}_L\|}$ for different values of $t$ and $b_m = 1$. Despite the matrix product parenthesis heuristic, pairwise and LTR are more efficient than RTN, stochastic and fixed-point. Indeed, we see that LTR and pairwise need about 4 bits to achieve an error close to $10^{-2}$ while RTN, fixed-point and stochastic must have more than 6 bits to achieve this error. In general, according to the exponential fit, pairwise and LTR need $1 - \frac{1}{4} \approx 30\%$ fewer bits than RTN, fixed-point and stochastic to achieve the same quantization error.

**Evaluation on the quantization of the FFT.** We also apply the five quantization algorithms to the (exact) Cooley-Tukey factorization $F = B_1 \cdots B_L \in \mathbb{C}^{n \times n}$ where $F$ is the DFT and look at the error made if we apply the quantized version of $F$ to a vector $x$ instead of applying the original high-resolution DFT.

Let $x \in \mathbb{R}^n$ be the signal and $y := Fx \in \mathbb{C}^n$ its Fourier transform. We also define $\hat{y} := \hat{F}x = \hat{B}_1 \cdots \hat{B}_L x$. The average of $\rho_{\text{fft}} := \frac{\|y - \hat{y}\|}{\|y\|}$ to over 10 standard Gaussian signals, with the five considered algorithms are displayed in Table 2 with $n = 256$, $t = 5$ and $b_m = 1$. Like for the butterfly quantization error, the LTR and pairwise methods are more efficient than the RTN, stochastic and fixed-point quantizations.

<table><caption>Table 2 – $\rho_{\text{fft}}$ with different quantization strategies for $n = 256$, $t = 5$ and $b_m = 1$.</caption><thead><tr><th></th><th>LTR</th><th>pairwise</th><th>RTN</th><th>Stochastic</th><th>Fixed</th></tr></thead><tbody><tr><td>ρ<sub>fft</sub> × 10<sup>-2</sup></td><td>0.697</td><td>0.293</td><td>2.353</td><td>2.747</td><td>3.263</td></tr></tbody></table>

# 5 Conclusion

In this paper, we proposed an efficient quantization algorithm exploiting rescaling-invariances for complex-valued rank-one matrices. We showed that this approach enables more accurate quantization than a naive element-wise rounding. Applying this algorithm to butterfly matrices, which play a key role in many fast transforms such as the Fast Fourier Transform, we demonstrated the possibility to reduce quantization error for a given number of bits, or alternatively to reduce by 30% the required number of bits for a given precision.

Our algorithm depends on a parameter that controls the accuracy but also impacts the computation time. Empirically, this parameter is not very limiting, but since the real-valued case leads to a fully optimal algorithm with bounded complexity, a natural challenge is to understand whether this remains possible in the complex-valued case via further mathematical analysis steps. Furthermore, throughout this paper, we have neglected underflow and overflow [3, 13]. A similar study taking these phenomena into account would be interesting.

This work opens the way to several perspectives: the first is the quantization of the product of matrices of any rank. The application to butterfly matrices leads us to consider applying our algorithm to other matrix decompositions, such as the Toeplitz matrices that models, for example, the covariance matrix of some time series [7]. Another challenge is to extend this work to quantize ReLU networks, which satisfy similar rescaling-invariances.

# References

[1] M. Chaumette, R. Gribonval, and E. Riccietti. Code for reproducible research - CROQuant: Complex Rank-One Quantization Algorithm. https://inria.hal.science/hal-05133647, June 2025.

[2] J. W. Cooley and J. W. Tukey. An algorithm for the machine calculation of complex Fourier series. *Math. Comp.*, 19(90):297–301, 1965.

[3] A. Cuyt, P. Kuterna, B. Verdonk, and D. Verschaeren. Underflow revisited. *Calcolo*, 39:169–179, 2002.

[4] E.-M. El Arar, D. Sohier, P. de Oliveira Castro, and E. Petit. Stochastic rounding variance and probabilistic bounds: A new approach. *SISC*, 45(5):C255–C275, 2023.

[5] M. Gerken. On fixed-point quantization schemes. In *38th MWSCAS*, volume 1, pages 350–353 vol.1, 1995.

[6] R. Graham, D. Knuth, and O. Patashnik. Concrete mathematics, a foundation for computer science. *Math. Gaz.*, 75(471):117–119, 1991.

[7] R. M. Gray. On unbounded Toeplitz matrices and nonstationary time series with an application to information theory. *Information and control*, 24(2):181–196, 1974.

[8] R. Gribonval, T. Mary, and E. Riccietti. Optimal quantization of rank-one matrices in floating-point arithmetic—with applications to butterfly factorizations. June 2023. working paper or preprint on HAL.

[9] Q.-T. Le, E. Riccietti, and R. Gribonval. Spurious valleys, np-hardness, and tractability of sparse matrix factorization with fixed support. *SIMAX*, 44(2):503–529, 2023.

[10] Q.-T. Le, L. Zheng, E. Riccietti, and R. Gribonval. Fast learning of fast transforms, with guarantees. In *ICASSP*, pages 3348–3352. IEEE, 2022.

[11] Q.-T. Le, L. Zheng, E. Riccietti, and R. Gribonval. Butterfly factorization with error guarantees, November 2024. arXiv:2411.04506.

[12] L. Le Magoarou and R. Gribonval. Flexible multilayer sparse approximations of matrices and applications. *JSTSP*, 10(4):688–700, 2016.

[13] W. T. Padgett and D. V. Anderson. *Quantization Effects -Round-Off Noise and Overflow*, pages 83–111. Springer, Cham, 2009.

[14] E. Quemener and M. Corvellec. Sidus—the solution for extreme deduplication of an operating system. *Linux Journal*, 2013(235):3, 2013.