Fast approximation of orthogonal matrices and application to PCA

TL;DR

This paper introduces a fast approximation method for orthogonal matrices using extended Givens transformations, enabling efficient PCA computations with a good balance of accuracy and speed.

Contribution

It proposes a novel approximation technique for orthogonal matrices based on extended Givens transformations and an efficient greedy algorithm, improving computational speed for spectral methods.

Findings

The method achieves a good trade-off between accuracy and computational speed.

It effectively approximates orthogonal matrices for spectral methods like PCA.

Experimental results demonstrate improved efficiency in PCA applications.

Abstract

We study the problem of approximating orthogonal matrices so that their application is numerically fast and yet accurate. We find an approximation by solving an optimization problem over a set of structured matrices, that we call extended orthogonal Givens transformations, including Givens rotations as a special case. We propose an efficient greedy algorithm to solve such a problem and show that it strikes a balance between approximation accuracy and speed of computation. The approach is relevant to spectral methods and we illustrate its application to PCA.

Tables

Table 1. Table 1 : Average classification accuracies for k-NN when using PCA projections and our approximations without spectrum update, for various datasets. We also show the speedup (FLOPS and actual running time) and the number of features selected in the calculations as a proportion out of the total d 𝑑 d (see also Remark 2). Results are averaged over 100 random realizations (train/test splits).

DATASET	FULL PCA	PROPOSED ALGORITHM
	accuracy	accuracy	speedup(FLOPS)	speedup(TIME)	selection
PENDIGITS	95 $\pm$ 0.7	91 $\pm$ 2.1	$\times$1.6	$\times$1.1	1
ISOLET	92 $\pm$ 0.4	90 $\pm$ 1.0	$\times$12	$\times$10.1	1
USPS	95 $\pm$ 1.2	94 $\pm$ 0.8	$\times$7.7	$\times$4.7	0.64
UCI	90 $\pm$ 1.9	87 $\pm$ 1.5	$\times$2.5	$\times$1.6	0.72
20NEWS	80 $\pm$ 3.1	77 $\pm$ 2.1	$\times$3.1	$\times$2.5	0.3
EMNIST digits	97 $\pm$ 2.5	95 $\pm$ 1.8	$\times$13	$\times$11.3	0.37
MNIST 8m	96 $\pm$ 2.0	94 $\pm$ 0.9	$\times$15	$\times$13.7	0.28

Table 2. Table 2 : Speed-up achieved with the proposed approximations 𝐔 ¯ p subscript ¯ 𝐔 𝑝 \mathbf{\bar{U}}_{p} against a vanilla C implementation of dense matrix-vector multiplication and BLAS Level 2 functions, in both cases with no parallelism.

DATASET:	ISOLET	USPS	EMNIST digits	MNIST 8m
C language (TIME)	$\times$10.1	$\times$4.7	$\times$11.3	$\times$13.7
BLAS (TIME)	$\times$4.8	$\times$2.5	$\times$5.7	$\times$6.8

Equations

$$\underset{\bar{\mathbf{U}},{\bar{\mathbf{\Sigma }}}_p}{\mathrm{minimize}}\Vert \mathbf{U}{\mathbf{\Sigma }}_p-\bar{\mathbf{U}}{\bar{\mathbf{\Sigma }}}_p{\Vert }_F^2\text{ subject to }\bar{\mathbf{U}}\in {\mathcal{F}}_g,$$

$${\mathbf{G}}_{ij}=\left[\begin{array}{ccccc}{\mathbf{I}}_{i-1}& & & & \\ & \ast & & \ast & \\ & & {\mathbf{I}}_{j-i-1}& & \\ & \ast & & \ast & \\ & & & & {\mathbf{I}}_{d-j}\end{array}\right],{\tilde{\mathbf{G}}}_{ij}\in \left\{\left[\begin{array}{cc}c& -s\\ s& c\end{array}\right],\left[\begin{array}{cc}c& s\\ s& -c\end{array}\right]\right\},\text{ such that }{c}^2+s^2=1,$$

$$\bar{\mathbf{U}}=\prod _{k=1}^g{\mathbf{G}}_{i_k{j}_k}={\mathbf{G}}_{i_1{j}_1}\dots {\mathbf{G}}_{i_g{j}_g}.$$

$$\Vert \mathbf{U}{\mathbf{\Sigma }}_p-{\bar{\mathbf{U}}\bar{\mathbf{\Sigma }}}_p{\Vert }_F^2=\Vert \underset{{\mathbf{L}}^{(k)}}{\underbrace{\prod _{t=1}^{k-1}{\mathbf{G}}_{i_t{j}_t}^T\mathbf{U}{\mathbf{\Sigma }}_p}}-{\mathbf{G}}_{i_k{j}_k}\underset{{\mathbf{N}}^{(k)}}{\underbrace{\prod _{t=k+1}^g{\mathbf{G}}_{i_t{j}_t}{\bar{\mathbf{\Sigma }}}_p}}{\Vert }_F^2=\Vert {\mathbf{L}}^{(k)}-{\mathbf{G}}_{i_k{j}_k}{\mathbf{N}}^{(k)}{\Vert }_F^2.$$

$$\ C_{ij}=\|\mathbf{Z}_{\{i,j\}}\|_{*}-\text{tr}(\mathbf{Z}_{\{i,j\}})=\left\{\begin{array}[]{lr}\sqrt{(Z_{ii}+Z_{jj})^{2}+(Z_{ij}-Z_{ji})^{2}}-Z_{ii}-Z_{jj},\text{if }\det(\mathbf{Z}_{\{i,j\}})\geq 0\\ \sqrt{(Z_{ii}-Z_{jj})^{2}+(Z_{ij}+Z_{ji})^{2}}-Z_{ii}-Z_{jj},\text{if }\det(\mathbf{Z}_{\{i,j\}})<0\end{array}\right.$$

$$(i^{\star },j^{\star })=\underset{(i,j),j>i}{\mathrm{argmax}}{C}_{ij}.$$

$$\mathbb{E}[\Vert \mathbf{U}-\bar{\mathbf{U}}{\Vert }_F^2]\le 2d-\sqrt{2\pi d}.$$

$$\mathbb{E}\left[\Vert \mathbf{U}-\bar{\mathbf{U}}{\Vert }_F^2\right]\le 2(d-\lfloor r\rfloor )-\frac{2\sqrt{2}}{\sqrt{\pi }}\sqrt{d-\lfloor r\rfloor },\text{ where }r=d-\frac{1+\sqrt{(2d-1{)}^2-8g}}{2}.$$

$$\text{off}({\mathbf{UG}}_{ij}^T{)}^2\le \text{off}(\mathbf{U}{)}^2+\frac{1}{2}((U_{ii}-U_{jj}{)}^2-(U_{ij}-U_{ji}{)}^2).$$

$$\mathbb{E}[C_{ij}(u)]\propto (1-u{)}^2,$$

$${\beta }_i=\arccos ({\tau }_i),i=1,\dots ,p,$$

$${\mathbf{XX}}^T\approx {\mathbf{U}}_p({\mathbf{\Sigma }}_p{\mathbf{\Sigma }}_p^T){\mathbf{U}}_p^T\text{ and }\mathbf{X}\approx {\mathbf{U}}_p{\mathbf{\Sigma }}_p{\mathbf{V}}_p^T.$$

$$\Vert {\mathbf{U}}_p{\mathbf{\Sigma }}_p{\mathbf{V}}_p^T-\bar{\mathbf{U}}{\bar{\mathbf{\Sigma }}}_p{\mathbf{V}}_p^T{\Vert }_F^2=\Vert {\mathbf{U}}_p{\mathbf{\Sigma }}_p-\bar{\mathbf{U}}{\bar{\mathbf{\Sigma }}}_p{\Vert }_F^2,$$

$$\underset{\bar{\mathbf{U}},\bar{\mathbf{\Lambda }}}{\mathrm{minimize}}\Vert {\mathbf{XX}}^T-\bar{\mathbf{U}}\bar{\mathbf{\Lambda }}\bar{\mathbf{U}}{\Vert }_F^2\text{ subject to }\bar{\mathbf{U}}\in {\mathcal{F}}_g.$$

$$\mathbf{U}=\left(\prod _{i=1}^d\prod _{j=i+1}^d{\mathbf{G}}_{ij}({\theta }_{ij})\right)\mathbf{D},$$

$$Z\propto \left(\prod _{k=2}^d{\cos }^{k-2}{\theta }_{1k}\right)\left(\prod _{k=3}^d{\cos }^{k-3}{\theta }_{2k}\right)\left(\prod _{k=d}^d{\cos }^{k-d}{\theta }_{(d-1)k}\right).$$

$$f(y,\alpha ,\beta )=\frac{y^{\alpha -1}(1-y{)}^{\beta -1}}{B(\alpha ,\beta )},B(\alpha ,\beta )=\frac{\Gamma (\alpha )\Gamma (\beta )}{\Gamma (\alpha +\beta )}\text{ with }y\in [0,1].$$

$$\Vert \mathbf{L}-{\mathbf{G}}_{ij}\mathbf{N}{\Vert }_F^2=\Vert \mathbf{L}{\Vert }_F^2+\Vert {\mathbf{G}}_{ij}\mathbf{N}{\Vert }_F^2-2\text{tr}({\mathbf{N}}^T{\mathbf{G}}_{ij}^T\mathbf{L})=\Vert {\sigma }_p{\Vert }_2^2+\Vert {\bar{\sigma }}_p{\Vert }_2^2-2\text{tr}({\mathbf{G}}_{ij}^T\mathbf{L}{\mathbf{N}}^T),$$

$$\text{tr}({\mathbf{G}}_{ij}^T\mathbf{L}{\mathbf{N}}^T)=\sum _{k=1,k\notin \{i,j\}}^d{Z}_{kk}+\text{tr}({\tilde{\mathbf{G}}}_{ij}^T{\mathbf{Z}}_{\{i,j\}})=\text{tr}(\mathbf{Z})+\text{tr}({\tilde{\mathbf{G}}}_{ij}^T{\mathbf{Z}}_{\{i,j\}})-\text{tr}({\mathbf{Z}}_{\{i,j\}}).$$

$${\tilde{\mathbf{G}}}_{ij}^{\star }={\mathbf{V}}_1{\mathbf{V}}_2^T,{\mathbf{Z}}_{\{i,j\}}={\mathbf{V}}_1{\mathbf{SV}}_2^T,$$

$$\text{max}\text{tr}({\mathbf{G}}_{ij}^T\mathbf{L}{\mathbf{N}}^T)=\text{tr}(\mathbf{Z})+\text{tr}(\mathbf{S})-\text{tr}({\mathbf{Z}}_{\{i,j\}})=\text{tr}(\mathbf{Z})+\Vert {\mathbf{Z}}_{\{i,j\}}{\Vert }_{\ast }-\text{tr}({\mathbf{Z}}_{\{i,j\}})=\text{tr}(\mathbf{Z})+C_{ij}.$$

$$C_{ij}=\Vert {\mathbf{Z}}_{\{i,j\}}{\Vert }_{\ast }-\text{tr}({\mathbf{Z}}_{\{i,j\}}).$$

$$\sum _{k=1}^{d/2}{C}_{i_k{j}_k}=\sum _{k=1}^{d/2}(\Vert {\mathbf{U}}_{\{i_k,j_k\}}{\Vert }_{\ast }-U_{i_k{i}_k}-U_{j_k{j}_k})=\sum _{k=1}^{d/2}\Vert {\mathbf{U}}_{\{i_k,j_k\}}{\Vert }_{\ast }-\text{tr}(\mathbf{U}).$$

$${\Vert \mathbf{U}-\prod _{k=1}^{d/2}{\mathbf{G}}_{i_k{j}_k}\Vert }_F^2=2d-2\text{tr}(\mathbf{U})-2\sum _{k=1}^{d/2}{C}_{i_k{j}_k}=2d-2\sum _{k=1}^{d/2}\Vert {\mathbf{U}}_{\{i_k{j}_k\}}{\Vert }_{\ast }.$$

$$|\text{tr}(\mathbf{U}_{\{i,j\}})|\!\!=\!\!|\text{tr}(\mathbf{V}_{1}\mathbf{SV}_{2}^{T})|\!\!=\!\!|\text{tr}(\mathbf{SV}_{2}^{T}\mathbf{V}_{1})|\!\!=\!\!|\text{tr}(\mathbf{S\Delta})|\!=\!\Biggl{\lvert}\sum_{t=1}^{2}\!s_{t}\Delta_{tt}\Biggr{\rvert}\!\leq\!\Delta_{\max}\!\sum_{t=1}^{2}s_{t}\!\!=\!\!\Delta_{\max}\|\mathbf{U}_{\{i,j\}}\|_{*},$$

$${\Delta }_{\max }=\max \{|{\Delta }_{11}|,|{\Delta }_{22}|\}.$$

$$\mathbb{E}[\Vert {\mathbf{U}}_{\{i_k,j_k\}}{\Vert }_{\ast }]\ge \frac{\pi }{2}\mathbb{E}[\text{tr}({\mathbf{U}}_{\{i_k,j_k\}})],$$

$$\mathbb{E}[{\Delta }_{\max }]=\frac{1}{{\pi }^2}{\iint }_0^{\pi }|\cos (x)\cos (y)+\sin (x)\sin (y)|dxdy=\frac{2}{\pi }.$$

$$\mathbb{E}\left[{\Vert \mathbf{U}-\prod _{k=1}^{d/2}{\mathbf{G}}_{i_k{j}_k}\Vert }_F^2\right]$$

$$\le 2d-\pi \mathbb{E}[\text{tr}(\mathbf{U})]=2d-\pi \mathbb{E}\left[\sum _{t=1}^d|U_{tt}|\right]=2d-\sqrt{2\pi d},$$

$${\mathbf{J}}_{i_g{j}_g}\dots {\mathbf{J}}_{i_1{j}_1}\mathbf{U}=\left[\begin{array}{cc}{\mathbf{D}}^{\prime \prime }& 0_{r\times (d-r)}\\ 0_{(d-r)\times r}& {\mathbf{U}}^{\prime }\end{array}\right],$$

Full text

1 Introduction

Orthonormal transformations play a key role in most matrix decomposition techniques and spectral methods [4]. As such, manipulating them in an efficient manner is essential to many practical applications. While general matrix-vector multiplications with orthogonal matrices take $O(d^{2})$ space and time, it is natural to ask whether faster approximate computations (say $O(d\log d)$) can be achieved while retaining enough accuracy.

Approximating an orthonormal matrix with just a few building blocks is hard in general. The standard decomposition technique meant to reduce complexity is a low-rank approximation. Unfortunately, for an orthonormal matrix, which is perfectly conditioned, this approach is meaningless.

In this work, we are inspired by the fact that several orthonormal/unitary transformations that exhibit low numerical complexity are known. The typical example is the discrete Fourier transform with its efficient implementation as the fast Fourier transform [47] together with other Fourier-related algorithms: fast Walsh-Hadamard transforms [17], fast cosine transforms [33], and fast Hartley transforms [9]. Other approaches include fast wavelet transforms [5], banded orthonormal matrices [40, 44] and fast Slepian transforms [27]. Decomposition of orthogonal matrices into $O(d)$ Householder reflectors or $O(d^{2})$ Givens rotations [20][Chapter 5.1] are known already. These basic building blocks have been further extended, for example we have fast Givens rotations [36] and two generalizations of the Givens rotations [6] and [35]. In theoretical physics, unitary decompositions parametrize symmetry groups [45], and they are compactly parametrized using $\sigma$-matrices [42] or symmetric positive definite matrices [3]. To the best of our knowledge, none of these factorizations focus on reducing the computational complexity of using the orthogonal/unitary transformations but rather they model properties of physical systems.

Our idea is to approximately factor any orthonormal matrix into a product of a fixed number of sparse matrices such that their application (to a vector) has linearithmic complexity. In this paper, we derive structured approximations to orthonormal matrices that can be found efficiently and applied remarkably fast. We pose the search for an efficient approximation as an optimization problem over a product of structured matrices, the extended orthogonal Givens transformations. These structures extend Givens rotations to also include reflectors with no computational drawback and suggest a decomposition of the main optimization problem into sub-problems that are easy to understand and solved via a greedy approach. The theoretical properties of the obtained solution are characterized in terms of approximation bounds while the empirical properties are studied extensively.

We illustrate our approach considering dimensionality reduction with principal component analysis (PCA). Here the goal is not to propose a new fast algorithm to compute the principal directions, plenty of efficient algorithms for this task [15][20][Chapter 10] [22]. Rather, we aim at constructing fast dimensionality reduction operators. While the calculation of the principal components is a one-off computation, a numerically efficient projection operation is critical since it is required multiple times in downstream applications. The problem of deriving fast projections has also been previously studied. Possible approaches include: fast wavelet transforms [49], sparse PCA [48, 52], structured transformations such as circulant matrices [24], Kronecker products [21], Givens rotations [12, 32, 34] or structured random projections [1, 18]. Compared to these works, we propose a new way to factorize any orthogonal matrix, including PCA directions, into simple orthogonal structures that we call extended orthogonal Givens transformations and which naturally lead to optimization problems that have closed-form solutions and are therefore efficiently computed.

We note that our approach provides new perspectives on the structure of the orthogonal group and how to coarsely approximate it, which might have an impact on other open research questions.

The paper is organized as follows: Section 2 describes the basic building blocks and algorithm that we propose, Section 3 gives the theoretical guarantees for our contributions, Section 4 details the application of our method to PCA projections and Section 5 shows the numerical experiments.

2 The proposed algorithm

Given a $d\times d$ orthonormal $\mathbf{U}$, the matrix-vector multiplication $\mathbf{Ux}$ takes $O(d^{2})$ operations. We want to build $\mathbf{\bar{U}}$ such that $\mathbf{U}\approx\mathbf{\bar{U}}$ and $\mathbf{\bar{U}x}$ takes $O(d\log d)$ operations. Parametrizations of orthonormal matrices [2] are known, but to be best of our knowledge, the problem of accurately approximating $\mathbf{U}$ as product of only a few $O(d\log d)$ transformations is open. Given a $d\times p$ diagonal $\mathbf{\Sigma}_{p}$, in the spirit of previous work minimizing Frobenius norm approximations [30, 31], we consider the problem

[TABLE]

where $\mathbf{\bar{U}}$ is $d\times d$ and $\mathbf{\bar{\Sigma}}_{p}$ is a $d\times p$ diagonal matrix. Choosing $p=d$ while $\mathbf{\Sigma}_{p}$ and $\mathbf{\bar{\Sigma}}_{p}$ to be the identity, we simply approximate $\mathbf{U}$. We will also use $\mathbf{U}_{p}$ do denote the first $p\leq d$ columns of $\mathbf{U}$. The above general formulation allows to also consider cases where different directions might have different importance. Then, $\mathcal{F}_{g}$ is a set of orthogonal matrices – defined next – that can be applied fast and allow to efficiently but approximately solve (1).

2.1 The basic building blocks

Classic matrix building blocks that are numerically efficiency include circulant/Toeplitz matrices or Kronecker products. These choices are inefficient as they depend on $O(d)$ free parameters but their matrix-vector product cost is $O(d\log d)$ or even $O(d\sqrt{d})$, i.e., they do not scale linearly with the number of parameters they have. Consider the sparse orthogonal matrices

[TABLE]

where the non-zero part (denoted by $*$ and $\mathbf{\tilde{G}}_{ij}$) on rows and columns $i$ and $j$. The transformation in (2), with the first option in $\mathbf{\tilde{G}}_{ij}$, is a Givens (or Jacobi) rotation. With the second option, we have a very sparse Householder reflector. These transformations were first used by Rusu and Thompson [37] to learn numerically efficient sparsifying dictionaries for sparse coding. The $\mathbf{G}_{ij}$s have the following advantages: i) they are orthogonal; ii) they are sparse and therefore fast to manipulate: matrix-vector multiplications $\mathbf{G}_{ij}\mathbf{x}$ take only $6$ operations; iii) there are two degrees of freedom to learn: $c$ (or $s$) and the binary choice; and iv) allowing both sub-matrices in $\mathbf{\tilde{G}}_{ij}$ enriches the structure and as we will see, leads to an easier (closed-form solutions) optimization problem.

We propose to consider matrices $\mathbf{\bar{U}}\in\mathcal{F}_{g}$ that are products of $g$ transformations from (2), that is

[TABLE]

Matrix-vector multiplication with $\mathbf{\bar{U}}$ takes $6g$ operations – when $g$ is $O(d\log d)$ this is significantly better than $O(d^{2})$, while the coding complexity of each $\mathbf{G}_{ij}$ is approximately $2\log_{2}d+C$: $2\log_{2}d-1$ bits to encode the choice of the two indices, a constant factor $C$ for the pair $(c,s)$ and 1 bit for the choice between the rotation and reflector. The coding complexity of $\mathbf{\bar{U}}$ scales linearly with $g$.

We note that Givens rotations have been used extensively to build numerically efficient transformations [11, 19, 30, 31, 32]. However, $2\times 2$ reflector was not used before. This may be because in linear algebra (e.g. in QR factorization) and in optimization [39] considering also the reflector has no additional benefit: the rotation alone introduces each zero in the QR factorization and the reflector does not have an exponential mapping on the orthogonal manifold, respectively. As we will show, considering both the rotation and the reflector has the advantage of providing a closed-form solution to our problem.

2.2 The proposed greedy algorithm

We propose to solve the optimization problem in (1) with a greedy approach: we keep $\mathbf{\bar{\Sigma}}_{p}$ and all variables fixed except for a single $\mathbf{G}_{i_{k}j_{k}}$ from $\mathbf{\bar{U}}$ and minimize the objective function. When optimizing w.r.t. $\mathbf{G}_{i_{k}j_{k}}$ it is convenient to write

[TABLE]

The next result characterizes the Givens transformation $\mathbf{G}_{i_{k}j_{k}}$ minimizing the above norm. We drop the dependence on $k$ for ease of notation.

Theorem 1** (Locally optimal $\mathbf{G}_{ij}$).**

Let $\mathbf{L}$ and $\mathbf{N}$ be two $d\times p$ matrices. Further, let $\mathbf{Z}=\mathbf{L}\mathbf{N}^{T}$ and $\mathbf{Z}_{\{i,j\}}=\begin{bmatrix}Z_{ii}&Z_{ij}\\ Z_{ji}&Z_{jj}\end{bmatrix}$ then we have

[TABLE]

Let $\mathbf{Z}_{\{i^{\star},j^{\star}\}}\!=\!\mathbf{V}_{1}\mathbf{SV}_{2}^{T}$ be the SVD of $\mathbf{Z}_{\{i^{\star},j^{\star}\}}$ , with

[TABLE]

Then, the extended orthogonal Givens transformation that minimizes $\|\mathbf{L}-\mathbf{G}_{ij}\mathbf{N}\|_{F}^{2}$ is given by $\mathbf{\tilde{G}}_{i^{\star}j^{\star}}^{\star}=\mathbf{V}_{1}\mathbf{V}_{2}^{T}.$

The above theorem derives a locally optimal way to construct an approximation $\mathbf{\bar{U}}$. We iteratively apply the result to find, for each component $k$ in (3), the extended orthogonal Givens transformation that best minimizes the objective function (1). The full procedure is in Algorithm 1 and can be viewed in two different ways: i) a coordinate minimization algorithm; or ii) a hierarchical decomposition where each stage is extremely sparse. The proposed algorithm is guaranteed to converge, in the sense of the objective function (1), to a stationary point. Indeed, no step in the algorithm can increase the objective function, since the sub-problems are minimized exactly: we choose the best indices and then perform the best $2\times 2$ transformation. We note three remarks on the properties of Algorithm 1.

Remark 1** (Complexity of Algorithm 1).**

The computational complexity of the iterative part of Algorithm 1 is $O(dg)$ and the initialization is dominated by the computation of all the scores $C_{ij}$ which takes $O(d^{2})$. Note that, the $C_{ij}$s are computed from scratch only once in the initialization phase. After that, at each step $k$ we need to recompute the $C_{ij}$ (redo the $2\times 2$ singular value decompositions) only for the indices $(i_{k},j_{k})$ currently used (the Givens transformations act on two coordinates at a time). All other scores are update by the same quantity: in (6), the $C_{ij}$ are the same except when $i$ or $j$ belong to the set $(i_{k},j_{k})$. This observation substantially reduces the running time.

Remark 2** (Complexity of applying $\mathbf{\bar{U}\mathbf{\bar{\Sigma}}}_{p}$).**

When $p<d$ the computational complexity of $6g$ operations is an upper bound. Since we keep only $p$ components, we need be careful not to perform operations whose result is thrown away by the mask $\mathbf{\bar{\Sigma}}_{p}$. Consider for example a transformation $\mathbf{G}_{1d}$ applied to a vector of size $d$ projected to a $p<d$ dimensional space. The three operations that take place on the $d^{\text{th}}$ component are unnecessary. Then, after computing $\mathbf{\bar{U}}$, a pass is made through each of the $g$ transforms to decide which of two coordinates the computations are necessary for the final result. As we will show, this further improves the numerical efficiency of our method.

Remark 3** (On the choice of indices).**

Algorithm 1 greedily chooses at each step $k$ the indices according to (6). Other factors might be considered: i) choosing indices based on previous choices so that only a select group of indices are used throughout the algorithm, or ii) make multiple choices at each step in order to speed up the algorithm.

3 Analysis of the proposed algorithm

We consider $p=d$, i.e., $\mathbf{\bar{\Sigma}}_{p}=\mathbf{I}_{d\times d}$ and therefore $\mathbf{Z}=\mathbf{U}$. We model the $\mathbf{U}$ as a random orthonormal matrix with Haar measure [25] updated so that the diagonal is positive. We perform this update because multiplication by a diagonal matrix with $\pm 1$ entries has no computational cost but it brings $\mathbf{U}$ closer to $\mathbf{I}_{d\times d}$. The goal of this section is to establish upper bounds for the distance between $\mathbf{U}$ and $\mathbf{\bar{U}}$, as a function of $d$ and $g$. We first comment on the inherent difficulty of the problem.

Remark 4** (The approximation gap).**

Since the orthogonal group has size $O(d^{2})$, by the pigeonhole principle a random orthogonal matrix as (3) and only $g\ll d^{2}$ degrees of freedom cannot be exactly approximated with less than $O(d^{2})$ operations. For our purposes, think $g$ either $O(d)$ or $O(d\log d)$. Our goal is to show that the fast structures we propose can perform well in practice and have theoretical bounds that guarantee worse case or average accuracy.

Next, we show two approximations bounds depending on the number of Givens transformations (2).

Theorem 2** (A special bound).**

Given a random $d\times d$ orthonormal $\mathbf{U}$, for large $d$, its approximation $\mathbf{\bar{U}}$ from (1) with $g=d/2$ transformations from (2) obeys

[TABLE]

Theorem 3** (A general bound).**

Given a random $d\times d$ orthonormal $\mathbf{U}$, for large $d$, its approximation $\mathbf{\bar{U}}$ from (1) with $g\leq d(d-1)/2$ transformations from (2) is bounded by

[TABLE]

Theorem 2 shows that, on average, the performance might degrade with increasing $d$. As stated in Remark 4, this is not surprising since the orthonormal group is much larger than the structure we are trying to approximate it with. The next result provides a bound for other values of $g$. In Theorem 3, taking $g=c_{1}d\log d$ for some positive constant $c_{1}$ we have that $r\approx c_{1}\log d$. This means that whenever $p\ll d$ we will roughly need $O(d)$ Givens transformations from (2) to improve the $\lfloor r\rfloor$ term. Since the proof of the theorem uses only rotations (and furthermore, in a particular order of indices $(i_{k},j_{k})$) we expect our algorithm to perform much better than the bound indicates as it allows for a richer structure (2) and uses greedy steps that maximally improve the accuracy at each step.

The previous theorems consider the Frobenius norm. In the Jacobi iterative process for diagonalizing a symmetric matrix with Givens rotations [20][Chapter 8.4] the progress of the procedure (convergence) is measured using the off-diagonal “norm” $\text{off}(\mathbf{U})=\sqrt{\sum_{t}^{d}\sum_{q\neq t}^{d}U_{tq}^{2}}$.

Theorem 4** (Convergence in the off-diagonal norm).**

Given a $d\times d$ orthonormal $\mathbf{U}$ and a single Givens transformation $\mathbf{G}_{ij}$, assuming $\text{det}(\mathbf{U}_{\{i,j\}})\geq 0$ we have

[TABLE]

This result shows that, unlike with the Jacobi iterations, monotonic convergence in this quantity is not guaranteed and depends on the relative differences between the diagonal and the off-diagonal entries of $\mathbf{U}_{\{i,j\}}$. Our method convergence monotonically to a stationary point when we measure the progress in the Frobenius norm.

Remark 5** (The effect of a single $\mathbf{G}_{ij}$).**

Given a $d\times d$ orthonormal $\mathbf{U}$ and a Givens transformation $\mathbf{G}_{ij}$ from (2) we have that: i) $\mathbf{U}\mathbf{G}_{ij}^{T}$ is closer to the identity matrix in the sense that $\mathbf{\tilde{G}}_{ij}$ makes a positive contribution to the diagonal elements, i.e., $\text{tr}(\mathbf{U}\mathbf{G}_{ij}^{T})=\text{tr}(\mathbf{U})+C_{ij}$; and ii) $\mathbb{E}[C_{ij}]\approx 0.6956d^{-1/2}$ if $\mathbf{U}$ is random with Haar measure and positive diagonal for large $d$.

The above remark suggests a metric to study the convergence of the proposed method: each Givens transformation adds the score $C_{ij}$ to the diagonal entries of the current approximation (and therefore ensures that $\mathbf{U\bar{U}}^{T}$ converges to the identity – the only diagonal orthonormal matrix). By choosing the maximum $C_{ij}$ we are taking the largest step in this direction.

Remark 6** (Evolution of $C_{i_{k}j_{k}}$ with $k$).**

Given a fixed $0<u<1$, consider the toy construction $\mathbf{U}_{\{i,j\}}=\begin{bmatrix}u&z_{2}\\ z_{1}&u\end{bmatrix}$, i.e., a $2\times 2$ sub-matrix of a $d\times d$ orthonormal matrix where diagonal elements are equal and off-diagonals are two independent truncated standard normal random variables in the interval $[-\sqrt{1-u^{2}},\sqrt{1-u^{2}}]$ (since rows and columns of $\mathbf{U}$ are $\ell_{2}$ normalized). Then, for large $d$, by direct calculation we have that the expected score $\mathbb{E}[C_{ij}(u)]$, i.e., $C_{ij}$ as a function of $u$, obeys

[TABLE]

*i.e., the expected $C_{ij}$ decreases on average quadratically with the increase in the diagonal elements.

The remark is intuitive: as $k$ increases $\mathbf{\bar{U}}$ is more accurate and $\mathbf{U\bar{U}}^{T}$ becomes diagonally dominant, i.e., $\mathbf{U\bar{U}}^{T}\to\mathbf{I}_{d}$ as $k\to O(d^{2})$, and we do expect to reach lower scores $C_{i_{k}j_{k}}$, i.e., few Givens transformations from (2) provide a rough estimation while very good approximations require $k\approx d^{2}$.*

3.1 Other ways to measure the approximation error

Throughout this paper we use the Frobenius norm to measure and study the approximation error we propose. In this section we discuss this choice and explore a few alternatives.

In the linear algebra literature, a natural way to measure approximation error is through the operator norm. This is especially true in the randomized linear algebra field (where matrix concentration inequalities which bound the operator norm play a central role). Moreover, the review manuscript [46] deals explicitly with some potential issues that might arise from using the Frobenius, instead of the operator, norm in matrix approximations: the discussion in Chapter 6 of [46] entitled “Warning: Frobenius–Norm Bounds”. The text highlights situations where the matrix to be approximated is low rank and corrupted by noise and/or scaling issues are present. While that discussion holds true, in our case we deal with a given perfectly conditioned matrix $\mathbf{U}$ and its approximation displays the same scaling ($\mathbf{\bar{U}}$ has normalized columns just as $\mathbf{U}$). Consider the following remark.

Remark 7. [Simplifying the Frobenius norm objective function] Consider for simplicity $p=d$ and $\mathbf{\bar{\Sigma}}_{d}=\mathbf{\Sigma}_{d}$ in (1) and see that $\|(\mathbf{U}-\mathbf{\bar{U}})\mathbf{\Sigma}_{d}\|_{F}^{2}=2\sum_{i=1}^{d}\sigma_{i}(1-\mathbf{u}_{i}^{T}\mathbf{\bar{u}}_{i})=2\sum_{i=1}^{d}\sigma_{i}(1-\cos(\theta_{i}))$, where $\theta_{i}$ is the angle between column vectors $\mathbf{u}_{i}$ and $\mathbf{\bar{u}}_{i}$ (the columns of $\mathbf{U}$ and $\mathbf{\bar{U}}$, respectively), and $\sigma_{i}>0$ are the diagonal elements of $\mathbf{\Sigma}_{d}$. $\square$

Therefore, the proposed optimization objective function minimizes the weighted sum of the cosines of the angles between the original columns and their approximations. As such, low-rank or scaling issues cannot arise. The only potential problem is that different columns might have very different approximation errors (i.e., there could exist $i$ such that $\cos(\theta_{i})\approx 1$ while there could be some $j$ for which $\cos(\theta_{j})\ll 1$). This issue can be mitigated by increasing $g$ in (3) or choosing carefully the indices $(i_{k},j_{k})$ where the proposed transformations operate.

Assume that we consider the operator norm as the approximation error, i.e., we want to minimize $\|\mathbf{U}-\mathbf{\bar{U}}\|_{2}$ in (1). We have the following two results.

Remark 8. [The spectrum of the error matrix] Consider $p=d$ and $\mathbf{\bar{\Sigma}}_{d}=\mathbf{\Sigma}_{d}=\mathbf{I}_{d}$ in (1) then for any orthonormal $\mathbf{U}$ and $\mathbf{\bar{U}}$ the error matrix $\mathbf{U}-\mathbf{\bar{U}}$ is normal and has all its eigenvalues on a circle of radius one centered at $(1,0)$ in the complex plane. As such, we have that $\|\mathbf{U}-\mathbf{\bar{U}}\|_{2}\leq 2$. $\square$

Theorem 5. [A bound on the operator norm of the error matrix] Consider $p=d$ and $\mathbf{\bar{\Sigma}}_{d}=\mathbf{\Sigma}_{d}=\mathbf{I}_{d}$ in (1) and that $\mathbf{u}^{T}_{i}\mathbf{\bar{u}}_{i}\geq 0$ for all $i=1,\dots,d$, then the operator norm of the error matrix obeys $\|\mathbf{U}-\mathbf{\bar{U}}\|_{2}\leq 1-\epsilon_{\text{min}}+\sqrt{(d-1)(1-\epsilon_{\text{min}}^{2})}$ where $\epsilon_{\text{min}}=\underset{i}{\min}\ \mathbf{u}^{T}_{i}\mathbf{\bar{u}}_{i}=\underset{i}{\min}\ \cos(\theta_{i})$. The bound in Remark 8 is met when $\epsilon_{\min}\geq(d-2)/d$. $\blacksquare$

Remark 8 describes the full spectrum of the error matrix. It is interesting to notice that the upper bound $\|\mathbf{U}-\mathbf{\bar{U}}\|_{2}\leq 2$ coincides with the expectation result from [13]: as $d\to\infty$ if both $\mathbf{U}$ and $\mathbf{\bar{U}}$ are chosen uniformly at random with Haar measure then almost surely $\|\mathbf{U}+\mathbf{\bar{U}}\|\to 2$. Theorem 5 describes an upper bound on the operator norm of the error matrix. The bound is tight only for very high values of $\epsilon_{\min}$ and therefore is meant to give a qualitative measure of the approximation. As in Remark 7, the key quantity is $\cos(\theta_{i})$ but now the bound depends on the worst approximation: while the Frobenius norm objective function minimizes the sum of the pairwise distances between the columns of $\mathbf{U}$ and $\mathbf{\bar{U}}$, when using the operator norm the approximation accuracy depends on the largest distance between the same pairwise columns. The assumption that $\mathbf{u}^{T}_{i}\mathbf{\bar{u}}_{i}\geq 0$ is not restrictive at all as we expect $\mathbf{U}^{T}\mathbf{\bar{U}}$ to be diagonally dominant with positive elements on the diagonal (see Remark 5 and the discussion after Remark 6). We note that in this context the Frobenius norm approach could be used to minimize a proxy for the operator norm. As in Remark 7 we have the choice of $\sigma_{i}$ we could update these values with each iteration of the proposed algorithm such that $\sigma_{i}^{\text{(new)}}\leftarrow\frac{\mathbf{u}_{i_{\max}}^{T}\mathbf{\bar{u}}_{i_{\max}}}{\mathbf{u}_{i}^{T}\mathbf{\bar{u}}_{i}}$ where $i_{\max}=\underset{i}{\arg\max}\ \mathbf{u}_{i}^{T}\mathbf{\bar{u}}_{i}$. This choice will encourage the algorithm to improve upon the worst pairwise column approximation (i.e., the highest $\mathbf{u}_{i}^{T}\mathbf{\bar{u}}_{i}$). This approach would be in the spirit of an iteratively reweighted least squares algorithm such as [14].

Finally, the last measure of accuracy we consider is the angle between two subspaces as described by [7]. This value can be non-zero only when $p<d$, a case which is interesting for PCA projections and which we detail in the next section. Following [29], given $\mathbf{U}_{p}$ and $\mathbf{\bar{U}}_{p}$, whose columns span two $d$-dimensional subspaces of size $p$, we compute the angles between the two subspaces as

[TABLE]

where $0\leq\tau_{i}\leq 1$ are the singular values of $\mathbf{U}_{p}^{T}\mathbf{\bar{U}}_{p}$ and where $0\leq\beta_{1}\leq\dots\leq\beta_{p}\leq\pi/2$. We take the principal angle to be the largest angle above, i.e., $\beta=\beta_{p}=\arccos(\tau_{p})$.

Finally, we would like to note that all the approximation errors we have previously discussed measure (in different ways) how well $\mathbf{U}_{p}^{T}\mathbf{\bar{U}}_{p}$ approaches the identity matrix.

4 Application: fast PCA projections

Consider a training set of $d$-dimensional points $\{\mathbf{x}_{i}\}_{i=1}^{N}$ and the $d\times N$ matrix $\mathbf{X}=\begin{bmatrix}\mathbf{x}_{1}&\dots&\mathbf{x}_{N}\end{bmatrix}$. Given $1\leq p<d$, PCA provides the optimal $p$-dimensional projection that minimally distorts, on average, the data points. The projection is given by the eigenvectors of the $p$ largest eigenvalues of $\mathbf{X}\mathbf{X}^{T}$ or, equivalently, the left singular vectors of the $p$ largest singular values of $\mathbf{X}$, that is

[TABLE]

Given the above decompositions we can approximate $\mathbf{X}$ by $\mathbf{\bar{X}}=\mathbf{\bar{U}}\mathbf{\bar{\Sigma}}_{p}\mathbf{V}_{p}^{T}$, i.e., we keep $\mathbf{V}_{p}$ but we modify the principal components and their singular values, such that we minimize the error given by

[TABLE]

where $\mathbf{U}_{p}$ is $d\times p$, the diagonal matrix $\mathbf{\Sigma}_{p}$ is $p\times p$, $\mathbf{\bar{U}}$ is of size $d\times d$, $\mathbf{\bar{\Sigma}}_{p}$ is $d\times p$ and is zero except for its main $p\times p$ diagonal. In this paper, we work with $\mathbf{X}$, as opposed to $\mathbf{XX}^{T}$, to keep the relationship with $\mathbf{U}_{p}$ linear, rather than quadratic. In the context of applying our approach to PCA, we use our decomposition on the principal components $\mathbf{U}_{p}$ which we assume are already calculated together with the associated singular values $\mathbf{\Sigma}_{p}$ which we may use as weights in (1).

Note that in (13), $\mathbf{V}_{p}$, which has size $N$, is not necessary and that the two-step procedure is equivalent to computing the projections $\mathbf{\bar{U}}$ directly from $\mathbf{X}$. Also note that based on (12), we could factor $\mathbf{X}\approx\mathbf{\bar{U}}\mathbf{\Sigma}\mathbf{\bar{V}}^{T}$ where $\mathbf{\bar{U}}$ and $\mathbf{\bar{V}}^{T}$ are approximations in $\mathcal{F}_{O(d\log d)}$ and $\mathcal{F}_{O(N\log N)}$, respectively. The difficulty here is the dependency of $\mathbf{\bar{V}}$ on $N\gg d$ which would require a large running time.

With $\mathbf{\bar{U}}$ fixed, for $\mathbf{\bar{\Sigma}}_{p}$ we have several strategies: i) set it to the identity, i.e., flatten the spectrum; ii) keep it to the original singular values $\mathbf{\Sigma}_{p}$; or iii) continuously update it to minimize the Frobenius norm, i.e., get the new “singular values” that are optimal with the approximation $\mathbf{\bar{U}}$. The first approach favors the accurate reconstruction of all components while the other approach favors mostly the few leading components only (depending on the decay rate of the corresponding singular values).

4.1 Comparison with the symmetric diagonalization by Givens rotations approach

Because of the locally optimal way the Givens transformations are chosen (see Theorem 1), our proposed factorization algorithm is computationally slower than the Jacobi diagonalization process which chooses the Givens rotations on indices $(i,j)$ corresponding to the largest off-diagonal entry of the covariance matrix. Furthermore, the Jacobi decomposition uses each rotation to zero the largest absolute value off-diagonal entry and because of this sub-optimal choice needs $O(d^{2}\log d)$ Givens rotations [10] to complete de diagonalization (more than the $\frac{d(d-1)}{2}$ needed to fully represent the orthonormal group).

In all other aspects, our approach provides advantages over the Jacobi approach: i) we define a clear objective function that we locally optimize exactly; ii) it is known that the Jacobi process converges slowly when the number of rotations is low [28], which is exactly the practically relevant scenario we have, i.e., $g\ll d^{2}$; iii) with the same computational complexity, i.e., $g$ terms in the factorization, our proposed approach is always more accurate since we include as a special case the Givens rotations.

4.2 Comparison with structured matrix factorization

Our approach requires the explicit availability of the orthonormal principal directions $\mathbf{U}_{p}$. Previous methods that factor using only Givens rotations are not applied directly on an orthogonal matrix. These methods rely on receiving as input a symmetric object (e.g., $\mathbf{XX}^{T}$) and then using a variant of Jacobi iterations for matrix diagonalization [31] or multiresolution factorizations [30] to find the good rotations that approximate the orthonormal eigenspace. Applying Givens transformations directly to $\mathbf{X}$ on the left, i.e., $\mathbf{G}_{ij}\mathbf{X}$, cannot lead to the computation of the PCA projections $\mathbf{U}$ but only to the polar decomposition. On the other hand, when applying Givens rotations on both sides of the covariance matrix, i.e., $\mathbf{G}_{ij}\mathbf{XX}^{T}\mathbf{G}_{ij}^{T}$, then the right eigenspace $\mathbf{V}$ cancels out in the product (12) and we are able to directly recover $\mathbf{U}$ (but we need $\mathbf{XX}^{T}$ explicitly). Finally, note that the diagonalization process approximates the full eigenspace $\mathbf{U}$ and cannot separate from the start the $p$ principal components $\mathbf{U}_{p}$ because they are solving the following problem

[TABLE]

This formulation is useful to approximate the whole symmetric matrix $\mathbf{XX}^{T}$ (or $\mathbf{U}$), but not necessarily the $p$ principal eigenspace $\mathbf{U}_{p}$. To get these we would need to complete the diagonalization process, find the $p$ largest entries on the diagonal of $\mathbf{\bar{\Lambda}}$ and then work backward to identify the rotations that contributed diagonalizing those largest elements. This procedure would be prohibitively expensive. Previous work, e.g. [30, 32], deals with approximating $\mathbf{XX}^{T}$ rather than computing PCA. In the same line of work, the one-sided Jacobi algorithm for SVD [16] can also be applied.

5 Experimental results

Given the rather pessimistic guarantees, we tackle problems: how well does Algorithm 1 recover random orthogonal matrices and principal components such that we benefit from the computational gains but do not significantly impact the approximation/classification accuracy. Source code available.111https://github.com/cristian-rusu-research/fast-orthonormal-approximation

5.1 Synthetic experiments

For fixed $d$ we generate random orthonormal matrices from the Haar measure [25]. Figure 2 shows the representation error $(2d)^{-1}\|\mathbf{U}-\mathbf{\bar{U}}\|_{F}^{2}$ for the proposed method. The plot shows that allowing for the Givens transformations $\mathbf{\tilde{G}}_{ij}$ in (2) brings a 17% relative benefit as compared to using only the Givens rotations while, for the same $g$, the computational complexity is the same. The circulant approximation performs worst because it has the lowest number of degrees of freedom, only $d$ (computationally, it is comparable with the Givens and proposed approaches for $g=100$). Lastly, we can observe that the bound is very pessimistic, especially for these values of $g$. In Figure 2 (left) we show for fixed number of transformations $g$ the progress that the proposed algorithm makes with each iteration. It is interesting to observe that the initialization steps (first $g$ steps) significantly decrease the approximation error while the other step make only moderate improvements. This indicates that convergence is slow and might take a large number of iterations with little progress made by the latter steps. Also in Figure 2 (right) we show the number of stages in each transformation. A stage is a set of extended orthogonal Givens transformations that can be applied in parallel, i.e., they do not share any indices $(i_{k},j_{k})$ among them. This is important from an implementation perspective as parallel processing can be exploited to speedup the transformations.

5.2 MNIST digits and fashion

We now turn to a classification problem. We use the MNIST digits and fashion datasets. The points have size $d=400$ (we trimmed the bordering whitespace) and we have $N=6\times 10^{4}$ training and $N_{\text{test}}=10^{4}$ test points. In all cases, we use the k-nearest neighbors (k-NN) algorithm with $k=10$, and we are looking to correctly classify the test points. Before k-NN we apply PCA and our proposed method. Results are shown in Figure 4. We deploy two variants of the proposed method: approximate the principal components as if they had equal importance and approximate the principal components while simultaneously also updating estimates of the singular values.

For comparison, we also show the sparse JL [26]. In this case, the target dimension is $p\in\{15,30\}$ while the random transformation of size $p\times d$ only has three non-zero entries per column. More non-zeros did not have any significant effect on the classification accuracy while increasing (doubling, in this case) the target dimension $p$ increases the classification accuracy by 10%. The results reported in the plots are averages for 100 realizations and the standard deviation is below 1%. Of course, the significant advantage of the JL approach over PCA is that no training is needed. The disadvantage is that if we choose greedily the target projection dimension, i.e., low $p$, the accuracy degrades significantly for JL. Results are identical for PCA when $p$ is $15$ of $30$.

Since we are dealing with image data, we also project using the discrete cosine transform (DCT). For the digits dataset, the performance is poor but, surprisingly, for the fashion dataset, this approach is competitive given the large speedup it reports (we used the sparse fast Fourier transform [23] as we want only the largest $p$ components). We have also performed the projection by fast wavelet transforms with ‘haar’ and several of the Daubechies ‘dbx’ filters but the results were always similar to that of the DCT. For clarity of exposition, we did not add these results to the figures.

Our proposed methods report a clear trade-off between the classification accuracy and the numerical complexity of the projections. If we insist on an accuracy level close (within 1–2%) of the full PCA then the speedup is only about x3. Reasonable accuracy is obtained for a speedup of x4–x5 after which the results degrade quickly. For $p=15$ better performance seems impossible via randomization. In this figure, the speedup is measured in terms of the number of operations (FLOPS).

Finally, in Figure 7 (left) we compare our proposed methods against the sparse PCA on MNIST digits. sPCA performs exceptionally well in terms of the classification accuracy given the computational budget (a similar result is replicated for MNIST fashion). On other datasets where the principal components capture some global features (not local like in our example) we expect this performance to degrade. The training time of sPCA exceeds by 60% the running time of PCA plus that of our method. We used the implementation of Wang et al. [48].

Because of ideas in Remark 1, we perform only the finally useful calculations and further reduce the computational cost on average by one third (these are accounted for already in the numbers in the plots).

5.3 Experiments on other datasets

The 20-newsgroups dataset consists of $18827$ articles from 20 newsgroups (approximately 1000 per class). The data set was tokenized using the rainbow package (www.cs.cmu.edu/ mccallum/bow/rainbow). Each article is represented by a word-count vector for the $d=2\times 10^{4}$ common words in the vocabulary. For this dataset we have $N_{\text{test}}=5648$, $p=200$, and as shown in Figure 7 (right) in this case we outperform sparse PCA.

We also apply our algorithm to several other popular dataset from the literature: PENDIGITS (www.ics.uci.edu/$\sim$mlearn/MLRepository.html) with 10 classes and $d=16$, $N=7494$, $N_{\text{test}}=3498$, $p=4$; ISOLET (archive.ics.uci.edu/ml/datasets/isolet) with 26 classes and $d=617$, $N=6238$, $N_{\text{test}}=1559$, $p=150$; USPS (github.com/darshanbagul/USPS_Digit_Classification) with 10 classes and $d=256$, $N=7291$, $N_{\text{test}}=2007$, $p=12$; UCI (ftp.ics.uci.edu/pub/machine-learning-databases/optdigits) with 10 classes and $d=64$, $N=3823$, $N_{\text{test}}=1797$, $p=6$; EMNIST digits (nist.gov/itl/iad/image-group/emnist-dataset) with 10 classes and $d=784$, $N=24\times 10^{4}$, $N_{\text{test}}=4\times 10^{4}$, $p=15$; MNIST 8m (leon.bottou.org/papers/loosli-canu-bottou-2006) with 10 classes and $d=784$, $N=6.4\times 10^{6}$, $N_{\text{test}}=1.7\times 10^{6}$, $p=15$. All the results are shown in Table 1. Here we provide two measures for the speedup: the FLOPS (number of arithmetic operations, additions and multiplications) and the actual running time (in seconds). For the time speedup, we first computed $\mathbf{\bar{U}}$ and then implemented the matrix-vector multiplication and compared it to the generic matrix-vector multiplication with $\mathbf{U}$, both compiled in C using the gcc compiler and –O3 flag. Once computed by the proposed algorithm, the transformation could also be hardcoded and a further speedup improvement could be achieved. The running time speedups are slightly below the FLOPS speedups due to overhead in modern CPUs (fetching data, register loading etc.). Aside the computational benefits of the proposed transformations we note that while complex instruction set computing machines are closing the speedup gap in terms of running time the price is higher energy consumption and more complex execution pipelines/circuitry.

In Figure 4, for the USPS dataset, we provide insights into the behavior of the proposed algorithm with each iteration and the number of stages in each transformation that we construct. Results are similar to the ones shown for the random orthogonal approximation. We observe again that the initialization step is very efficient in reducing the approximation error (especially when compared to the iterative process that follows) and that the number of stages that have to be applied sequentially can further improve the running time when implementing the proposed transformations.

Finally, one of the most famous applications of PCA is in the field of computer vision for the problem of human face recognition. The eigenfaces [41] approach was used successfully for face recognition and classification tasks. Here, we want to reproduce the famous eigenfaces by using the proposed methods. The original eigenfaces and their approximations (sparse eigenfaces [51]), with different $g$ and therefore different levels of detail, are shown in Figure 7.

5.4 Results on other approximation errors

In Section 3.1 we have explored several other ways to measure the approximation accuracy of the proposed algorithm. We note that in the proposed method we keep the Frobenius norm objective function but we also measure the operator norm, the angle between subspaces distance and the average/minimum/maximum correlations. All results for three datasets are shown in Figure 5, the choices of parameters are the same as in Section 5.3. The first observation is that the approximation errors increase with larger $d$ (as already clear from Remark 4 and Theorem 3), i.e, best results are obtained for UCI ($d=64$) and the worst for MNIST ($d=400$). Second, note that improving the Frobenius norm error we also improve all the other error measures. Thirdly, note the similar behavior between the operator norm and subspace distance errors. It is interesting to observe that the two plots seem to suggest experimentally that $\arccos(\sigma_{\min}(\mathbf{U}_{p}^{T}\mathbf{\bar{U}}_{p}))\approx 1-\sigma_{\max}(\mathbf{I}_{p}-\mathbf{U}_{p}^{T}\mathbf{\bar{U}}_{p})$ which hold over a large array of number of Givens transformations $g$. Lastly, the two right-most plots highlight the connection between the operator norm optimization and the maximization of the lowest coherence between columns and their approximations (as per Theorem 5).

Depending on the application at hand we can choose how to measure the approximation error.

5.5 Details of the implementation

The source code attached to this paper is written for the Matlab environment. Besides the implementation of Algorithm 1, we also provide the code to efficiently perform the matrix-vector multiplication with the proposed $\mathbf{\bar{U}}$ and $\mathbf{\bar{U}}_{p}$, respectively (also taking into account Remark 3). Still, due to the characteristics of Matlab, our implementation is not faster than the dense matrix-vector multiplication “*”. To provide an appropriate comparison, we implement the matrix-vector multiplication with our structures in the compiled lower-level programming language C. We provide comparisons of this implementation against two scenarios: a vanilla matrix-vector multiplication (using the appropriate ordering of the loops to exploit locality) and the Basic Linear Algebra Subprograms (BLAS) Level 2 routines for matrix-vector multiplication (SGEMV). We use the single thread variant of BLAS as for the dimensions $d$ we consider the overhead of the parallel implementation is significant.

Applying the proposed transformations on batch features, i.e., matrix-matrix multiplications, would require careful application of the proposed transformations (2) to fully use the computing architecture and maximize performance.

6 Conclusions

This paper proposes a new matrix factorization algorithm for orthogonal matrices based on a class of structured matrices called extended orthogonal Givens transformations. We show that there is a trade-off between the computational complexity and accuracy of the approximations created by our approach. We apply our method to the approximation of a fixed number of principal components and show that, with a minor decrease in performance, we can reach significant computational benefits.

Future research directions include strengthening the theoretical guarantees since they are way above what we observe experimentally. Further, it would be of interest to improve the complexity of the proposed algorithm either by a parallel implementation or using randomization (e.g. computing a random subset of the $O(d^{2})$ scores). As an immediate application, it would be interesting to apply our decomposition to the recently proposed unitary recurrent neural networks [50].

Acknowledgment

This material is based upon work supported by the Center for Brains, Minds and Machines (CBMM), funded by NSF STC award CCF-1231216, and the Italian Institute of Technology. Part of this work has been carried out at the Machine Learning Genoa (MaLGa) center, Universita di Genova (IT). Lorenzo Rosasco acknowledges the financial support of the European Research Council (grant SLING 819789), the AFOSR projects FA9550-17-1-0390 and BAA-AFRL-AFOSR-2016-0007 (European Office of Aerospace Research and Development). Cristian Rusu acknowledges support by the Romanian Ministry of Education and Research, CNCS-UEFISCDI, project number PN-III-P1-1.1-TE-2019-1843, within PNCDI III.

Supplementary materials

More on Remark 3. We say that the orthonormal group has $O(d^{2})$ degrees of freedom as it was established [2] that any orthonormal matrix $\mathbf{U}$ can be factored into a product as

[TABLE]

where $\mathbf{D}$ is a diagonal matrix with entries only in $\{\pm 1\}$ and the $\mathbf{G}_{ij}(\theta_{ij})$ are Givens rotations, with angles $\theta_{ij}$, i.e., we have $c=\cos\theta_{ij}$ and $s=\sin\theta_{ij}$ in (2). To generate a random orthonormal $\mathbf{U}$ we therefore need to generate random $\mathbf{D}$ (which are $\{\pm 1\}$ with equal probability) and $\frac{d(d-1)}{2}$ random angles $0\leq\theta_{ij}\leq\pi/2$. These angles are mutually independent and it is known that their joint density function is a random variable

[TABLE]

We define the beta random variable $\sqrt{y}=\cos\theta$ (and $\pm\sqrt{1-y}=\sin\theta$) with density

[TABLE]

Because we are interested in the computational complexity of a random orthonormal matrix we focus on the following three special cases: i) do nothing: $\mathbf{\tilde{G}}_{ij}=\begin{bmatrix}\pm 1&0\\ 0&\pm 1\end{bmatrix}$; ii) permute coordinates: $\mathbf{\tilde{G}}_{ij}=\begin{bmatrix}0&\pm 1\\ \pm 1&0\end{bmatrix}$; and iii) $\mathbf{\tilde{G}}_{ij}\in\left\{\frac{1}{\sqrt{2}}\begin{bmatrix}1&1\\ -1&1\end{bmatrix},\frac{1}{\sqrt{2}}\begin{bmatrix}1&-1\\ 1&1\end{bmatrix}\right\}$. The first two cases perform no operations while the last performs only 4 (as compared to 6 for a general rotation). Unfortunately, the joint density in (16) does not seem to have a simlpe closed-form expression. As such we show in Figure 8 numerical results of the probability distribution of $\cos\theta$ over all angles $\theta_{ij}$, which we observe numerically that approaches an exponential distribution $\lambda\exp(-\lambda c)$. Concentration around the three special cases does not occur and therefore a random orthonormal matrix will generally have computational complexity $O(d^{2})$. For example, if we discretized the continuum of $c=\cos\theta$ then the probability that a random $\mathbf{U}$ is an approximate permutation matrix and therefore basically exhibits no numerical complexity is $(1-\exp(-\lambda\epsilon))^{\frac{d(d-1)}{2}}$ for $0<\epsilon\ll 1$, i.e., the probability that all $\frac{d(d-1)}{2}$ rotations have $\cos\theta\leq\epsilon$.$\hfill\square$

Proof of Theorem 1. For simplicity of exposition we will drop the sub-index $k$ herein and therefore (4) develops to the following

[TABLE]

where the Frobenius norms reduces to the $\ell_{2}$ norms of the spectra and we have used the circular permutation property of the trace. It is convenient to denote $\mathbf{Z}=\mathbf{L}\mathbf{N}^{T}$ and the $2\times 2$ matrix $\mathbf{Z}_{\{i,j\}}=\begin{bmatrix}Z_{ii}&Z_{ij}\\ Z_{ji}&Z_{jj}\end{bmatrix}$. Given that $\mathbf{G}_{ij}$ performs operations only on rows $i$ and $j$, the trace is

[TABLE]

To minimize the quantity in (18) we have to maximize (19) which is known as a Procrustes problem [38] whose solution is given by the polar decomposition of $\mathbf{Z}_{\{i,j\}}$ detailed in [20][Chapter 9.4.3]. Therefore, we set the optimal transformation to

[TABLE]

where use the SVD of $\mathbf{Z}_{\{i,j\}}$ ($\mathbf{S}=\text{diag}(s_{1},s_{2})$ are the singular values). With this choice, we have

[TABLE]

We denote the nuclear norm $\|\mathbf{Z}_{\{i,j\}}\|_{*}$, i.e., the sum of the singular values $s_{1}$ and $s_{2}$ and we define

[TABLE]

Intuitively, the results (18), (19) follows after observing that: 1) the $\mathbf{G}_{ij}$ can be viewed as a perturbed identity matrix; 2) if $\mathbf{G}_{ij}$ is exactly $\mathbf{I}_{d\times d}$ then $\|\mathbf{L}-\mathbf{N}\|_{F}^{2}=\|\mathbf{\sigma}_{p}\|_{2}^{2}+\|\mathbf{\bar{\sigma}}_{p}\|_{2}^{2}-2\text{tr}(\mathbf{Z})$ while if $\mathbf{G}_{ij}$ is the optimal orthonormal transformation $\mathbf{Q}$ that minimizes (18) given by the Procrustes solution [38] then we have $\|\mathbf{L}-\mathbf{Q}\mathbf{N}\|_{F}^{2}=\|\mathbf{\sigma}_{p}\|_{2}^{2}+\|\mathbf{\bar{\sigma}}_{p}\|_{2}^{2}-2\|\mathbf{Z}\|_{*}$, where the last term is the nuclear norm of $\mathbf{Z}$; 3) therefore, we actually apply the identity transformation on all coordinates, i.e., the $\text{tr}(\mathbf{Z})$ term, while for the two chosen coordinates we apply the best (in the sense of reducing the error) orthogonal transformation whose contribution is the nuclear norm term $\|\mathbf{Z}_{\{i,j\}}\|_{*}$ and then correct for the trace term that was wrongly added initially in $\text{tr}(\mathbf{Z})$, by subtracting $\text{tr}(\mathbf{Z}_{\{i,j\}})$.

There are $d(d-1)/2$ quantities $C_{ij}$ but they can be computed efficiently by noting that the singular values of $\mathbf{Z}_{\{i,j\}}$ are $s_{1,2}\!=\!\sqrt{\!\frac{1}{2}\!\!\left(\!\|\mathbf{Z}_{\{i,j\}}\|_{F}^{2}\!\pm\!\sqrt{\!\|\mathbf{Z}_{\{i,j\}}\|_{F}^{4}\!-\!4\det(\mathbf{Z}_{\{i,j\}})^{2}}\right)}$. Observe that both singular values are of the form $\sigma_{1,2}=\sqrt{A\pm\sqrt{B}}$ which can be written as $\sqrt{X}\pm\sqrt{Y}$ where $X=\frac{A+\sqrt{A^{2}-B}}{2}$ and $Y=A-X$. Written like this we can see that $\|\mathbf{Z}_{\{i,j\}}\|_{*}=\sigma_{1}+\sigma_{2}=2X=\sqrt{\|\mathbf{Z}_{\{i,j\}}\|_{F}^{2}+2|\det(\mathbf{Z}_{\{i,j\}})|}$. Depending on the sign of the determinant we have either $\|\mathbf{Z}_{\{i,j\}}\|_{*}=\sqrt{(Z_{ii}+Z_{jj})^{2}+(Z_{ij}-Z_{ji})^{2}}$ or $\|\mathbf{Z}_{\{i,j\}}\|_{*}=\sqrt{(Z_{ii}-Z_{jj})^{2}+(Z_{ij}+Z_{ji})^{2}}$, respectively. These give the final formulas for $C_{ij}$ from (5). $\hfill\blacksquare$

Proof of Theorem 2. Assume $d$ is even and partition the set $\{1,\dots,d\}$ into $d/2$ pairs of indices $(i_{k},j_{k})$. We therefore have

[TABLE]

As a side note, to maximization of this partitioned quantity is related to the weighted maximum matching algorithm (of maximum-cardinality matchings) on the graph with $d$ nodes and with edge weights $C_{i_{k}j_{k}}$. With this choice of indices, the objective function becomes

[TABLE]

We use the singular value decomposition of a generic $\mathbf{U}_{\{i,j\}}=\mathbf{V}_{1}\mathbf{SV}_{2}^{T},\mathbf{S}=\text{diag}(\mathbf{s})$, and develop:

[TABLE]

where we have use the circular property of the trace and $\mathbf{\Delta}=\mathbf{V}_{2}^{T}\mathbf{V}_{1}$ where $\Delta_{tt}$ are its diagonal entries which obey $|\Delta_{tt}|\leq\Delta_{\max}$. We define the diagonal coherence as

[TABLE]

With (25), we can state that for the $k^{\text{th}}$ transformation that

[TABLE]

where we have used the fact that $\mathbf{V}_{1}$ and $\mathbf{V}_{2}$ have the structure $\mathbf{\tilde{G}}_{ij}$ in (2) and therefore

[TABLE]

Finally, given an orthonormal $\mathbf{U}$ of size $d\times d,\ d\geq 4$, we use (24) and (27) to bound

[TABLE]

where we have used that $\mathbb{E}\left[\sum_{t=1}^{d}|U_{tt}|\right]=\sqrt{2\pi^{-1}d}$ because the diagonal elements of $\mathbf{U}$ can be viewed as Gaussian random variables with zero mean and standard deviation $d^{-1/2}$ (as the columns of $\mathbf{U}$ are normalized in the $\ell_{2}$ norm) [43] and because the $\ell_{1}$ norm of a standard Gaussian random vector of size $d$ is $\sqrt{2\pi^{-1}}d$.

When $d$ is odd, we extend the matrix with a zero column/row and “1” on the diagonal and the argument follows in the same way. In (24) we could use the expected value calculated in (37) but we reach a worse, lower, constant in (28) for the $-\sqrt{d}$ term and therefore a worse overall bound.$\hfill\blacksquare$

Proof of Theorem 3. Given the orthonormal $\mathbf{U}$, by [20][Theorem 5.2.1], we can construct its QR factorization using a set of Givens rotations [20][Chapter 5.2.5]. After introducing zeros in the first $r$ columns of $\mathbf{U}$, by left multiplication with Givens rotations, we reach its following partial triangularization

[TABLE]

where the diagonal matrix $\mathbf{D}^{\prime\prime}$ of size $r\times r$ has entries $D^{\prime\prime}_{tt}\in\{\pm 1\}$ and $\mathbf{U}^{\prime}$ of size $(d-r)\times(d-r)$ is orthonormal. To introduce the zeros on the $t^{\text{th}}$ column we need $(d-t)$ Givens rotations and therefore to bring $\mathbf{U}$ to the structure in (29) we need $g=\frac{r}{2}(2d-r-1)$ Givens rotations which we have denoted $\mathbf{J}_{i_{k}j_{k}},\ k=1,\dots,g$. We are exploiting the fact that the triangularization of an orthogonal matrix leads to a $\pm 1$ diagonal. Then we might consider a good approximation to $\mathbf{U}$ the product $\mathbf{\bar{U}}=\mathbf{J}_{i_{1}j_{1}}^{T}\dots\mathbf{J}_{i_{g}j_{g}}^{T}\mathbf{D}$. where $\mathbf{D}=\begin{bmatrix}\mathbf{D}^{\prime\prime}&\mathbf{0}_{r\times(d-r)}\\ \mathbf{0}_{(d-r)\times r}&\mathbf{D}^{\prime}\end{bmatrix}$ with $D^{\prime}_{tt}=\text{sgn}(U^{\prime}_{tt})$ and $\mathbf{D}^{\prime\prime}$ is taken from (29). The goal of the diagonal matrix $\mathbf{D}$ is to ensure that the product in $\mathbf{\bar{U}}$ has a nonnegative diagonal. Then, given $g$ transforms we can bound

[TABLE]

If we consider $\mathbf{\tilde{G}}_{ij}$ in (2) instead of the rotations $\mathbf{J}_{i_{k}j_{k}}$ then the quantity on the right becomes an upper bound, since Givens rotations are a special case of $\mathbf{\tilde{G}}_{ij}$ – we can always initialize the $\mathbf{G}_{i_{k}j_{k}}$ of Algorithm 1 with the $\mathbf{J}_{i_{k}j_{k}}$ defined above and the iterative procedure is guaranteed not to worsen the factorization. Therefore, the result follows after using $\mathbb{E}\left[\text{tr}(\mathbf{D}^{\prime}\mathbf{U}^{\prime})\right]=\mathbb{E}\left[\sum_{t=1}^{d-\lfloor r\rfloor}|U^{\prime}_{tt}|\right]=\sqrt{2(d-\lfloor r\rfloor)\pi^{-1}}$.$\hfill\blacksquare$

Proof of Theorem 4. First, we introduce the off-diagonal “norm”, i.e., the square-root of the squared sum of the off-diagonal elements of an orthonormal matrix $\mathbf{U}\in\mathbb{R}^{d\times d}$ as

[TABLE]

Better approximations $\mathbf{\bar{U}}$ of $\mathbf{U}$ lead to lower $\text{off}(\mathbf{U}\mathbf{\bar{U}}^{T})$, as $\mathbf{U}\mathbf{\bar{U}}^{T}$ approaches the identity. If we use this measure, we reach the following result.

We use the fact that $C_{ij}$ is added to the diagonal of $\mathbf{U}$, but only to $U_{ii}$ and $U_{jj}$. The quantity $f(\gamma)=(U_{ii}+\gamma C_{ij})^{2}+(U_{jj}+(1-\gamma)C_{ij})^{2},\gamma\in\mathbb{R},$ is minimized for $C_{ij}\neq 0$ when $\gamma_{0}=\frac{1}{2}+\frac{U_{jj}-U_{ii}}{2C_{ij}}$ (and therefore $f(\gamma_{0})=\frac{1}{2}\|\mathbf{U}_{\{i,j\}}\|_{*}^{2}$) which leads to

[TABLE]

We now use de explicit formulas for the singular values of a $2\times 2$ matrix and the fact that

[TABLE]

and expand this expression to get in (32)

[TABLE]

Therefore, to guarantee $\text{off}(\mathbf{UG}_{ij}^{T})^{2}\leq\text{off}(\mathbf{U})^{2}$ we need $2(U_{ii}^{2}+U_{jj}^{2})\leq\|\mathbf{U}_{\{i,j\}}\|_{*}^{2}$ which is equivalent to

[TABLE]

In this paper we assume that $\mathbf{U}$ is taken randomly from the Haar measure [25] and then modified to have positive diagonal. Therefore, we have $U_{tt}\geq 0$ for all $t$ by construction, otherwise we would just consider the update $U_{tt}\leftarrow\text{sign}(U_{tt})U_{tt}$. Moreover, as the algorithm progresses we continue to have that $\text{det}(\mathbf{U}_{\{i,j\}})\geq 0$ because each $\mathbf{G}_{ij}$ adds a positive amount (the $C_{ij}$ value) to the diagonal elements $U_{ii}$ and $U_{jj}$ thus converging towards $\mathbf{I}_{d\times d}$ in the sense of (36).

In a similar way we can construct a lower bound. Assuming w.l.o.g. that $U_{ii}\geq U_{jj}\geq 0$, the quantity $f(\gamma)=(U_{ii}+\gamma C_{ij})^{2}+(U_{jj}+(1-\gamma)C_{ij})^{2}$ is maximized when $\gamma_{0}=1$ and therefore $f(\gamma_{0})=(U_{ii}+C_{ij})^{2}+U_{jj}^{2}$ which, similarly to (32), leads to

[TABLE]

More on Remark 4. Understanding the properties of these $C_{ij}$ is of crucial importance for our approach. First, notice the effect one generalized Givens transformation has: we have $\|\mathbf{U}-\mathbf{G}_{ij}\|_{F}^{2}=\|\mathbf{U}\mathbf{G}_{ij}^{T}-\mathbf{I}\|_{F}^{2}=2d-2\text{tr}(\mathbf{U}\mathbf{G}_{ij}^{T})$, which together with (18) leads to

[TABLE]

In this sense, the $\mathbf{G}_{ij}$ “pushes” $\mathbf{U}\mathbf{G}_{ij}^{T}$ towards the identity matrix by “contributing” $C_{ij}$ to the diagonal of $\mathbf{U}$, i.e., we are estimating the inverse of $\mathbf{U}$ which in this case is just the transpose.

Notice that $0\leq C_{ij}\leq 4$. The minimum is achieved for symmetric positive semidefinite matrices (because in this case the eigenvalues and singular values are the same and therefore the nuclear norm equals the trace) and the maximum for $\mathbf{U}_{\{i,j\}}=-\mathbf{I}_{2\times 2}$. This immediately leads to a local optimality condition for our approach: there is no $\mathbf{G}_{ij}$ to improve the approximation if all $\mathbf{U}_{\{i,j\}}$ are symmetric positive definite. Assume now we are given a random orthogonal matrix. Because the singular values $s_{1,2}$ depend on the entries of $\mathbf{U}_{\{i,j\}}$ which we model as Gaussian random variables with zero mean and standard deviation $d^{-1/2}$ [43], we have by direct calculation that $\mathbb{E}[s_{1}]\approx 1.7724d^{-1/2}$ and $\mathbb{E}[s_{2}]\approx 0.5190d^{-1/2}$ which leads to $\mathbb{E}[\|\mathbf{U}_{\{i,j\}}\|_{*}]\approx 2.2914d^{-1/2}$. The trace is the sum of two absolute value Gaussian random variables and therefore $\mathbb{E}[\text{tr}(\mathbf{U}_{\{i,j\}})]=2\sqrt{2(\pi d)^{-1}}$ which leads to

[TABLE]

More on Remark 5. To see how the scores $C_{ij}$ depend on the diagonal entries of our toy model, by direct calculation with the truncated standard Gaussian random variables we have that

[TABLE]

We show in Figure 9 the empirical results (mean and standard deviation of $C_{ij}$) on the toy matrix $\mathbf{U}_{\{i,j\}}=\begin{bmatrix}u&z_{2}\\ z_{1}&u\end{bmatrix}$ for $0<u<1$. The empirical mean follows the approximation in Remark 5 (and is tight for $u\leq 0.7$) while we notice that the variance is high for almost the whole interval. In Figure 10 we show the average (left) and maximum (right) costs $C_{ij}$ achieved for another toy model where the diagonal elements are distinct $\mathbf{U}_{\{i,j\}}=\begin{bmatrix}U_{ii}&z_{2}\\ z_{1}&U_{jj}\end{bmatrix}$ for $0\leq U_{ii},U_{jj}\leq 1$.

Finally, notice that when $\mathbf{U}_{\{i,j\}}=\begin{bmatrix}u&U_{ij}\\ U_{ij}&u\end{bmatrix}$ we have $C_{ij}=2(U_{ij}-u)$ if $u\leq U_{ij}$ and zero otherwise, and when $\mathbf{U}_{\{i,j\}}=\begin{bmatrix}u&-U_{ij}\\ U_{ij}&u\end{bmatrix}$ we have $C_{ij}=2\sqrt{u^{2}+U_{ij}^{2}}-2u$ indicating that skew symmetric sub-matrices have higher $C_{ij}$ than symmetric ones, in general.$\hfill\square$

Proof of Remark 7. We use the fact that columns of $\mathbf{U}$ and $\mathbf{\bar{U}}$ have unit norm, and the fact that the Frobenius norm is entrywise:

[TABLE]

Proof of Remark 8. First, note that given any proper norm which is invariant to orthonormal transformations, we have that $\|\mathbf{U}-\mathbf{\bar{U}}\|_{2}=\|\mathbf{I}-\mathbf{U}^{T}\mathbf{\bar{U}}\|_{2}$ and we denote the error matrix $\mathbf{E}=\mathbf{I}-\mathbf{U}^{T}\mathbf{\bar{U}}$. Now, start from the error matrix $\mathbf{E}=\mathbf{I}-\mathbf{U}^{T}\mathbf{\bar{U}}$ and use the fact that $\mathbf{U}^{T}\mathbf{\bar{U}}$ is orthonormal to notice by straightforward calculation that:

[TABLE]

Denote now a complex-valued eigenvalue of $\mathbf{E}$ by $z_{k}=a_{k}+ib_{k}$, then because of the equality above we have that $z_{k}^{*}z_{k}=2\Re(z_{k})$ which in turn can be written as $a_{k}^{2}+b_{k}^{2}=2a_{k}$ or $(a_{k}-1)^{2}+b_{k}^{2}=1$. This is the equation of a circle of radius one centered at $(1,0)$ in the complex plane. Another way to view this result is to observe that $\mathbf{U}^{T}\mathbf{\bar{U}}$ is orthonormal an therefore its spectrum is on the unit circle and then $\mathbf{I}$ shifts the whole spectrum to the right by one unit.

Moreover, we also have that $\mathbf{EE}^{T}=\mathbf{E}+\mathbf{E}^{T}$ and therefore $\mathbf{E}^{T}\mathbf{E}=\mathbf{EE}^{T}$ which means that $\mathbf{E}$ is a normal matrix. As such, the singular values of $\mathbf{E}$ are the absolute values of its eigenvalues and therefore $\|\mathbf{E}\|_{2}\leq 2$.$\hfill\square$

Proof of Theorem 5. The proof is based on the Gershgorin circle theorem (detailed in Chapter 7.2 of [20]) applied to the error matrix $\mathbf{E}=\mathbf{I}-\mathbf{U}^{T}\mathbf{\bar{U}}$. If we denote $\mathbf{Q}=\mathbf{U}^{T}\mathbf{\bar{U}}$, we then have for the $i^{\text{th}}$ eigenvalue of $\mathbf{E}$ that

[TABLE]

The last inequality on the right hand side comes from the $\ell_{1}-\ell_{2}$ inequality $\|\mathbf{x}\|_{1}\leq\sqrt{d-1}\|\mathbf{x}\|_{2}$ applied to the vectors of size $d-1$ (the rows of $\mathbf{Q}$ except for their diagonal elements and whose $\ell_{2}$ norm is $\sqrt{1-Q_{ii}^{2}}$). Returning to (41), for the term on the left hand side we have by the reverse triangle inequality ($|a-b|\geq||a|-|b||$) that

[TABLE]

Finally combining (41) and (42) we have that

[TABLE]

This inequality holds for $0\leq\mathbf{u}_{i}^{T}\mathbf{\bar{u}}_{i}\leq 1$ and $d\geq 2$. The expression on the right-hand side is monotonic decreasing for $\mathbf{u}_{i}^{T}\mathbf{\bar{u}}_{i}\geq 0$. As such, the largest upper bound happens for $i=\underset{i}{\arg\min}\ \mathbf{u}_{i}^{T}\mathbf{\bar{u}}_{i}$. And because, as shown in Remark 8, the error matrix is a normal matrix we get a bound on the operator norm (the singular values of a normal matrix are the absolute values of its eigenvalues). $\hfill\blacksquare$