{'fraction_non_alphanumeric': 0.03986284754618399, 'fraction_numerical': 0.015348012688934503, 'mean_word_length': 4.505714652626172, 'pattern_counts': {'":': 1, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Constructing fast approximate eigenspaces with application to the fast graph Fourier transforms

We investigate numerically efficient approximations of eigenspaces associated to symmetric and general matrices. The eigenspaces are factored into a fixed number of fundamental components that can be efficiently manipulated (we consider extended orthogonal Givens or scaling and shear transformations). The number of these components controls the trade-off between approximation accuracy and the computational complexity of projecting on the eigenspaces. We write minimization problems for the single fundamental components and provide closed-form solutions. Then we propose algorithms that iterative update all these components until convergence. We show results on random matrices and an application on the approximation of graph Fourier transforms for directed and undirected graphs.


Introduction
Matrix decomposition techniques [Stewart, 2000], and specifically eigenvalue decompositions [Golub and van der Vorst, 2000], are widely used in numerical linear algebra, scientific computing, machine learning, quantum computing and other scientific fields.
In general, given no assumptions on the structure of an eigenspace, the eigenvector matrix of a given linear operator of size n × n exhibits no advantageous numerical properties and therefore they require O(n 2 ) operations when performing matrix-vector multiplications. In this paper, we want to perform an approximate eigenvalue decomposition so that we accurately represent the original eigenspace with a new one that also exhibits favorable numerical complexity, for example, it requires only O(n log n) operations when performing matrix-vector multiplication with a generic vector. In such cases, a trade-off between the accuracy of the approximation and its numerical complexity exists.
Several previous works, such as [Lee et al., 2008] and [Kondor et al., 2014], have already introduced these ideas in the machine learning community with considerable success. Such approximations are particularly useful in situations where, once computed, the eigenspace is repeatedly used in matrix-vector calculations in downstream applications.
Eigendecomposition algorithms developed in the matrix computations literature are different for symmetric and unsymmetric matrices. In the symmetric case, the eigenspace is always full (non-defective), real-valued and furthermore, orthonormal [Golub and van Loan, 1996][Chapter 8]. We approximate these eigenspaces by using extended Givens transformations (which are themselves orthonormal and include as a particular case the well-known Givens, sometimes also called Jacobi, rotations [Givens, 1958]). In this case, given the spectrum or an estimation of it, we can provide a locally optimal iterative algorithm similar to Jacobi diagonalization for symmetric matrices [Jacobi, 1846]. The general, unsymmetric, case [Golub and van Loan, 1996][Chapter 7] is much more challenging as the given matrix might not even be diagonalizable and furthermore, even when it is, the factorization has to be done over the complex-valued field in general. As the eigenvector matrix is generally unstructured, in this case, we rely on a given number of scaling and shear transformations to approximate it [Rusu, 2018]. We formulate optimization problems for each of these basic components and show how to locally, optimally solve them with an iterative algorithm and closed-form solutions.
For both proposed algorithms, we show experimental results on the approximation of random unstructured symmetric matrices and then show an application to the construction of fast graph Fourier transforms on synthetic and real-world directed and undirected graphs.

Prior approaches
The literature has always distinguished between eigendecompositions of symmetric and unsymmetric matrices and we will do the same.
In the symmetric case, the diagonalization is done with orthonormal matrices, which are well understood in terms of their decomposition with Givens rotations or Householder reflectors (the QR algorithm, see Chapters 5.1 and 5.2 of [Golub and van Loan, 1996]). The starting point for some approaches in the literature is the Jacobi diagonalization process for symmetric matrices [Jacobi, 1846] which is an iterative procedure that uses Givens rotations to bring the symmetric matrix to a strongly diagonally dominant one. A truncated Jacobi procedure is used in [Le Magoarou et al., 2018] to compute fast graph Fourier transforms for undirected graphs (as undirected implies symmetry in the adjacency and Laplacian matrices). Other methods deal directly with the orthonormal eigenspace. For example, Treelets [Lee et al., 2008] and multiresolution [Kondor et al., 2014] structures use Givens rotations in a structured way to decompose the orthonormal components into hierarchies or multiple different scales, respectively. Another approach is to exploit manifold optimization techniques to find approximate factorizations of orthonormal matrices with few Givens rotations either by greedy coordinate descent [Shalit and Chechik, 2014] or by 1 -style optimization [Frerix and Bruna, 2019]. Recently, an approach that combines rotations and reflections was proposed with an application to fast principal component analysis (PCA) projections [Rusu and Rosasco, 2019]. While this latter work needs to precompute the orthonormal eigenspace, in this paper we show how to perform the same factorization given the dataset.
In the unsymmetric case, we rely on sparse structured components. For example, the incomplete LU [Meijerink and Vorst, 1977], the randomized LU [Shabat et al., 2018] factorizations, the additive low-rank plus multiresolution decomposition [Mudrakarta et al., 2019], and approximate Gaussian elimination [Kyng and Sachdeva, 2016] all rely on structured sparse matrices to construct efficient approximations of a given unstructured matrix.
In this paper, we use structured matrices to construct numerically efficient approximations of eigenspaces. We describe the fundamental building blocks of our factorizations and provide exact optimization problems with closed-form solutions to optimally, locally update these blocks efficiently.
3 Problem setup and formulation 3.1 The symmetric case Given a symmetric matrix S ∈ R n×n the main result that we use is its eigenvalue factorization as where we assume w.l.o.g. that the entries of s, which are the real-valued eigenvalues of S, are in descending algebraic order and U is the real-valued orthonormal eigenspace. Based on (1), we consider the problem where G g a set of orthonormal matrices such that matrix-vector multiplication with any matrix from this set is O(g), instead of the classic O(n 2 ). Let us now consider a particular set G g . Based on all the 2 × 2 orthonormal matrices we have the extended orthonormal Givens transformations [Rusu andThompson, 2017, Rusu andRosasco, 2019], which for simplicity we call a G-transform: where the non-zero entries located at rows/columns i and j, denoted as "*", are the two possible options in (3). The matrices in (3) are basic building blocks of the orthonormal group because every orthonormal matrix U ∈ R n×n can be diagonalized (and the diagonal entries are {±1}) using n(n−1) 2 Givens rotations (the matrix (4) with the first, unsymmetric, component in (3)) by the QR decomposition of U, see Chapter 5.2.5 of [Golub and van Loan, 1996]. Then, in this paper, anyŪ ∈ G g has the following structurē where all matrices G i k j k are G-transforms (4). The number g n 2 is given and fixed. With this structure, matrix-vector multiplicationŪx takes 6g operations while storingŪ takes approximately 2g log 2 n + gC bits, where C is the number of bits required for a double precision floating-point representation. Similar structures to (5) have been previously proposed by [Lee et al., 2008], [Kondor et al., 2014], and [Frerix and Bruna, 2019], but they all consider only Givens rotations, and no reflectors.

The unsymmetric case
Given a general diagonalizable C ∈ R n×n the main result that we use is its eigenvalue factorization as where c contains the complex-valued eigenvalues and T has the complex-valued eigenvectors. Based on (6), we consider the problem where T m a set of general matrices such that matrix-vector multiplication with any matrix from this set or its inverse is O(m), instead of O(n 2 ). Let us now consider a particular set T m . Based on 2 × 2 scaling and shear transformations and, similarly to (4), we define the T-transform: where the non-zero entries are the three possible options in (8). For the shear transformations we necessarily have j > i while for the the scaling transformations we abuse notation and impose i = j, i.e., T ii = T i in (9) is the identity matrix except for the i th diagonal element that is a. The matrices in (9) are building blocks of every diagonalizable C ∈ R n×n because, by Gaussian elimination (see Chapter 3.2.1 of [Golub and van Loan, 1996]), n 2 − n shear transformations (9) diagonalize C and then n scaling transformations (9) exactly represent the resulting diagonal. We choose these three specific matrix as the optimization in (7) takes place over the variableT and its inverse and these matrices have trivial inverses. Then, in this paper, anyT ∈ T m has the following structureT where all matrices T i k j k are T-transforms (9). Their number m n 2 is given and fixed. We assume that the factorization contains m 1 scalings and m 2 shears (m 1 + m 2 = m). With this structure, matrix-vector multiplicationTx takes m 1 +2m 2 operations while storingT takes approximately mC+(m 1 +2m 2 ) log 2 n bits, where C is the number of bits required for a double precision floating-point representation.
There are two important differences when compared to G-transforms. Neither the scaling nor the shears are orthogonal but they are more efficient: two computations and one degree of freedom per transform (as compared to the G-transform where we have 6 operations and one degree of freedom). Therefore, for the same computational cost, we expect T-transforms to provide more accurate approximations. The two types of transforms are connected since any 2 × 2 orthonormal transformation can be written as a product of three shears and scalings by the lifting scheme [Daubechies and Sweldens, 1998].

Proposed factorizations and algorithms
In this section we propose approximate solutions to the optimization problems (1) and (6). Therefore, we distinguish between the symmetric and unsymmetric cases. Furthermore, we analyze separately the initialization and iterative procedures that improve the approximation for each of the two problems. Both (2) and (6) echo the fast circulant matrix-vector multiplication which is possible because every circulant matrix of size n × n has a factorization as F H ΣF, Chapter 4.8 of [Golub and van Loan, 1996], where F is the Fourier matrix and Σ = diag(σ), σ ∈ C n . The idea is to replace the Fourier matrix with a new learned matrix (Ū as (5) orT as (10)), with similar computational properties to the Fourier [Cooley and Tukey, 1965].

Approximation of symmetric matrices
Based on the eigenvalue decomposition, the idea is to approximate U in (1) with a fast approximationŪ (5). The approximation of S would therefore bē S =Ūdiag(s)Ū T . Multiplications with S can be viewed as a sequence of fast multiplications byŪ T , diag(s) and finallyŪ. Matrix-vector multiplication with a diagonal matrix is fast, n operations, so therefore our goal is to constructŪ such that it also has advantageous numerical properties (for example, computations take less than 2n 2 operations). Therefore, based on (5), we propose an approximation as wheres are the estimated eigenvalues of S. These can be the actual eigenvalues of S if they are known, or else they can be randomly initialized or set to the diagonal elements of S -in the latter case, we should ensure entries ofs are distinct for reasons that will be clear later in this section. To find the best approximationS of S we can compute the best approximation to the spectrum by the following lemma.
Lemma 1 Let S and orthogonalŪ be fixed, thens the arg min of the expression S −Ūdiag(s)Ū T 2 F is given bȳ The complexity of computings is O(gn).
Let us now move to approximate the orthogonal eigenspace of S. Given an approximation as (11) where all G-transforms t = k + 1, . . . , g were initialized, we now study the problem of initializing G i k j k such that we minimize where we have defined the symmetric matrix Theorem 1 (Optimal initialization of each G-transform) Let S,s be fixed and all components G itjt = I n for t = k + 1, . . . , g while G itjt = I n for t = 1, . . . , k − 1, then the optimal k th component where we have denoted the quantity and the 2 × 2 symmetric S (k) Theorem 1 provides an efficient way to find the optimal G-transform that minimizes (13): both the indices and the transform values. Starting from k = g we can continue down to k = 1 and initialize in this fashion all g G-transforms in (5). Also, notice that the unified approach we propose (allowing for both the rotation and the reflection) simplifies the results, i.e., we are not looking to optimize an angle of rotation but we get the optimal local solution by an eigenvalue decomposition (the solution to a two-sided 2×2 Procrustes problem). The computational cost of (15) is dominated by the sweep of the indices (O(n 2 ) operations).
Remark 1 (Connection to the Jacobi method) The Jacobi method, see Chapter 8.4 of [Golub and van Loan, 1996], used to diagonalize symmetric matrices, only uses Givens rotations and selects indices (i k , j k ) that correspond to the off-diagonal element of S (k) with the highest magnitude, i.e., we have (15) with A = |S (k) ij |. The Jacobi algorithm is not concerned with the number of Givens rotations used to diagonalize and indeed, in general, more than n 2 rotations are used (see Chapter 8.4.3 of [Golub and van Loan, 1996] for a detailed discussion on how the number of rotations relates to the converge of the method). Furthermore, the Jacobi method uses only Givens rotations (while we have a richer structure for G i k j k given in (3)) and it does not explicitly have an objective function as (2), i.e., there is no reference matrix to be reconstructed in the sense of minimizing a Frobenius norm (the Jacobi method minimizes the squared sum of the off-diagonal entries of the approximation). Also, the Jacobi method does not need an estimate of the eigenvaluess. If we now ignore the eigenvalue information, we can consider jj . Therefore, the calculated score approximates the Jacobi approach when off-diagonal elements are large but the selection criterium for the indices is different as the iterative process makes progress and the working matrix becomes diagonally dominant. Finally, we note that A ij = 0 whenevers ii =s jj , i = j which agrees with theoretical convergence results on the Jacobi method that hold when assuming distinct eigenvalues [Henrici, 1958].
The proposed approach is also significantly different from other previous approaches. As opposed to the approach in [Kondor et al., 2014], by maximizing (15) we find the indices of the optimal 2×2 transform without actually explicitly having to compute it. This saves up computational time (the A are easy to compute: only 10 operations) and also leads to better approximation (as we also consider the reflector simultaneously with the Givens rotation in (3)). The approach in [Frerix and Bruna, 2019] uses again only Givens rotations to perform coordinate descent on the orthonormal manifold using a particular basis for the tangent space such that the exponential map is a Givens rotation. Noting that c s −s c has the same structure as c s s −c 0 1 1 0 , i.e., the reflection can be seen as a coordinate swap followed by a rotation, we can interpret our approach as a simultaneous dual tangent space descent. Unfortunately, efforts to integrate this view with manifold optimization, in general, seem difficult at this stage as many difficulties arise: for example, the logarithmic map of the reflector is complex-valued and therefore it is not clear how to choose a basis for the tangent space corresponding to the reflector. As rotations and reflections are disconnected components the unified approach used in this paper may work only for the objective function and constraints we consider (due to the existence of the Procrustes solutions based on eigendecompositions). Given an approximation as (11) where all G-transforms were initialized, we now study the problem of improving each individual G i k j k iteratively. We want to optimize each G-transform G i k j k sequentially such that we minimize where, due to the invariance of the Frobenius norm to multiplications by G ij and its transpose, we have defined the symmetric matrices Notice that Theorem 1 covered only the case when B (k) is a diagonal matrix. Next we state a general result that holds for the minimization of (17) and any A (k) and B (k) .
Theorem 2 (Optimal update of each G-transform) Let A (k) and B (k) be any symmetric n × n matrices, then the minimizer of the quantity in (17) has its non-trivial values given by where where B The new index f k runs through the two options (rotation and reflector) in (3). The matrices R (f k ) of size 2 × 2, the vectors g (f k ) of length 2 and w of length n 2 , and the matrices M (f k ) and N (f k ) all of size 4 × 4 depend only on the entries in A (k) , B (k) and are given explicitly in the supplementary materials.
Unlike with the initialization procedure, Theorem 2 shows that considering both the rotation and the reflector simultaneously does not lead to a unified optimization problem. Indeed, the index f k runs through both transformations from (3). Still, solving the second problem, for f k = 2, brings an extra computational load that is negligible as it shares most calculations with the first problem, for f k = 1.
The iterative process is computationally expensive as it covers all O(n 2 ) unique pairs of indices (i k , j k ) while the calculation of (20) is itself non-trivial and requires O(n 3 ) operations (substantially more expensive than the initialization (16)). If the running time is an important constraint, we can run the iterative process just as a "polishing step": keep the indices of the G-transforms fixed all the time and update only the values of the transformationsG k .

Approximation of unsymmetric matrices
Similarly to the symmetric case, based now on the eigenvalue decomposition (6), the idea is to approximate T with a numerically efficient approximationT. The approximation of C would therefore beC =Tdiag(c)T −1 . Multiplications withC can be viewed as a sequence of fast multiplications byT −1 , diag(c) and finallyT. Again, matrix-vector multiplication with a diagonal is fast and therefore the computational burden depends on the numerical properties ofT and its inverse. By using scaling and shear transformations (8) in the direct transformationT the numerical properties transfer also to its inverseT −1 since inverses of scalings and shears are themselves scalings and shears, respectively. Therefore, based on (10), we propose an approximation as wherec are the estimated eigenvalues of C, which we constraint to be realvalued. Just like in the symmetric case, there are several ways to setc: randomly, the diagonal values of C or the true eigenvalues, if they are known. To find the best approximationC of C we can compute the best approximation to the spectrum by the following lemma.
Lemma 2 Let C andT be fixed, thenc the arg min of the expression C − Tdiag(c)T −1 2 F is given byc where * is the Khatri-Rao product.
The computational complexity of gettingc is O(n 4 ), we solve a least squares problem of size n 2 × n where the columns ofT −T * T are Kronecker products. Alternatively, we can approximately solve the problem by some iterative method which exploits the Kronecker product structure to efficiently perform matrixvector multiplications.
Let us now move to approximate the eigenspace of T. Given an approximation as (22) where all T-transforms t = 1, . . . , k − 1 were initialized, we now study the problem of initializing T i k j k such that we minimize where we have defined the symmetric matrix Theorem 3 (Optimal initialization of each T-transform) Let C andc be fixed, let all components T itjt = I n for t = 1, . . . , k − 1 while T itjt = I n for t = k + 1, . . . , m, then the optimal k th component T i k j k that minimizes the quantity in (24) is given by where the quantities C , for an index f k that runs through all three options in (8), are rational functions in a k and are given explicitly in the supplementary materials.
Given an approximation as (22) where all T-transforms were initialized, we now study the problem of improving each individual T i k j k iteratively. Therefore, we want to optimize each T-transform T i k j k sequentially such that we minimize where we have defined the matrix Theorem 4 (Optimal update of each T-transform) Let A (k) and B (k) be any n × n matrices, then the optimal k th component T i k j k that minimizes the quantity in (27) is given by where the quantities D , for an index f k that runs through all three options in (8), are rational functions in a k and are given explicitly in the supplementary materials.
Although the initialization and iterative steps look very similar, (26) and (29), they are significantly different from a computational perspective. While Algorithm 1 Approximate eigenspaces factorization.
Input: The symmetric S or general C, the size of the approximation g or m, the update rule for the eigenvalues in {'original', 'update'} and the stopping criterion (default taken to be = 10 −2 ). Output: The linear transformationŪ and spectrums or linear transforma-tionT and spectrumc, the approximate solutions to (2) or (7), respectively. Setup: G i k j k = I n×n , k = 1, . . . , g and compute all scores A ij from to (15) or T i k j k = I n×n , k = 1, . . . , m and compute all scores C ij from to (26); if the update rule for the spectrum is 'update' and the true spectrum is not available thens = diag(S) or,c = diag(C), respectively.
Initialization: for k = g down to 1 initialize each G-transform G i k j k using Theorem 1 in the symmetric case or for k = 1 to g initialize each T-transform T i k j k using Theorem 3 in the general case.

Iterations:
• for k = 1 to g update each G-transform G i k j k (with all others fixed) according to Theorem 2 or for k = 1 to m update each T-transform T i k j k (with all others fixed) according to Theorem 4.
• if rule is 'update' calculate the new spectrums according to Lemma 1 or calculate the spectrumc according to Lemma 2. Similarly to the symmetric case, we now have a locally optimal way of choosing and updating our T-transforms. Also, the update step is again significantly more computationally expensive than the initialization and more are more expensive than the results for the symmetric case. The basic difficulty stems from the fact that we are dealing now with building blocks that are not orthogonal and therefore are not invariant in the Frobenius norm. The exact initialization of T-transforms takes O(n 3 ) while the update takes O(n 4 ) operations. Due to this high computational cost, simplification can be brought to the algorithm: for example, in the update steps we no longer search over every index pair (i k , j k ) but we keep these indices fixed and just calculate the locally optimal coefficient of the transformation a k . We call this step a polishing step and it reduces the computational complexity of updating the T-transforms to O(n 3 ).
We now make two remarks regarding the proposed decomposition for general matrices.
Remark 2 (Using T-transforms for the symmetric case) The ideas of this section can also be applied to the symmetric case. Consider an approxima-tion analogous to (11) based on T-transforms as If the factorization in (11) is based on the eigendecomposition of symmetric matrices, this one is similar to a Cholesky factorization. With this structure we reach minimization problems similar to (27) that have similar solutions to (67) following the same development as in (63), for brevity we omit these formulas. The obvious disadvantage of this approach is that we do not explicitly preserve the eigen-information of the original matrix: we lose control over the spectrum of the approximation and the T-transform do not approximate explicitly the eigenspace. Still, there are some advantages: i) T-transforms (2 operations per degree of freedom) are more efficient than G-transforms (6 operations per degree of freedom) and therefore we expect to get better approximation accuracy for the same numerical complexity or vice-versa, lower computational complexity for the same representation accuracy; and ii) direct and inverse matrix-vector multiplications withS are still efficient, e.g., take at most 4g operations instead of 12g + n withS. An alternative construction is to keep the eigendecomposition but use realvalued approximate eigenvalues and use T-transforms to approximate the orthonormal eigenspacē We lose the orthogonality of the eigenspace but we expect this representation to be more accurate just because, as already discussed, T-transforms have better numerical properties as compared to G-transforms. As every matrix in (3) can be written as four transformations from (8) (two shears and two scalings by [Daubechies and Sweldens, 1998]), we can use the approximation of S from (11) as an initialization to the factorization in (31) with m = 4g.
Remark 3 (An approximate Schur decomposition for S) Notice that in (2), instead of the diagonal containings we can use for example an upper (or lower) triangular matrix which is also sparse, say O(g) off-diagonal elements. The computational complexity of using such an approximation (directly or inversely) would still be O(g) while we expect the approximation accuracy to be better, lower overall S −S 2 F due to the extra degrees in freedom in the new triangular factor. This factorization would be similar to the Schur decomposition C = VJV −1 , where J is upper triangular and V is orthonormal (but all real-valued). This is in contrast with the eigenvalues decomposition of C which is done over the complex values, in general. Figure 1: Approximation accuracy (mean and std) for Laplacians of randomly generated graphs as a function of the number of transformations g going as αn log 2 n. All graphs are randomly generated using the default settings of the GSP box: community graphs (left), Erdos-Renyi random graphs with the probability of a connection between two nodes p = 0.3 (center) and a sensor graphs (right). Given n the number of nodes in the graph, results are shown for n = 512 (dotted), n = 256 (dashed) and n = 128 (solid) and all methods update also the spectrum of the estimation. Top row shows undirected graphs while bottom row show directed graph (created from undirected graphs, direction of the edge between the nodes is decided randomly with probability 0.5) results. Results are averaged over 100 realizations.

The proposed algorithm
Given that we now have optimal ways to initialize (Theorems 1 and 3) and update G and T transforms (Theorems 2 and 4), we are ready the describe the full proposed procedure, in Algorithm 1. For brevity, we describe a single procedure for both the symmetric and general matrices. Whenever we allow for spectrum updates ofs orc we explicitly mention so.
As previously discussed, every step of the algorithm is locally optimal and can only decrease the objective functions we consider. Thus, convergence to a stationary point is guaranteed. Compared to previous methods, there are two characteristics that we would like to highlight at this point: i) in the symmetric case, we use simultaneously both rotations and reflections to construct our approximation; and ii) for each sub-problem we define we can provide a closed-form solution based either on singular and eigenvalue decompositions or least squares. Also, because the proposed method relies on the calculation of scores that span indices i and j it is naturally parallelizable and is amenable to randomized linear algebra techniques.
Proofs of the lemmas and theorems from this section are collected in the supplementary materials.
approx. directly U approx. directly weighted U proposed approx. proposed approx. spectrum update Figure 4: Given random undirected Erdos-Renyi graphs on size n = 1024 with Laplacian L = Udiag(λ)U T we show the average approximation accuracy when we are approximating U directly, i.e., the eigendecomposition is explicitly given, or we are approximating L directly. When the eigendecomposition is available we approximate U and U diag(λ) (a weighted eigenspace) using [Rusu and Rosasco, 2019].

Experimental results
In this section, we describe numerical experimental results with the proposed algorithms and compare them to previous work. Following previous methods [Le Magoarou et al., 2018, Frerix andBruna, 2019], we measure the quality of the approximation using the Frobenius norm objective functions of the optimization problems we consider. Source code for Algorithm 1 is available online 1 . In all the experiments we use 'update' for the spectrum estimation and only polishing in the iterative part of the algorithm (keep indices found in the initialization fixed and optimize only the values of the transforms -monotonicity in the objective function is still preserved). We show an application to the calculation of the fast graph Fourier transforms. Given a graph with n vertices we compute its n×n Laplacian L = D−A where D is the diagonal degree matrix and A is the n × n adjacency matrix, i.e., A ij = 1 if a directed edge exists between nodes i and j and A ij = 0 otherwise. Given the eigenvalue decomposition of the Laplacian, L = UΛU −1 we call U the graph Fourier transformation of the graph. Our goal is now to build approximations of U which enjoy the same numerical complexity as the classic time-domain Fourier transformation, i.e., O(n log n). We distinguish between undirected and directed graphs. With undirected graphs, the adjacency matrix A is symmetric and therefore the Laplacian is symmetric positive semidefinite allowing for an eigendecomposition with an orthonormal basis U as L = UΛU T . In this case, we use G-transforms to approximate the eigenspace. For directed graphs, the Laplacian does not have an orthonormal structure and therefore we use the more general T-transforms in the factorization. In Figure 1 we show approximation results for different types of graphs of different sizes (number of vertices n).
In Figure 2 we compared the proposed method against the previous state of the art for the computation of fast graph Fourier transforms. The results are shown for four graphs: Minnesota graph from [Defferrard et al., 2015] with n = 2642 and 3304 edges, HumanProtein graph from [Rual et al., 2005] with n = 3133 and 6726 edges, Email graph from [Guimerà et al., 2004] with n = 1133 and 5451 edges and the Facebook graph from [Leskovec and Mcauley, 2012] with n = 2888 and 2981 edges. Our proposed method performs best in all these situations.
In this paper, we assume we do not have information about the eigenspace of L, but have access to L itself. Previous work [Rusu andRosasco, 2019, Frerix andBruna, 2019] considered the possibility of performing first the eigendecomposition and then working directly with the eigenspace U. For the same graphs, we show in Figure 3 the evolution of the accuracy of the overall Laplacian L as a function of the number of basic transformations in their factorization. We emphasize again that this metric is different from the one used in Figure 2 on the accuracy of the eigenspace U.
In Figure 4 we analyze the approximation accuracy of estimating the Laplacian L for random undirected Erdos-Renyi graphs with n = 1024 and compare it against the approach in [Rusu and Rosasco, 2019] which needs directly the eigenspace U. We show that our proposed method, especially with the updated spectrum, is the most appropriate to build an accurate approximation of L.
The supplementary materials have further numerical experiments, comparisons, and measurements of the running time of the fast transformations (not just number of operations).

Conclusions
In this paper, we have described algorithms for the efficient approximate computation of eigenspaces which can be efficiently manipulated. We show an application to the computation of the fast graph Fourier transform (both for directed and undirected graphs) and compare against previous approaches from the literature, which we outperform. An open problem, that we cannot address at this time, is the setup of an appropriate theoretical framework where the proposed factorizations and algorithms can be analyzed and their limitations understood.

Supplementary materials
Eigenvalues of 2 × 2 symmetric matrices Given the symmetric 2 × 2 matrix S {i,j} = S ii S ij S ij S jj its two eigenvalues are given by

Proof of Lemma 1
The result follows directly by using the invariance of the Frobenius norm under orthogonal transformations: Then, because the Frobenius norm is entry-wise we have the minimizers = diag(Ū T SŪ).

Proof of Theorem 1
Given that we have initial values for all components G itjt for t = g, . . . , k + 1 (while G itjt = I n for t = 1, . . . , k − 1) and we now want to also initialize the k th component such that we minimize objective function of (2). That expression can be written as where we have defined the cost For convenience, we defined the 2 × 2 symmetric matrices and Z (k) = diag(s)S (k) ∈ R n×n . In the development of (34) we have used the trace definition of the Frobenius norm X 2 F = tr(X T X), the fact that the Frobenius norm is invariant to orthonormal transformation (in particular G-transformation) G i k j k X 2 F = G T i k j k X 2 F = X 2 F and that G i k j k operates only on rows and columns i k and j k . The problem of maximizing the quantity in (35) is known as the two-side orthonormal Procrustes problem [Schonemann, 1968] whose solution in our case is given byG We will use that V k V T k = V T k V k = I. We assume that in the eigenvalue decomposition of S (k) i k j k where the diagonal matrix D k = diag(d k ) contains the eigenvalues in algebraic descending order. The same ordering is also assumed in s. Therefore, by the rearrangement inequality, see Section 10.2, Theorem 368 of [Hardy et al., 1952], and with j k > i k the trace quantity is maximized and it reduces to Therefore, the overall cost (35) reduces to where we have denoted and we noticed that Z (k) {i k i k } in d k are computed by the formulas in (32). Therefore, the minimizer of (34) is (i k , j k ) = arg max

Proof of Theorem 2
With the definitions (18) and (19), the objective function (17) can be expressed to with the cost that is where we have defined A [i k ,j k ] are both matrices of size 2 × (n − 2) composed of only rows i k and j k , but with the columns i k and j k eliminated from both, from A (k) and B (k) respectively. Finally, we have used the operation denotes the entry-wise matrix-matrix product. The quantity in (44) seems difficult to minimize, in the sense that a solution based on an eigenvalue decomposition, such as in (38), does not seems possible. This is because the cost contains a term (the first trace, similar to the one in (35)) whose maximum is given as the solution to the two-side orthonormal Procrustes problem [Schonemann, 1968] by eigenvalue decompositions of A {i k ,j k } but also another term (the second trace) whose maximum is given as the solution to the one-sided orthogonal Procrustes problem by the singular value decomposition of another quantity, B The computational simplification in (35) is possible because one of the matrices is diagonal and therefore the second trace term does not appear in the cost. It is for this reason that we have to analyze separately the initialization and iteratively procedures. To find the minimizer of (17), both in terms of the indices and the values of the orthonormal G i k j k , take the following where (A (k) ) :,t and (B (k) ) :,t are the t th columns of A (k) and B (k) , respectively. We have introduced the matrix are the standard basis vectors for R n and ⊗ is the Kronecker product. For ease, we also define P (1) and P (2) , the two options for P. We have these two variants for P due to the dual structure of (3). Unlike the previous section, where the dual structure enabled the initialization by an eigenvalue decomposition, here we actually need to solve two different, but related, problems. In the development of (46) we have used again the invariance of norms to orthonormal transformations and the fact that the Frobenius norm is element-wise X 2 F = vec(X) 2 F and that vec(ABC) = (C T ⊗ A)vec(B). The rest of the section is dedicated to finding the minimizer of (46) in a numerically efficient manner. Therefore, the minimization of (46) is equivalent to a constrained least squares problem that can be solved efficiently using the singular value decomposition, see Chapter 12.1 of [Golub and van Loan, 1996] and [Gander et al., 1989]. We have to solve the problem twice, for P (1) and P (2) in (47), and keep the best result. The vector w and the matrix P are never explicitly constructed, but we build the products For these, we have the explicit formulas • for R (1) = P (1) T P (1) we have: R (1) • for g (2) = P (2) T w we have: g Here we have also introduces the quantities: that compute to the squared 2 norms of the columns of A (k) and B (k) , respectively. With this setup, the minimizer of (46) with some fixed indices i k and j k is where λ (f k ) is chosen such that the solution x (f k ) has unit 2 norm according to [Gander et al., 1989], the regularization parameter is computed by i.e., the λ i are the generalized real-valued eigenvalues of the 4 × 4 matrices using (50), the overall minimizer of (43) is found by searching for the indices where B The quantities R (f k ) , g (f k ) are already computed and w 2

Proof of Lemma 2
The result follows by using the trace product formula: Then, this is a simple least-squares problem inc where * is the Khatri-Rao product, the Kronecker products between the corresponding columns ofT −T andT.

Proof of Theorem 3
Consider the scenario where we have initial values for all components T itjt for t = 1, . . . , k − 1 (while T itjt = I n for t = k + 1, . . . , m) and we want to also initialize the k th component such that we minimize the quantity with the matrix where we have defined the cost values as and we have used the quantities: The goal is to search for Because all three C (k) i k j k (a k ) are polynomials in a k of degree four and five their minimization therefore reduces to finding the roots of their derivatives and evaluating the values at those points searching for the minimum value. We have used the following explicit inverse formulas: Proof of Theorem 4 The expression (27) can be developed to where we have used twice the fact that (A ⊗ C)(B ⊗ D) = (AB ⊗ CD), twice that vec(ABC) = (C T ⊗A)vec(B) and x ∈ R 2 is a function only of a k . We have :,j ). But this large vector of size n 2 and other large matrices (Kronecker products of size n 2 ) are never explicitly constructed but the objective function is minimized by exploiting the structures (9) and (8). The cost values for the scaling and the two shears respectively are For completeness, we give the explicit formulas for all the quantities in (64) and (65) • for f k = 1 we have: R (1) ; • for f k = 2 we have: R • for f k = 3 we have: R ii k H j k j k , g and we have used the quantities: Similarly to the initialization step, the goal is to search for (f k , i k , j k , a k ) = arg min Again, all three D (k) i k j k (a k ) are polynomials in a k of degree four and five their minimization therefore reduces to finding the roots of their derivatives and evaluating the values at those points searching for the minimum value.

Additional experimental results
In Figure 5, we show the accuracy of the approximation for randomly generated symmetric (both positive definite and indefinite), and general matrices. Details Results are shown for n = 512 (dotted), n = 256 (dashed) and n = 128 (solid) and all methods update also the spectrum of the estimation. Note that the achieved accuracy is better for the positive definite case. In black, for comparison, the results of r-rank approximations: for the symmetric case r = 3αn log 2 n while for the unsymmetric case r = αn log 2 n (for these values we match the numerical complexity of the transformationsŪ andT, we count 2rn operations for matrix-vector multiplication with the r-rank matrix). Results are averaged over 100 realizations.  Figure 6: Average speedup achieved for matrix-vector multiplication between the full eigenspaces versus their approximations using Algorithm 1 for the graphs from Figure 2. We show the FLOP count (count of the number of operations: 6αn log 2 n for the G-transformations and 2αn log 2 n for the Ttransformations) as compared to regular matrix-vector multiplication (2n 2 ) and the actual matrix-vector multiplication runtime as compated to the LAPACK, Level 2 BLAS, implementation (SGEMV). Application of the butterflies is implemented in the C programming language (scripting languages such as Matlab perform very poorly if all g or m basic transformations are applied sequentially). No parallelism is used in any of these experiments. Code runs on a 2.3GHz Quad-Core Intel Core i5 system with 16GB LPDDR3 memory. about the experimental setup are given in the figure caption. The proposed algorithm uses again only a polishing step and not a full transform update, for computational efficiency.
Finally, in Figure 6 we show the speedup achieved by the proposed transformations. Throughout the paper, we define numerically efficient transformations as those that present a low number of additions and multiplications in matrixvector operations, i.e., FLOP (floating-point operations) count. In this figure, we also show the speedup in terms of actual running time and compare it with the FLOP count. We stress that the speedup does not refer to the running time of the proposed Algorithm 1, but the application of the transformations reached by this algorithm. For this figure, the butterfly transformations (G and T transformations) are implemented in the C programming language (a Matlab implementation of the application of these butterflies is hopelessly slow as compared to just matrix-vector multiplication in Matlab, i.e., the "*" operation calls compiled BLAS functions as opposed to parsing and running Matlab scripts).