Title: KALEIDOSCOPE: AN EFFICIENT, LEARNABLE REPRE-SENTATION FOR ALL STRUCTURED LINEAR MAPS

Abstract: Modern neural network architectures use structured linear transformations, such as low-rank matrices, sparse matrices, permutations, and the Fourier transform, to improve inference speed and reduce memory usage compared to general linear maps. However, choosing which of the myriad structured transformations to use (and its associated parameterization) is a laborious task that requires trading off speed, space, and accuracy. We consider a different approach: we introduce a family of matrices called kaleidoscope matrices (K-matrices) that provably capture any structured matrix with near-optimal space (parameter) and time (arithmetic operation) complexity. We empirically validate that K-matrices can be automatically learned within end-to-end pipelines to replace hand-crafted procedures, in order to improve model quality. For example, replacing channel shuffles in ShuffleNet improves classification accuracy on ImageNet by up to 5%. K-matrices can also simplify hand-engineered pipelines-we replace filter bank feature computation in speech data preprocessing with a learnable kaleidoscope layer, resulting in only 0.4% loss in accuracy on the TIMIT speech recognition task. In addition, K-matrices can capture latent structure in models: for a challenging permuted image classification task, a K-matrix based representation of permutations is able to learn the right latent structure and improves accuracy of a downstream convolutional model by over 9%. We provide a practically efficient implementation of our approach, and use K-matrices in a Transformer network to attain 36% faster end-to-end inference speed on a language translation task.

1. INTRODUCTION:

Structured linear maps are fundamental and ubiquitous in modern machine learning. Their efficiency in speed (fast algorithms) and space (few parameters) can reduce computation and memory usage. The class of structured linear maps includes fixed specialized transforms such as the discrete Fourier transform (DFT) and Hadamard transform used in signal processing (Cooley et al., 1969) , convolutions for image, language, and speech modeling (Gu et al., 2018) , and low-rank and sparse matrices for efficient storage and inference on edge devices (Yu et al., 2017) . Forms of structure such as sparsity have been at the forefront of recent advances in ML (Frankle & Carbin, 2019) , and are critical for on-device and energy-efficient models, two application areas of tremendous recent interest (Tsidulko, 2019; Schwartz et al., 2019) .

There are a plethora of classes of structured linear maps, each with a significantly different representation, algorithm, and implementation. They have different tradeoffs in terms of inference speed, training speed, and accuracy, and the conventional wisdom is that no one class works uniformly well across all applications. As a result, ML practitioners currently hand-pick specific classes of structured linear maps for each of their applications. This is a difficult and labor-intensive task.

Ideally, these problems should be addressed with a universal representation for structured linear maps: (i) Such a parameterization should be expressive enough to capture important classes of structure, with a nearly tight parameter count and runtime: the space required to represent the linear map should be close to optimal, and the resulting algorithm for matrix vector multiplication should be close to the fastest possible algorithm. (ii) The parameterization should be differentiable in order to be learned as a component of end-to-end ML pipelines, enabling it to easily be used as a drop-in replacement for manually engineered structured components. (iii) The parameterization should admit practically efficient algorithms for training and inference, in terms of both speed and memory.

Currently, no class of structured linear maps satisfies all of these criteria. Most existing classes of structured matrices-such as the class of low-rank matrices-fail to tightly capture other important types of structure. For example, the DFT has an efficient structured representation of size O(n log n), yet cannot be well-approximated by a low-rank transform of size n 2 . Another important type of structure is sparsity; lots of exciting recent work has focused on the design of sparse neural networks. For instance, sparse networks of comparable quality to their dense counterparts-yet an order of magnitude fewer parameters-may be created via pruning (Han et al., 2016) or by identifying "winning lottery tickets" (Frankle & Carbin, 2019) . In parallel, recent theoretical results by De Sa et al. (2018) show that sparsity and the notion of structure in linear maps are fundamentally linked: any given matrix can be factored into a product of sparse matrices with total parameter count equal to the efficiency (i.e. minimum arithmetic circuit complexity) of the matrix. In other words, the representation of linear maps as products of sparse matrices tightly captures all forms of structure. Unfortunately, it is difficult to actually learn these sparse factorizations, because it requires finding the sparsity patterns of the factors-a discrete, nondifferentiable search problem. Thus, current methods for training sparse neural networks are either expensive (Frankle & Carbin, 2019) or rely on highly hand-tuned heuristics for evolving the sparsity patterns throughout training (Dettmers & Zettlemoyer, 2019) .

By contrast, we propose a representation of linear maps as products of sparse matrices with specific predefined sparsity patterns (Section 2), and show that it does satisfy our desiderata: it retains the expressiveness of unstructured sparsity, while being differentiably learnable and efficient like other structured representations. Concretely, our representation is based on products of a particular building block known as a butterfly matrix (Parker, 1995; Dao et al., 2019) ; we term such products kaleidoscope matrices (K-matrices for short).1 (i) Our main theoretical contribution (Section 2.3) concerns the expressiveness of this representation: we show that any structured linear map (i.e. one that can be applied using s n 2 arithmetic operations) can be represented as a K-matrix, with a nearly tight number of parameters and algorithmic complexity (both on the order of s up to logarithmic factors). (ii) The kaleidoscope representation is fully differentiable; thus, all the parameters of a K-matrix can be learned using standard optimization algorithms such as SGD. (iii) Because of their simple, regular structure, K-matrices are practical and easy to use. We provide memory-and runtime-efficient implementations of K-matrix multiplication on CPU and GPU for training and inference, with a simple PyTorch interface.

We empirically validate that, due to their expressiveness, learnability, and efficiency, we can use K-matrices as a drop-in replacement for linear components in deep learning models. In Section 3.1, we use K-matrices to replace hand-crafted structure in two different settings. We simplify the six steps of filter bank computation in speech preprocessing into a single learnable K-matrix step, with only an 0.4% accuracy drop on the TIMIT speech recognition task. We use K-matrices to replace channel shuffles in ShuffleNet, improving ImageNet classification accuracy by up to 5%. In Section 3.2, we show that K-matrices can successfully recover latent structure; a K-matrix is used to learn latent permutations in a permuted image dataset (Permuted CIFAR), resulting in 9 points higher accuracy in a downstream CNN model. In Section 3.3, we show that our efficient K-matrix multiplication implementation can be applied to speed up real-world tasks: we replace linear layers with K-matrices in a DynamicConv-Transformer network to attain 36% faster end-to-end inference speed with a 1.0 drop in BLEU score on the IWSLT14 German→English translation task.

2. A NEARLY-TIGHT PARAMETERIZATION OF ALL STRUCTURED MATRICES:

We first present some background on the characterization of all structured matrices (i.e. those with subquadratic multiplication algorithms) as products of sparse factors, along with the definition of butterfly matrices. We then propose a differentiable family of kaleidoscope matrices, composed of products of butterfly matrices, and prove their expressivity: all structured matrices can be represented in this form, with almost optimal parameter count and runtime.

2.1. BACKGROUND: SPARSE FACTORIZATION, BUTTERFLY MATRICES:

Sparse factorization One method of constructing matrices with theoretically fast matrix-vector multiplication algorithms is as a product of sparse matrices, so that multiplication by an arbitrary vector has cost proportional to the total number of nonzeros (NNZ) of the matrices in the product. Surprisingly, the converse is also true. De Sa et al. (2018) introduce the concept of sparse product width (SPW), which roughly corresponds to the total NNZ in a factorization of a matrix, and show that it is an asymptotically optimal descriptor of the algorithmic complexity of matrix-vector multiplication (Bürgisser et al., 2013) . We use a similar argument in the proof of our main theorem (Section 2.3). However, attempting to learn such a factorization of a given matrix is difficult, as the sparsity constraint is not continuous. Moreover, because of the possibly irregular sparsity patterns, it is difficult to realize the theoretical speedups in practice (Gray et al., 2017; Gahvari et al., 2007) .

Butterfly matrices Butterfly matrices, encoding the recursive divide-and-conquer structure of the fast Fourier transform (FFT) algorithm, have long been used in numerical linear algebra (Parker, 1995; Li et al., 2015) and machine learning (Mathieu & LeCun, 2014; Jing et al., 2017; Munkhoeva et al., 2018; Dao et al., 2019; Choromanski et al., 2019) . Here we define butterfly matrices, which we use as a building block for our hierarchy of kaleidoscope matrices. Definition 2.1. A butterfly factor of size k ≥ 2 (denoted as B k ) is a matrix of the form

B k = D 1 D 2 D 3 D 4

where each D i is a k 2 × k 2 diagonal matrix. We restrict k to be a power of 2.

Definition 2.2. A butterfly factor matrix of size n with block size k (denoted as

B (n)

k ) is a block diagonal matrix of n k (possibly different) butterfly factors of size k:

B (n) k = diag [B k ] 1 , [B k ] 2 , . . . , [B k ] n k

Definition 2.3. A butterfly matrix of size n (denoted as B (n) ) is a matrix that can be expressed as a product of butterfly factor matrices:

B (n) = B (n) n B (n) n 2 . . . B (n)

2 . Equivalently, we may define B (n) recursively as a matrix that can be expressed in the following form:

B (n) = B (n) n [B ( n 2 ) ] 1 0 0 [B ( n 2 ) ] 2 (Note that [B ( n 2 ) ] 1 and [B ( n 2 )

] 2 may be different.)

2.2. THE KALEIDOSCOPE HIERARCHY:

Using the building block of butterfly matrices, we formally define the kaleidoscope (BB * ) hierarchy and prove its expressiveness. This class of matrices serves as a fully differentiable alternative to products of sparse matrices (Section 2.1), with similar expressivity. In Appendix J, we show where various common structured matrix classes are located within this hierarchy.

The building block for this hierarchy is the product of a butterfly matrix and the (conjugate) transpose of another butterfly matrix (which is simply a product of butterfly factors taken in the opposite order).

Figure 1 visualizes the sparsity patterns of the butterfly factors in BB * , where the red and blue dots represent the allowed locations of nonzero entries. Definition 2.4 (Kaleidoscope hierarchy, kaleidoscope matrices).

• Define B as the set of all matrices that can be expressed in the form B (n) (for some n).

• Define BB * as the set of matrices M of the form The kaleidoscope hierarchy, or (BB * ) hierarchy, refers to the families of matrices (BB * ) 1 e ⊆ (BB * ) 2 e ⊆ . . . , for a fixed expansion factor e. Each butterfly matrix can represent the identity matrix, so (BB * ) w e ⊆ (BB * ) w+1 e . We show that the inclusion is proper in Appendix E. This hierarchy generalizes the BP hierarchy proposed by Dao et al. (2019) , as shown in Appendix J.

M = M 1 M * 2 for some M 1 , M 2 ∈ B.

Efficiency in space and speed Each matrix in (BB * ) w e is a product of 2w total butterfly matrices and transposes of butterfly matrices, each of which is in turn a product of log(ne) factors with 2ne nonzeros (NNZ) each. Therefore, each matrix in (BB * ) w e has 4wne log(ne) parameters and a matrix-vector multiplication algorithm of complexity O(wne log ne) (by multiplying the vector with each sparse factor sequentially). We prove this more formally in Appendix E. For the applications in Section 3, w and e are small constants (up to 2), so those K-matrices have O(n log n) parameters and runtime.

2.3. ALL LOW-DEPTH STRUCTURED MATRICES ARE IN THE KALEIDOSCOPE HIERARCHY:

We now present our main theoretical result: the fact that general linear transformations, expressed as low-depth linear arithmetic circuits, are captured in the BB * hierarchy with low width. Arithmetic circuits are commonly used to formalize algebraic algorithmic complexity (Bürgisser et al., 2013) ; we include a primer on this in Appendix M. The quantities of interest are the total number of gates in the circuit, representing the total number of steps required to perform the algorithm for a serial processor, and the depth, representing the minimum number of steps required for a parallel processor. Theorem 1. Let M be an n × n matrix such that multiplication of M times an arbitrary vector v can be represented as a linear arithmetic circuit with s total gates and depth d.

Then, M ∈ (BB * ) O(d) O( s n ) .

The representation of such a matrix M in the BB * hierarchy has O(ds log s) parameters and yields a O(ds log s) multiplication algorithm, compared to the O(s) parameters and runtime of the circuit representation. To the best of our knowledge, the most general classes of efficient matrices that have been studied (De Sa et al., 2018) have depth d on the order of log n or poly log n. In these cases, the representation with K-matrices matches the best known bounds up to polylogarithmic factors.

The crux of the proof of Theorem 1 (shown in Appendix F) is the construction of an almost tight representation of any sparse matrix as a K-matrix (i.e. a product of butterfly matrices): specifically, we show that any n × n sparse matrix with s nonzeros is in

(BB * ) O( s n ) O(1)

(Theorem 3, Appendix I). We then leverage the expressivity result of products of sparse matrices to represent all arithmetic circuits (similar to the sparse product width result of De Sa et al. (2018) referenced in Section 2.1) to complete the proof of Theorem 1. This intermediate result is also a novel characterization of sparse matrices. For a matrix with s NNZ, the kaleidoscope representation has O(s log n) parameters and runtime, instead of the optimal O(s) parameters and runtime; so, we trade off an extra logarithmic factor in space and time for full differentiability (thanks to the fixed sparsity patterns in the representation). The intuition behind the result is as follows: a sparse matrix with s NNZ can be written as a sum of s/n matrices each with at most n NNZ. Any n × n matrix with at most n NNZ, up to permuting the rows and columns, is a product of two butterfly matrices (Lemma I.1). Sorting networks (Knuth, 1997) imply that permutation matrices are in (BB * ) O(log n) , but we tighten the result to show that they are in fact in BB * (Theorem 2, Appendix G). We thus obtain a kaleidoscope representation for each summand matrix with O(n log n) parameters. By the addition closure property of the BB * hierarchy (Lemma H.5), each sparse matrix with s NNZ then has a kaleidoscope representation with O(s log n) parameters.

Tight representation for structured linear maps common in ML Even though Theorem 1 suggests that the kaleidoscope representation can be loose by logarithmic factors, many structured linear maps common in ML can be represented in this hierarchy with an optimal number of parameters and runtime compared to the best known parameterizations, up to constant factors. Appendix J includes several examples such as discrete transforms (the DFT, discrete cosine transform (DCT), discrete sine transform (DST), and Hadamard transform), convolution (i.e. circulant matrices), Toeplitz matrices (Gray, 2006) , structured matrices for kernel approximation ((HD) 3 (Yu et al., 2016) ) and compact neural network design (Fastfood (Le et al., 2013) , ACDC (Moczulski et al., 2016) ). There have been other large classes of structured matrices proposed in the machine learning literature, such as Toeplitz-like (Sindhwani et al., 2015) and low displacement rank (LDR) (Thomas et al., 2018) , but they are not known to be able to capture these common structures as tightly as K-matrices can. More detailed discussions are in Appendix A.

2.4. EXTENSIONS:

ReLU networks with low-depth structured weight matrices In Appendix L, we prove that finding an efficient circuit for a ReLU network can be reduced to finding efficient circuits for each of its weight matrices, with at most a constant factor greater size and run-time (i.e. number of gates). We also show that ReLU networks with kaleidoscope weight matrices have near-linear VC dimension in the number of parameters, matching the bound for networks with unconstrained weight matrices (Bartlett et al., 1999; Harvey et al., 2017) and LDR (Thomas et al., 2018) . This yields a corresponding sample complexity bound.

Orthogonal kaleidoscope hierarchy Orthogonal butterfly matrices are one commonly used variant due to their improved stability (Parker, 1995) , where each butterfly factor is constrained to be orthogonal:

C S -S C with C, S being diagonal and C 2 + S 2 = I. Similar to the BB * hierarchy, in Appendix K, we define the OBB hierarchy consisting of products of orthogonal butterfly matrices and diagonal matrices, and show that this hierarchy has the same expressiveness as the BB * hierarchy.

3. EMPIRICAL EVALUATION:

We validate three claims that suggest that kaleidoscopes are a promising technique to learn different types of structure in modern architectures.

1. Section 3.1: for applications in speech and lightweight computer vision relying on highly hand-crafted structured transformations, we show that we can recover-and even improve-the quality of such architectures by simply replacing existing hand-structured components with K-matrices, with only a small overhead in memory and computation. 2. In Section 3.2, for a challenging task with latent structure (Permuted CIFAR-10), a K-matrixbased relaxation of permutations is able to learn the right latent permutation, yielding 9 points better accuracy in a downstream CNN compared to standard RNN and CNN baselines used on such permuted image classification tasks. 3. In Section 3.3, we show that, although not yet highly optimized, our current implementation of K-matrices can improve the inference throughput of DynamicConv Transformer, a state-ofthe-art fast machine translation model, by 36%, with only a relatively small drop in translation quality.

In all of the above applications, as K-matrices are fully differentiable, we simply train them jointly with the rest of the model using standard learning algorithms (such as SGD). Full details for all of the experiments (precise architectures, hyperparameters, etc.) are in Appendix B2 .

3.1. REPLACING HAND-CRAFTED STRUCTURES:

We validate that kaleidoscope matrices can recover or improve on the performance of hand-crafted structure in ML models. For example, a single learnable kaleidoscope layer can be used to replace the hand-engineered filter bank speech preprocessing pipeline with only 0.4% loss in accuracy on the TIMIT speech recognition task (Section 3.1.1). Replacing channel shuffles in ShuffleNet with learnable K-matrices improves classification accuracy on ImageNet by up to 5.0% (Section 3.1.2).

3.1.1. SPEECH PREPROCESSING:

Figure 2 : Comparison of the standard MFSC featurization pipeline with our "kaleidoscope" pipeline.

We show that K-matrices can remove the need for hand-tuning by significantly simplifying speech recognition data preprocessing pipelines. In particular, we can entirely replace the complex handcrafted MFSC featurization commonly used in speech recognition tasks with a fully learnable kaleidoscope layer, with only 0.4% drop in accuracy on the TIMIT speech recognition benchmark.

Results are presented in Table 1 . Our approach is competitive with the accuracy of standard models that use hand-crafted features, and significantly outperforms current approaches for learning from raw audio input. Modern speech recognition models currently rely on carefully hand-crafted features extracted from the audio, which are then fed into an acoustic model. By contrast, learning directly from the raw audio-i.e. end-to-end learning from the audio waveform without any manual featurization-obviates the need for this complicated and often expensive preprocessing step. There have been recent attempts to learn directly from raw audio, such as SincNet (Ravanelli & Bengio, 2018) ; however, they often rely on specialized architectures designed by domain experts. Instead, we use a standard RNN speech recognition architecture, but use a learnable kaleidoscope layer to replace the featurization steps.

The baseline architecture takes as input filter bank (MFSC) features, which are a popular standard featurization for speech recognition (Paliwal, 1999) and involve several steps hand-crafted specifically for this domain. These features are extracted from the raw audio waveform, and fed as the input into a Bi-LSTM model. We significantly simplify this pipeline by replacing the featurization step with a trainable kaleidoscope layer that is trained end-to-end together with the Bi-LSTM. The original pipeline and our modified kaleidoscope version are depicted in Figure 2 .

The computation of MFSC features involves a series of painstakingly hand-designed steps (further described in Appendix B.1), each involving their own hyperparameters: (i) the waveform is framed (split into chunks), (ii) the waveform is dithered (noise is added), (iii) pre-emphasis is applied, (iv) the Hamming window is applied, (v) the FFT is applied and the power spectrum is computed, (vi) the result is mapped to the mel scale (which involves applying a particular linear transformation and then taking the logarithm of the result), (vii) cepstral mean and variance normalization is applied.

We replace the last six steps (ii-vii) of this featurization process with a learnable kaleidoscope layer; specifically, after windowing, we multiply the input by a K-matrix, and then compute the logarithm of the power spectrum; the output is fed into the Bi-LSTM model.

3.1.2. REPLACING CNN CHANNEL SHUFFLE:

We evaluate how K-matrices can improve the quality of hand-crafted, lightweight architectures for computer vision tasks, without the need for hand-tuning. We select ShuffleNet (Zhang et al., 2018) , which is a state-of-the-art lightweight CNN architecture that uses a manually designed "channel shuffle" permutation matrix to improve performance. By replacing this fixed permutation with a learnable K-matrix, we achieve up to 5% further improvement in classification accuracy, without hand-tuned components and with a modest space penalty of up to 10%. Results are given in Table 2 .

Table 2 : Top-1 classification accuracy of ShuffleNet on ImageNet validation set (parameter counts in parentheses). We compare our approach (col. 3) with our reimplementation of 'vanilla' ShuffleNet (col. 1) and a recent approach based on the Hadamard transform (col. 2). 4 We report results for different network width multipliers (# channels). The last column shows the differences in accuracy and parameter count between our approach and vanilla ShuffleNet; using a learnable K-matrix in place of each fixed permutation (shuffle) or Hadamard matrix improves accuracy by up to 5%. Grouped convolution (Krizhevsky et al., 2012) is often used to reduce parameter count and speed up inference compared to standard convolution, but, by default, channels in different groups cannot exchange information. To remedy this, ShuffleNet uses a permutation matrix to shuffle the channels after each grouped convolution. Zhao et al. (2019) propose to instead use the Hadamard transform before and after each grouped convolution to mix the channels. In place of these hand-engineered solutions, we use a K-matrix before and after each grouped convolution, and learn these end-to-end together with the rest of the network. As shown in Table 2 , across a range of sizes, replacing the channel shuffles with K-matrices results in improved performance at comparable parameter counts.

3.2. LEARNING A LATENT PERMUTATION:

We show that K-matrices can be used in a challenging task for which existing classes of structured linear maps have not been found suitable. We investigate the problem of image classification on a permuted image dataset (Permuted CIFAR-10). This problem is challenging due to the discrete nature of learning the latent permutation of the dataset; we present a differentiable relaxation for this using a K-matrix as a key component. Results are presented in Table 3 ; compared to methods that do Published as a conference paper at ICLR 2020 not have a permutation learning step, our approach gets 9 points higher accuracy (84.4% to 93.6%), coming within 2 points of the accuracy on the un-permuted dataset (94.9%). In this task, we use a permuted image classification dataset (Permuted CIFAR-10), wherein a fixed global permutation is applied to the pixels of every image in the original input set. Typically, only fully-connected (FC) and recurrent models are applied to such datasets (Le et al., 2015) , because the permutation destroys locality in the image, presenting a difficulty for CNNs. However, CNNs are much better-suited for standard image tasks. We thus expect that learning the permutation and then applying a standard CNN should outperform these baselines. As mentioned in Section 2, the kaleidoscope hierarchy provides a nearly tight parameterization of permutations; this makes them a natural fit for the permutation learning step.

Experimentally, we use a K-matrix to represent a distribution over permutations, which converges to a single permutation at the end of training. The correct latent structure is learned by applying samples from this distribution to the permuted training images, and minimizing an auxiliary smoothness-based loss that encourages the reconstructed images to be more "natural" (i.e. vary smoothly pixel-to-pixel).

The learned permutation is evaluated by training a ResNet18 with the K-matrix permutation layer inserted at the beginning. Full details of our approach are provided in Appendix B.3.

In Table 3 , we compare our approach to a ResNet18 without this extra K-matrix layer, a ResNet18 with an extra dense matrix at the beginning instead of a K-matrix, and other baselines. As generic representations such as unstructured matrices do not have the requisite properties to fit in the pipeline, these baselines fail to effectively learn the latent permutation. We emphasize that a K-matrix provides this ability to recover latent structure despite not being specialized for permutations. 

3.3. SPEEDING UP INFERENCE:

We evaluate the inference speed benefit of using K-matrices on a real language translation model. We choose the state-of-the-art DynamicConv Transformer translation model (Wu et al., 2019) , which offers 20% inference speedup over the standard Transformer model, and replace dense matrices in the decoder's linear layers with K-matrices, which leads to a further 36% inference speedup (Table 4 ).

As outlined in Section 2.3, K-matrices admit a simple and fast O(n log n) matrix-vector multiplication algorithm. We provide fast implementations of this algorithm in C++ and CUDA, with an interface to PyTorch (Paszke et al., 2017) , and use this implementation in our experiments. We use K-matrices to replace all the linear layers in the decoder of DynamicConv (since 90% of inference time is spent in the decoder). As shown in Table 4 , on the IWSLT-14 German-English translation task, this yields a 25% smaller model with 36% faster inference time on CPU, at the cost of 1.0 drop in BLEU score.5 (Our model also nearly matches the state-of-the-art BLEU performance of 2 years ago obtained by the Transformer model (Vaswani et al., 2017) , despite being over 60% faster for inference than the Transformer.) The majority (55%) of inference time is spent in matrix-vector multiplication; our implementation of K-matrix-vector multiplication is about 2 times faster than the optimized implementation of dense matrix-vector multiplication in the Intel MKL library. Direct comparisons of K-matrix multiplication with this and other highly-optimized routines such as the FFT are further detailed in Appendix C.

4. CONCLUSION:

We address the problem of having to manually choose among the numerous classes of structured linear maps by proposing the universal (expressive, efficient, and learnable) family of kaleidoscope matrices. We prove that K-matrices can represent any structured linear maps with near-optimal space and time complexity. Empirical validations suggest that K-matrices are a promising and flexible way to employ structure in modern ML; they can be used to reduce the need for hand-engineering, capture challenging latent structure, and improve efficiency in models. We are excited about future work on further hardware-optimized implementations of K-matrices, to fully realize the size and speed benefits of structured matrices on a broad array of real-world applications.

A RELATED WORK

A.1 STRUCTURED MATRICES IN MACHINE LEARNING:

Structured linear maps such as the DFT, the Hadamard transform and convolution are a workhorse of machine learning, with diverse applications including data preprocessing, random projection, featurization, and model compression. For example, the DFT is a crucial step in the standard filter bank speech preprocessing pipeline (Jurafsky & Martin, 2014) , and is commonly used when dealing with time series data in general (Panaretos & Tavakoli, 2013) . Fast random projection and kernel approximation methods rely on the fast Hadamard transform (Le et al., 2013; Yu et al., 2016) and convolution (Yu et al., 2015) , and convolution is a critical component of modern image processing architectures (Krizhevsky et al., 2012) as well as being useful in speech recognition (Zeghidour et al., 2018) and natural language processing (Wu et al., 2019) . Large learnable classes of structured matrices such as Toeplitz-like matrices (Sindhwani et al., 2015) and low-displacement rank (LDR) matrices (Thomas et al., 2018) have been used for model compression. However, despite their theoretical speedup, these structured matrix classes lack efficient implementations, especially on GPUs. Therefore, their use has largely been confined to small models (e.g. single hidden layer neural nets) and small datasets (e.g. CIFAR-10).

Butterfly matrices encode the recursive divide-and-conquer structure of the fast Fourier transform (FFT) algorithm. They were first used in numerical linear algebra for fast preconditioning (Parker, 1995) . The butterfly factorization is then generalized to encompass complementary low-rank matrices commonly encountered in solving differential and integral equations (Rokhlin & Tygert, 2006; Tygert, 2008; 2010b; a; Li et al., 2015; 2018) . In machine learning, butterfly matrices have been use to approximate the Hessian for fast optimization (Mathieu & LeCun, 2014) , and to perform fast random projection (Jing et al., 2017; Munkhoeva et al., 2018; Choromanski et al., 2019) . Dao et al. (2019) show that butterfly matrices can be used to learn fast algorithms for discrete transforms such as the Fourier transform, cosine/sine transform, Hadamard transform, and convolution.

A.2 SPARSE MATRICES:

Several classes of structured linear transforms are ubiquitous in modern deep learning architectures;

particularly widespread examples include convolution and multiheaded attention. Recently, attempts to impose sparsity on the neural network weights have been gaining traction. State-of-the art approaches of this type typically accomplish this by pruning small weights (either gradually during training (Zhu & Gupta, 2017) , or post-training (Han et al., 2016) ) or by training a dense network and then identifying "winning lottery tickets"-sparse subnetworks which may then be retrained from scratch with appropriate initialization (Frankle & Carbin, 2019) . Importantly, these approaches start from a dense network, and therefore training is expensive. There is also a more nascent line of work that aims to train unstructured sparse neural networks directly (Mocanu et al., 2018; Mostafa & Wang, 2019; Dettmers & Zettlemoyer, 2019; Evci et al., 2019) . These approaches maintain a constant network sparsity level throughout training, and use heuristics to evolve the sparsity pattern during training. One drawback is that the indices of the nonzero entries need to be stored in addition to the entry values themselves, which increases the memory required to store the sparse weight tensors. Another drawback is that these approaches to learn the sparsity pattern are based on intricate heuristics, which can be brittle. We note that these heuristic sparsification techniques could potentially be combined with our approach, to further sparsify the K-matrix factors.

A.3 SPEECH RECOGNITION FROM RAW AUDIO:

Numerous works focus on the problem of speech recognition from raw audio input, i.e. without manual featurization. SincNet (Ravanelli & Bengio, 2018 ) is a CNN-based architecture parameterized with sinc functions, designed so that the first convolutional layer imitates a band-pass filter. Zeghidour et al. (2018) formulate a learnable version of a filter bank featurization; their filters are initialized as an approximation of MFSC features and then fine-tuned jointly with the rest of the model. Sainath et al. (2015) proposed a powerful combined convolutional LSTM (CLDNN)-based model for learning from raw audio, using a large amount of training data. The WaveNet generative architecture (van den Oord et al., 2016) , based on dilated convolutions, has been adapted to speech recognition and can be trained on raw audio. Other approaches that can learn from raw audio can be found in (Palaz et al., 2013; Collobert et al., 2016; Ghahremani et al., 2016) . To our knowledge, the 14.6% PER achieved by our kaleidoscope + LSTM model on the TIMIT test set is the lowest error rate obtained by a model trained directly on the raw audio.

A.4 LEARNING PERMUTATIONS Permutation matrices find use in tasks such as matching and sorting (among many others). Techniques to obtain posterior distributions over permutations have been developed, such as the exponential weights algorithm (Helmbold & Warmuth, 2009) and the Gumbel-Sinkhorn network (Mena et al., 2018) .

Classifying images with permuted pixels is a standard task to benchmark the ability of RNNs to learn long range dependencies. Le et al. (2015) propose the Permuted MNIST task, in which the model has to classify digit images with all the pixels permuted. Many new RNN architectures, with unitary or orthogonal weight matrices to avoid gradient explosion or vanishing, have been proposed and tested on this task (Le et al., 2015; Arjovsky et al., 2016; Wisdom et al., 2016; Mhammedi et al., 2017; Trinh et al., 2018) . Standard gated RNN architectures such as LSTM and GRU have also been found to be competitive with these new RNN architectures on this task (Bai et al., 2018) .

B ADDITIONAL EXPERIMENTAL DETAILS B.1 SPEECH PREPROCESSING:

In this section, we fully describe our settings and procedures for the speech preprocessing experiments in Section 3.1.1, and present additional auxiliary baselines and results.

B.1.1 EXPERIMENTAL SETUP:

We evaluate our speech recognition models on the TIMIT speech corpus (Garofolo et al., 1993) , a standard benchmark for speech recognition. The input is audio (16-bit, 16 kHz .wav format), and the target is the transcription into a sequence of phonemes (units of spoken sound). Our evaluation metric is the phoneme error rate (PER) between the true phoneme sequence and the phoneme sequence predicted by our model. We use PyTorch (Paszke et al., 2017) , the Kaldi speech recognition toolkit (Povey et al., 2011) , and the PyTorch-Kaldi toolkit (Ravanelli et al., 2019) for developing PyTorch speech recognition models for all our experiments and evaluations.

B.1.2 MODEL AND EVALUATION:

Our baseline Bi-LSTM architecture is taken from the PyTorch-Kaldi repository.6 This is a strong baseline model that, to the best of our knowledge, matches state-of-the-art performance for models that use a single type of input featurization (Ravanelli et al., 2019) . The original Bi-LSTM model takes as input filter bank features. These are computed as follows: (i) the waveform is framed (split into chunks of 25 ms each that overlap by 10 ms each), (ii) the waveform is dithered (zero-mean Gaussian random noise is added), (iii) pre-emphasis is applied to amplify high frequencies, (iv) the Hamming window function (Harris, 1978) is applied, (v) the FFT is applied, and the power spectrum of the resulting (complex-valued) output is computed, (vi) the power spectrum (which has dimension 512) is mapped to the "mel scale" (which is a scale intended to mimic human auditory perception (Stevens et al., 1937) ) by multiplication with a specific banded matrix of dimension 512 × 23, and the entrywise logarithm of the output is taken (the 23 outputs are called the filters), and (vii) cepstral mean and variance normalization (Liu et al., 1993) is applied. Numerical hyperparameters of this procedure include the dither noise scale, the pre-emphasis coefficient, the Hamming window size, the number of mel filters, and more; we kept all these the same as the Kaldi/PyTorch-Kaldi defaults.

In contrast, our "K-matrix version" of the model takes as input the raw waveform, split into chunks the same way as before but with no normalization, dithering, or other preprocessing, which is then fed into a complex-valued kaleidoscope [(BB * ) 2 ] matrix. Similarly to the nonlinear steps in computing filter bank features, the logarithm of the power spectrum of the output (which has dimension 512) is then computed. This output is fed into the Bi-LSTM; the Bi-LSTM and kaleidoscope layer are trained together in standard end-to-end fashion. The Bi-LSTM architecture is not modified aside from changing the input dimension from 23 to 512; this (along with the ≈ 75K parameters in the kaleidoscope layer itself) results in approximately a 1.1M increase in the total number of parameters compared to the model that takes in MFSC features (a modest 8% relative increase). Total training time for our kaleidoscope-based architecture is 7% greater than that required for the model that uses MFSC features, not counting the time required to precompute the MFSC features; the FLOPs for inference-time are approximately 15% greater (mostly due to the larger dimension of the input to the Bi-LSTM; the kaleidoscope layer accounts for less than 0.5% of the total FLOPs).

As baselines, we also compare to inserting other types of linear transformations before the Bi-LSTM: fixed linear transformations (such as the fixed FFT, or no transform at all [i.e. the identity]), other trainable structured layers (low-rank, circulant, and sparse [using the sparse training algorithm of Dettmers & Zettlemoyer (2019)]), and a trainable unstructured (dense) linear layer. The kaleidoscope layer performs the best out of all such approaches. The fact that it outperforms even a dense linear layer with more parameters is particularly notable, as it suggests that the structural bias imposed by the K-matrix representation is beneficial for performance on this task. Full results are given in Table 5 . (Ravanelli et al., 2018) 13.8 12.3M

In our experiments, we grid search the initial learning rate for the "preprocessing layer" (if applicable) in {5e-5, 1e-4, 2e-4, 4e-4, 8e-4, 1.6e-3}, and fix all other hyperparameters (including the initial learning rates for the other parts of the network) to their default values in the PyTorch-Kaldi repository.

The model and any preprocessing layers are trained end-to-end with the RMSProp optimizer for 24 epochs (as per the defaults in PyTorch-Kaldi). For each model, we use the validation set to select the best preprocessing learning rate, while the final error rates are reported on the separate held-out test set. For all structured matrix baselines except circulant (which always has n parameters for an n × n matrix), the number of parameters in the structured matrices is set to equal the number of parameters in the butterfly layer, while the unconstrained matrix is simply a standard dense complexvalued square matrix. For all experiments with a trainable "preprocessing layer," we initialize the preprocessing matrix to represent the FFT (or approximate it as closely as possible [i.e. minimize the Frobenius error to the true FFT matrix], in the case of low-rank, sparse, and circulant), which we found to outperform random initialization.

B.1.3 EXTENSION: COMBINING MFSC AND KALEIDOSCOPE:

As an additional experiment, we sought to investigate whether combining the hand-engineered MFSC featurization pipeline and a learnable kaleidoscope layer (instead of replacing the former with the latter) could lead to accuracy gains. Specifically, in this experiment we first used the standard filter bank featurization pipeline described above, and trained end-to-end as usual. Then, we replaced the FFT step with a K-matrix initialized to the FFT, and made the weights of the Hamming window function and the mel filter bank matrix learnable as well (similarly to (Zeghidour et al., 2018) ). We fine-tuned the resulting architecture for an additional 10 epochs. The final test PER% attained by this "hybrid" model is 14.0 ± 0.3; the model has 14.4M parameters-a negligible increase over the 14.3M in the original architecture. Thus, by combining the manually encoded domain knowledge in the filter bank featurization and allowing this structure to be learnable rather than fixed, we are able to nearly match the state-of-the-art 13.8% accuracy on TIMIT. While this "hybrid" model certainly involves some hand-engineering, the state-of-the-art results use a concatenation of three different speech audio featurizations-MFSC, MFCC, and fMLLR-as the neural network input, along with a customized RNN architecture (LiGRU) specifically designed for speech recognition, and thus require a more complicated pipeline that is arguably even more hand-crafted. In our architecture, we use a kaleidoscope matrix in OBB (product of an orthogonal butterfly matrix, a diagonal matrix, and the transpose of another butterfly matrix) before and after each grouped 1x1 convolution. We place the second K-matrix after the batch norm and ReLU, to more closely mimic the original ShuffleNet design. The structure for each block is: K-matrix → 1x1 group conv → Batch norm, ReLU → K-matrix → 3x3 depthwise conv → Batch norm → 1x1 group conv. The K-matrices are trained along with the rest of the network, rather than being fixed.

B.2.2 EXPERIMENTAL SETUP:

We evaluate the CNN architectures on the image classification task of the standard ImageNet dataset (Russakovsky et al., 2015) . We use the standard data augmentation, training, and evaluation pipeline as in (Xie et al., 2017) . We train with SGD on 8 GPUs for 90 epochs, with a total batch size of 2048 and initial learning rate 0.8. For the 1.0 ShuffleNet g8 architecture, we reduce the total batch size to 1792 to fit into GPU memory, and correspondingly linearly scale the initial learning rate to 0.7. Other hyperparameters (e.g. learning rate schedule, weight decay, etc.) are kept the same as in the ShuffleNet paper (Zhang et al., 2018) . We use the training script from NVIDIA's deep learning examples repository.7 

B.2.3 ADDITIONAL RESULTS:

In Table 6 , we report top-5 classification accuracy on ImageNet, to complement the top-1 accuracies in Table 2 . In each setting, the total training time of our K-matrix approach is within 20% of the total training time of vanilla ShuffleNet.

In Figure 4 , we plot the loss and accuracy on the training set and validation set when we train 1.0 ShuffleNet g8, with either a fixed permutation (Shuffle) or a K-matrix for channel shuffling. Even though each K-matrix is a product of multiple (sparse) matrices, the model with K-matrices takes about the same number of training steps to converge as the baseline model does. One possible reason is that we constrain the K-matrices to be orthogonal (Section 2.4), thus avoiding vanishing or exploding gradients. The permuted CIFAR-10 dataset is constructed by applying a fixed permutation to every input. We choose to use the 2-D bit-reversal permutation,8 i.e., the bit reversal permutation on 32 elements is applied to the rows and to the columns. This permutation was chosen because it is locality-destroying: if two indices i, j are close, they must differ in a lower-order bit, so that the bit-reversed indices i , j are far. This makes it a particularly challenging test case for architectures that rely on spatial locality such as "vanilla" CNNs.

B.3.2 MODEL AND TRAINING:

We describe the model architectures used in Section 3.1 (those reported in Table 3 ).

Our model (K + CNN):

The model represents a fixed permutation P , parametrized as a K-matrix, to learn to recover the true permutation, followed by a standard ResNet18 architecture (He et al., 2016) .

Because of the simple decomposable nature of the butterfly factors (Section 2.1), our parameterization is easily extensible with additional techniques:

(i) We constrain each butterfly factor matrix in the K-matrix to be doubly-stochastic. For example, each 2 × 2 block in the butterfly factor matrix of block size 2 has the form a 1 -a 1 -a a , where a ∈ [0, 1]. We treat this block as a distribution over permutations, generating the identity 1 0 0 1 with probability a and the swap 0 1 1 0 with probability 1-a. Butterfly factor matrices with larger block sizes are constrained to be doubly-stochastic in a similar manner. In this way, a permutation is sampled for each butterfly factor matrix, and these permutations are composed to get the final permutation that is applied to the image. (ii) For each minibatch, the examples P x by applying permutation samples on the (permuted) inputs are fed into an additional unsupervised reconstruction loss 0≤i,j<n

(P x)[i + 1, j] -(P x)[i, j] (P x)[i, j + 1] -(P x)[i, j] 2 (1)

measuring total variation smoothness of the de-noised inputs. Such loss functions are often used in image denoising (Rudin et al., 1992) . A final regularization loss was placed on the entropy of P , which was annealed over time to encourage P to converge toward a sharper doubly-stochastic matrix (in other words, a permutation).

The model is trained with just the reconstruction loss to convergence before the standard ResNet is trained on top.

These techniques are applicable to the K-matrix as well as specialized methods for representing permutations such as Gumbel-Sinkhorn (Mena et al., 2018) and are important for recovering the true permutation. However, they are not applicable to a general linear layer, which showcases the flexibility of K-matrices for representing generic structure despite not being specially tailored for this task. We also remark that other classes of structured linear maps such as low-rank, circulant, and so on, are even less suited to this task than dense matrices, as they are incapable of representing all permutations.

Baseline architectures:

1. Fully connected (FC): This is a 3-layer MLP, with hidden size 1024 and ReLU nonlinearity in-between the fully connected layers. 2. Recurrent neural network (RNN): We use a gated recurrent unit (GRU) model (Cho et al., 2014) , with hidden size 1024. Many RNN architectures have been proposed to capture long-range dependency on permuted image dataset such as Permuted MNIST (Arjovsky et al., 2016) . Standard gated architectures such as LSTM and GRU have shown competitive performance on the Permuted MNIST dataset, and we choose GRU as a baseline since it has been reported to slightly outperform LSTM (Bai et al., 2018) . 3. CNN: We use the standard ResNet18 architecture, adapted to smaller image size of the CIFAR-10 dataset (changing stride from 2 to 1 of the first convolutional layer, and removing max-pooling layer that follows). 4. Dense + CNN: We add an additional linear layer (i.e. a dense matrix) of size 1024 × 1024 before the ResNet18 architecture. This dense layer can in theory represent a permutation, but cannot benefit from the additional techniques described above. 5. Baseline CNN (unpermuted): We use the standard ResNet18 architecture applied to the unpermuted CIFAR-10 dataset.

All models are trained for 200 total epochs, with the Adam optimizer. We use the standard learning rate schedule and weight decay from Mostafa & Wang (2019) . We use Hyperband (Li et al., 2017) to tune other hyperparameters such as the initial learning rate and annealing temperature.

B.4 SPEEDING UP DYNAMICCONV'S INFERENCE:



B.4.1 MODEL ARCHITECTURE:

We start with the DynamicConv Transformer architecture (Wu et al., 2019) , which is a variant of the Transformer architecture (Vaswani et al., 2017) where the self-attention in each layer is replaced with a light-weight DynamicConv module. We use the implementation from the Fairseq library (Ott et al., 2019) , 9 with PyTorch version 1.2.

The architecture of each layer of the decoder is: Linear → DynamicConv → Linear → LayerNorm → Encoder-decoder attention → LayerNorm → Linear → ReLU → Linear → ReLU → LayerNorm.

In every layer of the decoder, we replace the dense weight matrix in each of the four Linear layers with a K-matrix from the B class (i.e. a butterfly matrix).

B.4.2 TRAINING AND EVALUATION:

The models are trained from scratch using the training script from the Fairseq repository, with the same hyperparameters (optimizer, learning rate, number of updates, etc.) used in the DynamicConv paper (Wu et al., 2019) . We note that the DynamicConv model with K-matrices in the decoder trains slightly faster than the default DynamicConv model (both models are trained for 50,000 updates, which requires approximately 7% less time for the K-matrix model than for the default model).

To evaluate inference speed, we run the decoding script on the IWSLT-14 De-En test set in singlethreaded mode on a server Intel Xeon CPU E5-2690 v4 at 2.60GHz, and measure wall-clock time.

The test set contains 6750 sentences, with 149241 tokens. Following Wu et al. (2019) , we set the batch size to 1 and beam size to 1 for this evaluation.

B.4.3 ADDITIONAL COMPARISON WITH OTHER STRUCTURED MATRICES:

We additionally compare the speed-quality tradeoff of K-matrices with other classes of structured matrices, when used to replace the fully-connected layers of DynamicConv's decoder. We consider the following additional classes of structured matrices: low-rank, circulant, Toeplitz-like (Sindhwani et al., 2015) , ACDC (Moczulski et al., 2016) , Fastfood (Le et al., 2013) , and sparse. For classes with a variable number of parameters (e.g. low-rank, sparse), we set the number of parameters to match that of K-matrices. For sparse matrices, besides the result for an ensemble of 10 models (the default setting in the Fairseq repository), we also report the result for a single model, as that could have faster inference time (since ensembling/averaging sparse matrices produces a less sparse matrix).

In Figure 5 , we plot the tradeoff between translation quality (measured by BLEU score) and inference speed (sentences per second). Most classes of structured matrices produce similar translation quality (between 34.1 and 34.4 BLEU score). K-matrices have the second fastest inference time, only 7% slower than low-rank matrices. We note that low-rank matrices benefit from very well-tuned BLAS routines (matrix-matrix multiplication). Even though our implementation of K-matrix multiplication

is not yet highly optimized, it is already quite close to the speed of low-rank matrix multiplication at an equivalent parameter count.

C SPEED BENCHMARK AND IMPLEMENTATION DETAILS:

Each K-matrix (for fixed width and expansion), has an O(n log n) matrix-vector multiplication algorithm: sequentially multiply the input vector with each of the sparse factors. Our implementation of this simple algorithm is surprisingly competitive with optimized subroutines, both on GPU (e.g. for training) and on CPU (e.g. for inference). In Figure 6 , we compare the speed of multiplying by a K-matrix in class B (i.e. a butterfly matrix) against a specialized implementation of the FFT. We normalize the speed by the speed of dense matrix-matrix multiply (on GPU) or dense matrix-vector multiply (on CPU). On GPU, with input sizes n = 1024 and batch size 2048, the training time (forward and backward) of K-matrices matrix is about 3x faster than dense matrix multiply (GEMM from cuBLAS). For inference on CPU, the kaleidoscope fast multiplication can be one or two orders of magnitude faster than GEMV. Over a range of matrix sizes, our implementation is within a factor of 2-4x of specialized implementations of the FFT, a highly optimized kernel.

Our implementation is also memory-efficient. In the forward pass through the O(log n) sparse factors, we do not store the intermediate results, but recompute them during the backward pass. Therefore the activation memory required is O(bn) for an input batch size of b.

D SYNTHETIC MATRIX RECOVERY:

We directly validate Theorem 1 on well-known types of structured matrices used in machine learning. Given a structured matrix M, we attempt to represent M as closely as possible using K-matrices as well as the standard classes of structured matrices: sparse and low-rank. In Table 7 , we quantify the expressivity of each of these three methods, as measured by their ability to approximate a range of different structures. Results for "global minimum" of kaleidoscope matrices are obtained from the theoretical expressiveness results in Section I and Section J. Low-rank and sparse approximation have closed form solutions: truncating the SVD and keeping the largest-magnitude entries, respectively. We also report the results using SGD for kaleidoscope matrices to validate that good approximation with K-matrices can be obtained even from standard first-order optimization algorithms. Even with imperfect optimization, kaleidoscope matrices can still capture out-of-class target matrices better than low-rank and sparse matrices.

The target matrices are kaleidoscope, low-rank, sparse, convolution (i.e. circulant matrices), Fastfood (Le et al., 2013) , and entrywise random IID Gaussian matrix (to show the typical magnitude of the error). All target matrices M were randomly initialized such that E[M T M] = I. To find a kaleidoscope approximation with SGD, we used Hyperband to tune its learning rate (between 0.001 and 0.5).

E PROPERTIES OF THE BB * HIERARCHY:

Here, we justify why the definitions in Section 2.2 give rise to a hierarchy. We first make some basic observations about the parameterization. Observation E.1. An n × n matrix M ∈ BB * has 4n log n parameters.

Proof. M can be expressed as a product of 2 log n butterfly factor matrices of size n × n. Each of these factor matrices has 2 parameters per row, for a total of 2n parameters each. Hence, the total number of parameters is 4n log n.

Observation E.2. Let M be an n × n matrix in (BB * ) w e . Then, given an arbitrary vector v of length n, we can compute Mv with O(wne log(ne)) field operations.

Proof. Since M ∈ (BB * ) w e , we can decompose it as SE 1 E 2 . . . E w S T , where S is as given in Definition 2.4, and each E i is an en × en matrix in BB * . Therefore, to compute Mv, we can use associativity of matrix multiplication to multiply the vector by one of these matrices at a time. Since all of these factors are sparse, we use the naïve sparse matrix-vector multiplication algorithm (begin with a 0-vector and perform the corresponding multiplication and addition for each nonzero matrix entry). S (and thus S T ) have n NNZ. Therefore, matrix-vector multiplication by S or S T requires O(n) operations, which is dominated by the butterfly matrix-vector multiplication. Each E i can be further decomposed into 2 log(ne) matrices with at most 2ne non-zero entries each (by Observation E.1). Therefore, matrix vector multiplication by each E i requires O(ne log(ne)). Since there are w such E i , we require a total of O(wne log(ne)) operations. Now, we are ready to show that our definition of classes (BB * ) w e forms a natural hierarchy. First, we must argue that all matrices are contained within the hierarchy.

Lemma E.3. Let M be an arbitrary n × n matrix. Then M ∈ (BB * ) (2n-2) .

Proof. Corollary E.3 in Appendix K shows that any n × n matrix can be written in the form

M 1 M 1 * . . . M n-1 M n-1 * MM n M n * . . . M 2n-2 M n-2 *

, where M i , M i are orthogonal butterfly matrices and M is a diagonal matrix. We can combine D with M n to form another (possibly not orthogonal) butterfly matrix. This yields a decomposition of M as products of (possibly not orthogonal) butterfly matrices and their (conjugate) transposes, completing the proof.

Next, we argue that, up to a certain point, this hierarchy is strict.

Lemma E.4. For every fixed c ≥ 1, there is an n × n matrix M n (with n sufficiently large) such that

M n ∈ (BB * ) c+1 but M n ∈ (BB * ) c .

Proof. Given c, fix n to be a power of 2 such that c < n 4 log 2 n . For sake of contradiction, assume that every n × n matrix in (BB * ) c+1 is also in (BB * ) c . Let A be an arbitrary n × n matrix. From Lemma E.3, A ∈ (BB * ) (2n-2) . From our assumption, we can replace the first c + 1 BB * factors of A with c (potentially different) BB * factors and still recover A. We can repeat this process until we are left with c BB * factors, implying that A ∈ (BB * ) c . From Observation E.1, we require 4cn log n < n 2 (by our choice of n) parameters to completely describe A. This is a contradiction since A is an arbitrary n × n matrix, and therefore has n 2 arbitrary parameters. Hence, there must be some n × n matrix in (BB * ) c+1 that is not in (BB * ) c .

F ARITHMETIC CIRCUITS IN BB * HIERARCHY:

In this appendix, we prove our main theoretical result, namely, our ability to capture general transformations, expressed as low-depth linear arithmetic circuits, in the BB * hierarchy. This result is recorded in Theorem 1. Theorem 1. Let M be an n×n matrix such that matrix-vector multiplication of M times an arbitrary vector v can be represented as a be a linear arithmetic circuit C comprised of s gates (including inputs) and having depth d. Then, M ∈ (BB * )

O(d) O( s n ) .

To prove Theorem 1, we make use of the following two theorems. Theorem 2. Let P be an n × n permutation matrix (with n a power of 2). Then P ∈ BB * . Theorem 3. Let S be an n × n matrix of s NNZ. Then S ∈ (BB * )

4 s n 4 .

Theorem 2 is proven in Appendix G, and Theorem 3 is proven in Appendix I.

Proof of Theorem 1. We will represent C as a product of d matrices, each of size s × s , where s is the smallest power of 2 that is greater than or equal to s.

To introduce some notation, define w 1 , . . . w d such that w k represents the number of gates in the k'th layer of C (note that s = n + d k=1 w k ). Also, define z 1 , . . . z d such that z 1 = n and z k = w k-1 + z k-1 (z k is the number of gates that have already been used by the time we get to layer k).

Let g i denote the i'th gate (and its output) of C (0 ≤ i < s), defined such that:

g i = v i 0 ≤ i < n α j g i1 + β i g i2 n ≤ i < s

where i 1 , i 2 are indices of gates in earlier layers.

For the k'th layer of C, we define the s × s matrix M k such that it performs the computations of the gates in that layer. Define the i'th row of M k to be:

M k [i :] =    e T i 0 ≤ i < z k α i e T i1 + β i e T i2 z k ≤ i < z k + w k 0 i ≥ z k + w k For any 0 ≤ k ≤ d, let v k be vector v k = M k . . . M 2 M 1 v 0 .

We'd like to argue that v d contains the outputs of all gates in C (i.e, the n values that make up Mv). To do this we argue, by induction on k, that v k is the vector whose first z k+1 entries are g 0 , g 1 , . . . , g (z k -1) , and whose remaining entries are 0. The base case, k = 0 is trivial. Assuming this holds for the case k -1, and consider multiplying v k-1 by M k . The first z k rows of M k duplicate the first z k entries of v k-1 The next w k rows perform the computation of gates g z k , . . . , g (z k+1 -1) . Finally, the remaining rows pad the output vector with zeros. Therefore, v k is exactly as desired.

The final matrix product will contain all n elements of the output. By left multiplying by some permutation matrix P, we can reorder this vector such that the first n entries are exactly Mv. Hence, we are left to argue the position of PM d . . . M 2 M 1 within the BB * hierarchy. Each M k is a matrix with total 2w k + z k < 2s NNZ. From Theorem 3, we can, therefore, represent M k as a product of O( 1) matrices (of size 2s ) in BB * . From Theorem 2, P ∈ BB * . Note that s ≤ s < 2s, so s = Θ(s).

Our final decomposition will have O(d) BB * factors, and requires an expansion from size n to size 2s , or an expansion factor of O( s n ). Therefore, M ∈ (BB * )

O(d) O( s n )

, as desired. Remark F.1. By applying Observation E.2, we see that Theorem 1 gives an O(sd log s) matrix vector multiplication algorithm for M.

G PERMUTATIONS IN BB *:

In this appendix, we prove Theorem 2. In addition, we will also show that permutations are in B * B, where the set B * B is defined analogously to BB * (i.e. matrices of the form

M = M * 1 M 2 for some M 1 , M 2 ∈ B).

To prove Theorem 2, we decompose permutation matrix P into P = LR, with L ∈ B and R ∈ B * . Throughout the proof, we make use of the following definition.

Definition G.1. Let L be an n × n permutation matrix (n a power of 2). We say that L meets the 2 j balance condition if L can be divided into chunks of 2 j (with each chunk having all columns i such that i 2 j has the same value) such that for every 0 ≤ m < 2 j , each chunk has exactly one L[:, k] = e π k with π k ≡ m ( mod 2 j ). We say that L is modular-balanced if it meets the 2 j balance condition for each 2 ≤ 2 j ≤ n. Proof. We proceed by induction on n. The base case n = 2 is trivial. As our inductive hypothesis, we assume that all modular-balanced matrices of size n 2 × n 2 are butterfly matrices of size n 2 . From Definition 2.3, it is sufficient to show that L can be decomposed as:

0 1 2 3 4 5 6 7 (0) (1) (2) (3) (0) (1) (2) (3) L = B n L

L = B n L 1 0 0 L 2 L ,

where B n is a butterfly factor of size n and each L j is an n 2 × n 2 modular-balanced matrix.

Define L 1 and L 2 such that:

L 1 [i, j] = L[i, j] + L i + n 2 , j L 2 [i, j] = L i, j + n 2 + L i + n 2 , j + n 2 .

Note that since L is a permutation matrix (and thus has exactly one non-zero entry per column), at most one term of each of these sums can be non-zero.

For sake of contradiction, assume L 1 is not modular-balanced. Then, for some 2 j ≤ n 2 , there are two columns c 1 , c 2 such that c1 2 j = c2 2 j and such that indices of the non-zero entries of L 1 in columns c 1 and c 2 are the same modulo 2 j . However, from the definition of L 1 , this implies that the indices of the non-zero entries of L in columns c 1 and c 2 are also the same modulo 2 j , contradicting L being modular-balanced. Hence, L 1 is modular-balanced. An analogous argument (that instead considers columns

c 1 + n 2 , c 2 + n 2 of L)

shows that L 2 is also modular-balanced. To complete the proof, we must argue that B n is a butterfly factor of size n. Since each L i is modular-balanced, it is a permutation matrix. Therefore, L has exactly 1 non-zero entry in each of the first n 2 rows and columns from L 1 and exactly 1 non-zero entry in each of the second n 2 rows and columns from L 2 . Hence, L is a permutation matrix. Since both L and L are permutation matrices, B = L (L )

-1 must also be a permutation matrix. Therefore, we can view B as performing a permutation of the rows of L to get L.

Consider the i'th row of L , with 0 ≤ i < n 2 . There are two possible cases. Case 1:

L [i, :] = L[i, :]

In this case, the column of L with a non-zero entry in row i is in the left n 2 columns. The column of L with a non-zero entry in row i + n 2 must, therefore, be in the right n 2 columns, otherwise L would not satisfy the n 2 balance condition. Therefore,

L i + n 2 , : = L i + n 2 , : , so we set B[i, i] = B i + n 2 , i + n 2 = 1. Case 2: L [i, :] = L[i, :]

By the definition of L , L [i, :] = vL i + n 2 , : . In this case, the column of L with a non-zero entry in row i + n 2 must be in the left n 2 columns. By the n 2 balance condition of L, the column of L with a non-zero entry in row i must be in the right n 2 columns. Therefore, L i + n 2 , : = L [i, :], so we set B i, i + n 2 = B i + n 2 , i = 1. In both cases, the non-zero entries of B fall into the correct diagonal bands (the main diagonal, and the bands n 2 away). Hence, B is a butterfly factor of size n. Now, we consider the process of transforming P into a modular-balanced matrix. We make use of the following lemma. Lemma G.2. Let M be a k × k matrix with 1 non-zero entry per column, such that for each 0 ≤ m < k 2 , there are exactly 2 columns with non-zero entry in a row with index ≡ m mod k 2 . Then, there is a butterfly factor B k such that MB k = M , where M meets the k 2 balance condition.

EQUATION

Lemma H.2. Let A, B ∈ F n×n . If A ∈ (BB * ) w1 e and B ∈ (BB * ) w2 e then AB ∈ (BB * ) w1+w2 e . Proof. Let E A , E B ∈ F en×en be defined such that A = SE A S T , B = SE B S T (with S as in Definition 2.4). Then AB = S I n 0 0 0 en × en E A I n 0 0 0 en × en E B S T I n 0 0 0 E A ∈ (BB * ) w1 , I n 0 0 0 E B ∈ (BB * ) w2 by Lemma H.1. Hence, AB ∈ (BB * ) w1+w2 e by Definition 2.4. Lemma H.3. Let A 1 , . . . , A m ∈ F k×k . If A 1 , . . . , A m ∈ (BB * ) w e then Diag(A 1 , . . . , A m ) ∈ (BB * ) w+2 e . Proof. For each 1 ≤ i ≤ m, let E Ai ∈ F ek×ek be defined such that A i = SE Ai S T (with S as in Definition 2.4). Then     A 1 0 . . . 0 0 A 2 . . . 0 . . . . . . . . . 0 0 0 . . . A m     = SP     E A1 0 . . . 0 0 E A2 . . . 0 . . . . . . . . 0 0 0 . . . E Am     P T S T

where P is a permutation that that moves the first k rows of each E Ai (in order) into the top mk rows. From Theorem 2, P ∈ BB * , (and so is P T , also a permutation). Within the RHS block matrix, the decompositions of each E Ai can be done in parallel, requiring total width w. Hence,

Diag(A 1 , . . . , A m ) ∈ (BB * ) w+2 e

, as desired.

Remark H.4. If e = 1 in Lemma H.3, then P is unnecessary. Hence,

Diag(A 1 , . . . , A m ) ∈ (BB * ) w . Lemma H.5. Let A 1 , . . . , A m be k × k matrices in (BB * ) w e then m i=1 A i ∈ (BB * ) mw 4e .

Proof. For each 1 ≤ i ≤ m, let E Ai ∈ F ek×ek be defined such that A i = SE Ai S T (with S as in Definition 2.4). Note that E Ai ∈ (BB * ) w . Consider matrices of the form:

   I ek E Ai 0 0 0 I ek 0 0 0 0 0 0 0 0 0 0    M i ∈ F 4ek×4ek = I 2ek I 2ek 0 0 L    I ek 0 0 0 0 I ek 0 0 0 0 E Ai 0 0 0 0 0    S     I ek 0 0 I ek 0 0 0 0 0 0 0 0 0 I ek I ek 0     P 1 I 2ek 0 I 2ek 0 R .

Here, L and R compute the sum of the 2ek × 2ek matrices on the diagonal of SP 1 , where P 1 is a permutation swapping E Ai to the 4 th ek-block column. Note that S is the diagonalization of four matrices in (BB * ) w , so S ∈ (BB * ) w by Remark H.4. In addition, since each block in S is a butterfly matrix of size ek, S only uses butterfly factors up to size ek, so the outer factor matrices of sizes 4ek and 2ek in S are unused. Also note that L and R are butterfly factor matrices of size 4ek (or B (4ek) 4ek ), and P 1 is a butterfly factor matrix of size 2ek (or B (4ek) 2ek ). This allows us to fold the surrounding matrices L, P 1 , R into S, so M i ∈ (BB * ) w .

Through repeated application (m times) of the identity

EQUATION

From Lemma H.2, M ∈ (BB * ) mw . Finally, note that

m i=1 A i = SMP 2 S T

, where P 2 is a permutation that moves the first k columns of the second block-column of M to the left. P 2 can be folded into the final summation factor M m as follows:

    I ek 0 0 I ek 0 0 0 0 0 0 0 0 0 I ek I ek 0     P 1 I 2ek 0 I 2ek 0 R    0 I ek 0 0 I ek 0 0 0 0 0 I ek 0 0 0 0 I ek    P 2 =     0 I ek I ek 0 0 0 0 0 0 0 0 0 I ek 0 0 I ek     P 1 I 2ek 0 I 2ek 0 R (3) Hence, m i=1 A i ∈ (BB * ) mw

4e , as desired.

Lemma H.6. Let M be an invertible n × n matrix such that M ∈ B. Then M -1 ∈ B * .

Proof. We prove this in a series of steps.

First, let B k be an invertible butterfly factor of size k. Consider the method of computing B -1 k by performing Gaussian elimination on the matrix [B k |I k ] to obtain the matrix I k |B -1 k . By the form of B, non-zero entries within a row or column are always exactly k 2 positions apart. Therefore, the only row operations needed this Gaussian elimination are:

• Scaling a row by a constant factor c = 0

• Addition of a row to another row exactly k 2 rows apart Performing these operations on I k will only allow non-zeros on the main diagonal and k 2 diagonals away from the main diagonal. Hence, B -1 k is also a butterfly factor of size k.

Next, let B

(n)

k be an invertible butterfly factor matrix of size n and block size k. Its inverse is the block diagonal matrix formed by the inverses of each of its constituent butterfly factors. From above,

B (n) k -1

is also a butterfly factor matrix of size n and block size k.

Finally, consider M ∈ B.

M -1 = B (n) n B (n) n 2 . . . B (n) 2 -1 = B (n) 2 -1 B (n) 4 -1 . . . B (n) n -1 = B 2 (n) B 4 (n) . . . B n (n) ∈ B *

Finally, we include a closure result for the Kronecker product, another common matrix composition operation. Although Lemma H.7 is not directly used in the subsequent proofs, it allows for examples the results for the DFT to be lifted to higher-dimensional Fourier transforms. We also note that the closure bound in Lemma H.7 can be tightened in such cases (cf Ṙemark H.4). We use this lemma and the addition closure lemma to prove Theorem 3.

Proof of Theorem 3. We note that any s sparse matrix is the sum of s n matrices of at most n NNZ, and we appeal to Lemma H.5.

In the rest of the section we will prove Lemma I.1. We begin by defining two classes of matrices that will be used in our decomposition.

Definition I.1. An n × n matrix H is a horizontal step matrix if for every 0 ≤ i, i < n and 0 ≤ j ≤ j < n, if H[i, j] = 0 and H[i , j ] = 0, then j -j ≥ (i -i) mod n.

An n × n matrix V is a vertical step matrix if V * is a horizontal step matrix.

With this definition, the horizontal step matrix obeys a "Lipschitz-like" condition. Each column of a horizontal step matrix can have at most one non-zero entry, and given two non-zero columns k apart, the non-zero entry in the right column must be between 0 and k rows below the non-zero entry in the left column. Note that to show that a matrix is a horizontal step matrix, it is sufficient to argue that this condition holds for each pair of neighboring non-zero columns.

Similarly, each row of a vertical step matrix can have at most one non-zero entry, and given two non-zero rows k apart, the non-zero entry in the lower row must be between 0 and k columns to the right of the non-zero entry in the upper row. Lemma I.2. Let H be an n × n horizontal step matrix. Then H ∈ B.

Proof. We proceed by induction on n. The base case n = 2 is trivial. As our inductive hypothesis, we assume that all horizontal step matrices of size n 2 × n 2 are butterfly matrices of size n 2 . From Definition 2.3, it is sufficient to show that H can be decomposed as:

EQUATION

where Then, define H 1 and H 2 by:

H 1 = H 11 + H 21 H 2 = H 12 + H 22

For sake of contradiction, assume that H 1 is not a horizontal step matrix. Then, there are 0 ≤ i, i < n 2 , 0 ≤ j ≤ j < n 2 such that H 1 [i, j] = 0, H 1 [i , j ] = 0, and j -j < (i -i) mod n 2 . From our definition of H 1 , the non-zero entries in columns j and j of H are either (i -i) mod n 2 or n 2 + (i -i) mod n 2 , both of which are greater than j -j, rows apart. This contradicts H being a horizontal step matrix. Hence, H 1 must be a horizontal step matrix, as must H 2 from an analogous argument.

Next, we define D 1 , D 2 , D 3 , D 4 by:

D 1 [k, k] = 1 H 21 [k, :] = 0 0 otherwise D 2 [k, k] = 1 H 22 [k, :] = 0 0 otherwise D 3 [k, k] = 1 H 11 [k, :] = 0 0 otherwise. D 4 [k, k] = 1 H 12 [k, :] = 0 0 otherwise.

To finish the proof, we argue the correctness of the decomposition by equating arbitrary entries of each of the 4 corner submatrices. We begin with the upper left submatrix.

D 1 H 1 [i, j] = n 2 k=0 D 1 [i, k] • H 1 [k, j]

by definition of matrix multiplication

= D 1 [i, i] • H 1 [i, j] D 1 is a diagonal matrix = 1 (H21[i,:]=0) • (H 11 [i, j] + H 21 [i, j]

) by definition of D 1 and H 1

Here, we consider two cases:

Case 1: H 21 [i, j] = 0

Since H is a horizontal step matrix (and hence may have at most one non-zero entry per column), it follows that H 11 [i, j] = 0. In this case, the indicator function evaluates to 0, so D 1 H 1 [i, j] = 0 = H 11 [i, j], as desired.

Case 2: 

H 21 [i, j] = 0 If H 11 [i, j] = 0, then D 1 H 1 [i, j] = 0 = H 11 [i, j].

= 0, so D 1 H 1 [i, j] = H 11 [i, j], as desired.

In all cases, D 1 H 1 [i, j] = H 11 [i, j], so our decomposition correctly recovers the upper left corner of H. Analogous arguments show that the other three corners are also correctly recovered. Hence, our decomposition is correct, and by induction, H ∈ B.

Corollary I.3. Let V be a vertical step matrix. Then V ∈ B * .

1 2 3 4 S = 1 1 1 1 P 1 1 2 3 4 S 1 1 1 1 P 3 || 1 2 3 4 H 1 1 1 1 V || 1 1 1 1 P 2 1 1 1 1 V Figure 9: Decomposition of 4 × 4 sparse matrix S into P 1 HP 2 VP 3

Now, we use step matrices to prove Lemma I.1.

Proof of Lemma I.1. Given S, we decompose it as S = P 1 HP 2 VP 3 , where each P is a permutation matrix, H is a horizontal step matrix, and V is a vertical step matrix. For an example of this, see Figure 9 .

We first decompose S as S = P 1 S P 3 , where P 1 is the permutation that moves all 0 rows of S to the bottom and P 3 is the permutation that moves all 0 columns of S to the right.

Next, we further decompose S into S = HV as follows. Since S has s ≤ n NNZ, we can parameterize S by θ = {(c k , i k , j k ) : 0 ≤ k < s} such that S [i k , j k ] = c k , with the non-zero entries indexed in row-major order. Define matrix H by:

H[:, k] = c k • e i k 0 ≤ k < s 0 otherwise.

Define matrix V by:

V [k, :] = e T j k 0 ≤ k < s 0 otherwise.

pattern of Q 1 and the orthogonality of its columns. Moreover, the norm of each column of Q 1 R is 1. Therefore, Q 1 R is an orthogonal butterfly factor matrix, so M = QR ∈ (OBB) n-1 , as desired.

We now return to the proof of Lemma K.1

Proof of Lemma K.1. Given u ∈ C n (n a power of 2), let

u 0 = u[: n/2] ∈ C n/2 , u 1 = u[n/2 :] ∈ C n/2

denote the first and second halves of u.

To show that H ∈ OBB with inner diagonal matrix I, we proceed by induction. The base case for n = 2 is trivial. It suffices to show that there exist unitary butterfly factors L, R such that LHR has the form

I n/2 -2v 0 v * 0 0 0 I n/2 -2v 1 v * 1 for some unit vectors v 0 , v 1 ∈ C n/2 . Define (v 0 [i], v 1 [i]) =    u0[i] √ |u0[i]| 2 +|u1[i]| 2 , u1[i] √ |u0[i]| 2 +|u1[i]| 2 if |u 0 [i]| 2 + |u 1 [i]| 2 = 0

(1, 0) otherwise .

(5)

EQUATION

)

We choose

L = Diag(v * 0 ) Diag(v * 1 ) Diag(v 1 ) Diag(-v 0 ) and R = L * . L, R are (permuted) direct sums of blocks of the form v 0 [i] * v 1 [i] * v 1 [i] -v 0 [i]

, which are orthogonal by construction (via (5)). Hence, L ∈ B

(n) n and R ∈ ( B * ) (n) n . Further, LHR = Diag(v * 0 ) Diag(v * 1 ) Diag(v 1 ) Diag(-v 0 ) I -2 u 0 u 1 u 0 u 1 * Diag(v * 0 ) Diag(v * 1 ) Diag(v 1 ) Diag(-v 0 ) * = I -2 Diag(v * 0 ) Diag(v * 1 ) Diag(v 1 ) Diag(-v 0 ) u 0 u 1 u 0 u 1 * Diag(v * 0 ) Diag(v * 1 ) Diag(v 1 ) Diag(-v 0 ) * = I -2 v * 0 • u 0 + v * 1 • u 1 v 1 • u 0 -v 0 • u 1 w v * 0 • u 0 + v * 1 • u 1 v 1 • u 0 -v 0 • u 1 w *

, where • denotes the Hadamard product. From ( 6)

w[i] = |u 0 [i]| 2 + |u 1 [i]| 2 i ∈ [n/2] 0 i ∈ [n/2 : n]

Denoting the first half of this vector by w 0 ∈ C n/2 , we have

LHR = I -2w 0 w * 0 0 0 I ,

where w 0 2 = u 2 = 1. The result follows inductively.

As an immediate corollary, we can use Singular Value Decomposition to obtain a factorization for an arbitrary n × n matrix. Corollary K.3. Let M be an arbitrary n × n matrix. Then, M ∈ (OBB) 2n-1 , where all but one matrix in the decomposition is orthogonal (unitary).

Proof. By employing Singular Value Decomposition, we can decompose M as M = UΣV * , where U, V * are orthogonal and Σ is diagonal. By Lemma K.2, U, V * ∈ (OBB) n-1 , and trivially Σ ∈ OBB. Hence, M ∈ (OBB) 2n-1 . Note that Σ is the only matrix in the decomposition that is not orthogonal (unitary).

1 1 1 1 1 1 1 1 P 8 H 8 1 1 1 1 1 1 1 1 P * 8 = H 8 Figure 10: Block diagonalization of H 8

We argue that each of these 2×2 blocks can be decomposed into a diagonal matrix times an orthogonal matrix. Note that the butterfly factor matrices constructed in the proof of Lemma I.2 each have at most one non-zero entry per column. Hence, there are 4 cases to consider. Note that matrices with at most one non-zero entry are exhausted by Cases 1 and 2.

Case 1:

a 0 0 b = a 0 0 b D 1 0 0 1 O Case 2: 0 a b 0 = a 0 0 b D 0 1 1 0 O Case 3: a b 0 0 = √ a 2 + b 2 0 0 0 D a √ a 2 +b 2 b √ a 2 +b 2 b √ a 2 +b 2 -a √ a 2 +b 2 O , a, b = 0 Case 4: 0 0 a b = 0 0 0 √ a 2 + b 2 D b √ a 2 +b 2 -a √ a 2 +b 2 a √ a 2 +b 2 b √ a 2 +b 2 O , a, b = 0

In the last two cases, O is a 2 × 2 rotation matrix, which is commonly known to be orthogonal. Assume that we perform the above decomposition on all of the blocks of H 2 k in parallel, therefore expressing H 2 k = D O . We now have

H 2 k = P * 2 k D O P 2 k . By Lemma K.6, we can rewrite this as H 2 k = D P * 2 k O P 2 k . Note that P * 2 k O P 2 k

is the product of three orthogonal matrices, and thus orthogonal. Additionally, the construction of P 2 k ensures that P * 2 k O P 2 k is butterfly factor. 11 Hence, H 2 k can be expressed as the product of a diagonal matrix and an orthogonal butterfly factor, as desired. Now, we show that this decomposition of butterfly factors implies Lemma K.7. By performing this decomposition in parallel on each butterfly factor, we conclude that any butterfly factor matrix H

(n) 2 k of H can be decomposed as H (n) 2 k = D 2 k O (n) 2 k . 12

11 Conjugation by P 2 k is an isomorphism from 2 k × 2 k butterfly factors onto block diagonal matrices with 2 k-1 , 2 × 2 blocks. Therefore, conjugation by P -1 2 k = P * 2 k maps a block diagonal matrix to a butterfly factor.

We complete the argument by induction on n. The base case n = 2 holds by the observation about butterfly factor matrices above. Assume that any horizontal step matrix of size n 2 × n 2 can be expressed as a diagonal matrix times an orthogonal butterfly matrix. Now, consider the n × n horizontal step matrix H. From Lemma I.2, H can be expressed as

H = B (n) n H 1 0 0 H 2 ,

where H 1 , H 2 are n 2 × n 2 horizontal step matrices. By our inductive hypothesis,

EQUATION

n D 1 is a butterfly factor, and therefore can be expressed as

D n O (n) n . Therefore, H = D n O (n) n O 1 0 0 O 2 = D n O,

with O ∈ B, as desired.

Just as with the BB * hierarchy, the decomposition of vertical step matrices falls out as an immediate corollary to the horizontal step matrix proof. Corollary K.8. Let V be a vertical step matrix. Then we can decompose V = O * D, where D is a diagonal matrix and O * ∈ B * .

K.3.3 SPARSE MATRICES:

Now that we have argued about the decomposition of permutation and step matrices in the OBB hierarchy, we can leverage the construction from Lemma I.1 to argue about matrices with at most n NNZ.

Corollary K.9. Let S be an n × n matrix with at most n NNZ. Then, S ∈ (OBB) 4 .

Proof. We use the construction from Lemma I.1, along with Lemma K.7 and Corollary K.8, to express S as: 

S = O 1 O 1 P 1 D 2 O 2 H O 3 O 3 P 2 O 4 D 4 V O 5 O 5 P 3 , with each O i ∈ B, each O j ∈ B * ,

S = O 1 D 2 O 1 M 1 O 2 Õ3 M 2 Õ3 O 4 M 3 O 5 D 4 O 5 M 4

.

Note that each M ∈ OBB. Hence, S ∈ (OBB) 4 , as desired.

Just as in Appendix I, we would like to extend this orthogonal-based construction to capture matrices of general sparsity. To accomplish this, we introduce an addition closure lemma analogous to Lemma K.10 for the OBB hierarchy. Lemma K.10. Let A 1 , . . . , A m be k × k matrices in (OBB) w e then m i=1 A i ∈ (OBB) mw 4e .

With Lemma K.10, we arrive at the following Corollary on general orthogonal sparsity.

Corollary K.11. Let S be an n × n matrix with s NNZ. Then, S ∈ (OBB)

4 s n 4 .

Proof. Just as in the proof of Theorem 3, we accomplish this using a sum of s n matrices of at most n NNZ. For handling the sum of matrices, we need to appeal to Lemma K.10.

To conclude the argument, we give the proof of Lemma K.10.

Proof of Lemma K.10. For each 1 ≤ i ≤ m, let E Ai ∈ F ek×ek be defined such that A i = SE Ai S * (with S as in Definition 2.4). Note that E Ai ∈ (OBB) w . Consider matrices of the form:   

I ek 0 0 E Ai 0 I ek 0 0 I ek 0 0 -E Ai 0 I ek 0 0    M i ∈ F 4ek×4ek = √ 2 1 √ 2 I 2ek 1 √ 2 I 2ek 1 √ 2 I 2ek -1 √ 2 I 2ek O ∈ B (4ek) 4ek    I ek 0 0 0 0 I ek 0 0 0 0 E Ai 0 0 0 0 0    K    I ek 0 0 0 0 I ek 0 0 0 0 0 I ek 0 0 I ek 0    P ∈ B (4ek) 2ek

Note that K, a block diagonal matrix composed of matrices in (OBB) w , is itself in (OBB) w since

K = w j=1    I ek 0 0 0 0 I ek 0 0 0 0 O j 0 0 0 0 I ek    L j ∈ B    I ek 0 0 0 0 I ek 0 0 0 0 D j 0 0 0 0 0    Diagonal    I ek 0 0 0 0 I ek 0 0 0 0 O j 0 0 0 0 I ek    R j ∈ B * ,

where each O j is a ek × ek matrix in B, and each O j is a ek × ek matrix in B * . L w (the leftmost factor) is a block diagonal matrix composed of 4 ek × ek matrices in B. Therefore, we can fold O into this factor (since a butterfly factor in B

4ek

was not yet used in L w ) to conclude that OL w ∈ B. Similarly, since no btterfly factor from B (4ek) 2ek has been used in R 1 , we may fold P into R 1 to conclude that R 1 P ∈ B * . Finally, we address the scalar multiple of √ 2 by multiplying all entries of any diagonal matrix in the decomposition of K by √ 2. Hence, we may conclude that M i ∈ (OBB) w .

Through repeated application (m times) of the identity   

EQUATION

we see that

m i=1 M i =    I ek m i=2 E Ai 0 E A1 0 I ek 0 0 I ek -E Am + m-1 i=2 E Ai 0 E A1 0 I ek 0 0    M ∈ F 4en×4en . Therefore, M ∈ (OBB) mw . Next, we note that m i=1 A i = SM    0 I ek 0 I ek I ek 0 I ek 0 0 I ek 0 -I ek I ek 0 -I ek 0    Q S T .

We would like to show that we can fold Q into the rightmost OBB factor of M. The rightmost matrix in the decomposition of M is P. .

We now prove an asymptotic bound on the VC dimension of a ReLU network whose weight matrices are kaleidoscope matrices with bounded width and expansion. Proposition 3. Let F be the class of ReLU neural networks consisting of L layers, where each layer is a K-matrix with width and expansion bounded by some constant C. Suppose that the network has W total parameters. Let sign F denote the corresponding classification functions: {x → sign f (x) : f ∈ F}. Then this class has VC dimension: VCdim(sign F) = O(LW log W ).

We leverage the result from Thomas et al. (2018) for the case where the entries of the weight matrices interact multiplicatively, but with polynomially bounded degrees. This proof is similar to the VC bound for ReLU networks whose weight matrices are butterfly matrices (Dao et al., 2019) .

Proof. To use Theorem 3 of Thomas et al. (2018) , we simply need to check that the entries of the linear layer, as polynomials of the parameters, has degree at most c 1 m c2 l for some universal constant c 1 , c 2 > 0, where m l is the size of output of the l-th layer. If the network weight matrices are K-matrices with bounded width and expansion, each weight matrix is a product of at most c 3 log m l sparse factors, for some universal constant c 3 > 0. This means that the degree is polynomially bounded, which satisfies the condition of the theorem. Therefore the VC dimension is bounded to be almost linear in the number of parameters:

VCdim(sign F) = O(LW log W ).

M ARITHMETIC CIRCUIT PRIMER:

We give a quick overview of arithmetic circuits. This is a model of computation that has been studied for numerous computational problems (and is the basic model for algebraic complexity theory). For our purposes, we will exclusively focus on arithmetic circuits for the matrix-vector multiplication problem. For a more detailed exposition, the reader is referred to the standard book on this topic (Bürgisser et al., 2013) . Definition M.1 (Arithmetic Circuits). An arithmetic circuit that computes y = Ax (for A ∈ F m×n ) has n input gates (corresponding to x[0], . . . , x[n -1]) and m output gates (corresponding to y[0], . . . , y[m -1]). All the internal gates correspond to addition, subtraction, multiplication and division13 over the underlying field F. The circuit is also allowed to use constants from F for 'free.'

The definition of the internal gates can depend on A (as well as x of course). In other words, one can 'bake' the knowledge about A into the circuit.

The size s of a circuit is n plus the number of addition, multiplication, subtraction and division gates used in the circuit. The depth d of a circuit is the minimum number of layers such that all gates in a given layer take as its input gates from previous layers.14 

One drawback of arithmetic circuits (especially for infinite fields e.g. F = R, which is our preferred choice in this work) is that they assume operations over F can be performed exactly. In particular, it ignores precision issues involved with real arithmetic. Nonetheless, this model turns out to be a very useful model in reasoning about the complexity of doing matrix-vector multiplication for any family of matrices.

Perhaps the strongest argument in support of arithmetic circuits is that a large (if not an overwhelming) majority of matrix-vector multiplication algorithm also imply an arithmetic circuit of size comparable to the runtime of the algorithm (and the depth of the circuit roughly correponds to the time taken to compute it by a parallel algorithm). For example consider the obvious algorithm to compute Ax (i.e. for each i ∈ [m], compute y[i] as the sum n-1 i=0 A[i, j]x[j]). It is easy to see that this algorithm implies an arithmetic circuit of size O(nm) and depth O(log n). 15One thing to note about the arithmetic circuit above is that all the multiplications involve at least one input that is a constant from F (recall that we can assume that the entries of A are constants that can be used to build the circuit). This leads to the following important sub-class of arithmetic circuits: Definition M.2 (Linear Arithmetic Circuits). An arithmetic circuit is called a linear arithmetic circuit if it only uses addition, subtraction and multiplication. Further, every multiplcation has a fixed constant from F as at least one of its two inputs. In other words, all gates in the circuit are linear functions of their inputs (i.e. of the form ax + by for fixed constants a, b ∈ F).

Intuitively for the matrix-vector multiplication, it makes sense to consider linear arithmetic circuits since the final function we want to compute Ax is indeed a linear function of its inputs. For inifinite fields (e.g. F = R or F = C), it turns out that this is essentially without loss of generality: Theorem 4 ( (Bürgisser et al., 2013) ). Let F be an infinite field. Any (general) arithmetic circuit to compute Ax over F of size s and depth d can be converted into a linear arithmetic circuit of size O(s) and depth O(d).

The above result implies that for asymptotic considerations, linear arithmetic circuits for matrix-vector multiplication are equivalent to general arithmetic circuits. 16One important property of linear arithmetic circuits of depth d, which we will use in our arguments, is that such a circuit can be equivalently represented as product of d sparse matrices (see the proof of Theorem 1 for the precise derivation17 ).

As mentioned earlier, a vast majority of efficient matrix vector multiplication algorithms are equivalent to small (both in size and depth) linear arithmetic circuit. For example the FFT can be thought of as an efficient arithmetic circuit to compute the Discrete Fourier Transform (indeed when one converts the linear arithmetic circuit for FFT into a matrix decomposition,18 then each matrix in the decomposition is a butterfly factor, with each block matrix in each factor being the same). For an illustration of this consider the DFT with n = 4 as illustrated in Figure 11 . Finally, Figure 13 is representation of the arithmetic circuit of Figure 12 as a product of a butterfly matrix and (the bit-reversal) permutation. We note that our generic arithmetic circuit to decomposition into BB * is not as tight as in Figure 13 . One reason for the vast majority of existing efficient matrix vector algorithms leading to (linear) arithmetic circuits is that they generally are divide and conquer algorithms that use polynomial operations such as polynomial multiplication or evaluation (both of which themselves are divide and conquer algorithms that use FFT as a blackbox) or polynomial addition. Each of these pieces are well known to have small (depth and size) linear arithmetic circuits (since FFT has these properties). Finally, the divide and conquer structure of the algorithms leads to the circuit being of low depth. See the book of Pan (Pan, 2001) for a more elaborate description of this connection.

1 1 1 1 1 -i -1 i 1 -1 1 -1 1 i -1 -i

EQUATION

In fact, the recent work of De Sa et al. (De Sa et al., 2018) makes this fact explicit and presents the most general known structure on matrices that imply near-linear size linear arithmetic circuits for the corresponding matrix vector multiplication. Their work combines two separate classes of structures matrices-orthogonal polynomial transforms (Driscoll et al., 1997; Szegö, 1967) as well as matrices with low displacement rank (Kailath et al., 1979; Olshevsky & Shokrollahi, 2000) -and presents a linear class of linear arithmetic circuits to solve their matrix vector multiplication problem. We note that structured matrices with low displacement rank have been used to replace fully connected layers in some neural network architectures (Sainath et al., 2013; Thomas et al., 2018) .

:

A group of butterflies is known as a kaleidoscope.

Code that implements Kaleidoscope matrix multiplication is available at https://github.com/ HazyResearch/learning-circuits

The current state-of-the-art results fromRavanelli et al. (2018) use a concatenation of three different speech audio featurizations-MFSC, MFCC, and fMLLR-as the neural network input, along with a customized RNN architecture (LiGRU) specifically designed for speech recognition.

Despite our best effort, we were unable to reproduce the original accuracy reported byZhang et al. (2018), a problem similarly faced byZhao et al. (2019) andLyu et al. (2019).Zhao et al. (2019) use block Hadamard transform and pre-activation ShuffleNet, so their results are not directly comparable with those reported here.

BLEU score is a measure of translation quality; higher is better.

This open-source repository can be found at https://github.com/mravanelli/ pytorch-kaldi.

https://github.com/NVIDIA/DeepLearningExamples/tree/master/PyTorch/ Classification/RN50v1.5

The bit-reversal permutation reverses the order of the bits in the binary representation of the indices. For example, indices [0, 1, ..., 7] with binary representations[000, 001, ..., 111] are mapped to[000, 100, ..., 111], which corresponds to[0, 4, 2, 6, 1, 5, 3, 7]

This library can be found at https://github.com/pytorch/fairseq

Q is the (orthogonal/unitary) product of n -1 Householder reflections. R, the remaining upper triangular matrix after performing these reflections, is itself orthogonal/unitary, and therefore diagonal.

Note that a block diagonal matrix composed of orthogonal matrices is, itself, orthogonal.

Here we assume all the gates have two inputs.

The input layer corresponding to the input gates does not contriubte to the depth.

The claim on the depth follow from the fact that each of the sums n-1 i=0 A[i][j]x[j] can be computed in parallel. Further, the sum for each i ∈ [m] can be done in log 2 m depth by first computing the partial sumsA[i][2j ]x[2j ] + A[i][2j + 1]x[2j + 1] for all j ∈ [n/2]in parallel and recursively computing pair-wise sums till we are done.

This follows from the fact that by definition any linear arithmetic circuit is also an arithmetic circuit; the other direction follows from Theorem 4.

To the best of our knowledge, this connection was explicitly made by De Sa et al. (2018) though the connection seems to be folklore.

Using the conversion mentioned in the paragraph above.

Figures and Tables:

Figure 1: Visualization of the fixed sparsity pattern of the building blocks in BB * , in the case n = 16. The red and blue dots represent all the possible locations of the nonzero entries.

Figure 3 describes the pipeline and displays examples of permuted and unpermuted images.

Figure3: (a) (Left) Schematic describing permutation learning approach. The inputs are multiplied by a K-matrix and then fed into a CNN, from which the classification loss is computed. Separately, the input is permuted by a permutation matrix sampled from the distribution described by the Kmatrix, and a "smoothness" loss(Rudin et al., 1992) is computed from the result, as described in Appendix B.3. (b) (Right) Left panel: original (unpermuted) example images. Center panel: the permuted versions. Right panel: these images after then applying the permutation recovered by the K-matrix. The K-matrix is able to nearly unscramble the images into their unpermuted versions.

Figure 4: Loss and top-1 accuracy of 1.0 ShuffleNet g8 with either a fixed permutation (Shuffle) or a K-matrix for channel shuffling. The K-matrix model takes about the same number of training steps to converge as does the baseline model.

Figure 5: Tradeoff between translation quality (measured by BLEU score) and inference speed (sentences per second). K-matrices have the second fastest inference speed, only 7% slower than low-rank matrices.

Figure 6: Speedup of FFT and Kaleidoscope against dense matrix-matrix multiply (GEMM) for training, and against dense matrix-vector multiply (GEMV) for inference.

Figure 7: First step of decomposition of modular-balanced matrix L. Here, the red entries must be permuted into the main diagonal blocks.

Figure 8: First step of balancing 8 × 8 bit reversal permutation (a component of the 8 × 8 DFT). Red signifies edges that must be flipped.

Let A, B ∈ F n×n . If A ∈ (BB * ) w1 e and B ∈ (BB * ) w2 e then A⊗B ∈ (BB * ) w1+w2+6 e . Proof. Note that A ⊗ B = (A ⊗ I)(I ⊗ B) = P -1 (I ⊗ A)P(I ⊗ B), for some permutation P. By Lemma H.3, I⊗A and I⊗B are in (BB * ) w1+2 e , (BB * ) w2+2 e respectively. The result follows from combining with P ∈ BB * and Lemma H.2.I SPARSE MATRICES IN BB * HIERARCHYIn this appendix, we prove Theorem 3. First, we consider matrices with at most n NNZ.Lemma I.1. let S be an n × n matrix with at most n NNZ. Then, S ∈ (BB * ) 4 .

Figure 11: DFT of order 4.

Figure 12 represent the arithmetic circuit corresponding to FFT with n = 4.

Figure 13: Decomposition of DFT of Figure 11 via the arithmetic circuit of Figure 12.



Table content: <table><tr><td colspan="2">: TIMIT phoneme error rate (PER%) for different methods. Our kaleidoscope, raw-input</td></tr><tr><td colspan="2">version of the model (row 3) performs competitively with the original model trained on MFSC features</td></tr><tr><td colspan="2">(row 1), with only an 0.4% drop in PER. It significantly outperforms existing approaches that learn</td></tr><tr><td colspan="2">from raw audio, i.e. without handcrafted featurization (e.g. SincNet [row 2], which to our knowledge</td></tr><tr><td colspan="2">attains the previous state-of-the-art for learning from raw audio), and is only 0.8% less accurate than</td></tr><tr><td colspan="2">the overall state-of-the-art on TIMIT. 3 Additional comparisons are given in Appendix B.1.</td></tr><tr><td>Method</td><td>Test set PER% Raw audio input</td></tr><tr><td>MFSC features + LSTM</td><td>14.2</td></tr><tr><td>SincNet (Ravanelli et al., 2019)</td><td>17.2</td></tr><tr><td>Kaleidoscope + LSTM</td><td>14.6</td></tr></table>

Permuted CIFAR-10 validation set classification accuracy (%). Our kaleidoscope layer is able to nearly perfectly recover the latent structure, allowing a downstream CNN to approach the accuracy of a standard ResNet18 on the unpermuted dataset (last column).

Table content: <table><tr><td>Model</td><td colspan="5">FC RNN CNN Dense + CNN K + CNN Unpermuted</td></tr><tr><td colspan="2">Accuracy 61.2 57.8</td><td>73.7</td><td>84.4</td><td>93.6</td><td>94.9</td></tr></table>

Inference speed on the IWSLT-14 German-English translation task (test set). Using Kmatrices instead of dense matrices in the DynamicConv decoder linear layers results in 36% faster inference speed (measured on a single-threaded CPU with a batch size of 1 and beam size of 1).

Table content: <table><tr><td>Model</td><td colspan="4"># params BLEU Sentences/sec Tokens/sec</td></tr><tr><td>Transformer (Vaswani et al., 2017)</td><td>43M</td><td>34.4</td><td>3.0</td><td>66.4</td></tr><tr><td>DynamicConv Transformer (Wu et al., 2019)</td><td>39M</td><td>35.2</td><td>3.6</td><td>80.2</td></tr><tr><td>DynamicConv Transformer w/ K-matrices (ours)</td><td>30M</td><td>34.2</td><td>4.9</td><td>103.4</td></tr></table>

TIMIT phoneme error rate (PER%, ± standard deviation across 5 random seeds).

Table content: <table><tr><td>Model</td><td colspan="2">Test set PER% # Parameters</td></tr><tr><td>Low rank + LSTM</td><td>23.6 ± 0.9</td><td>15.5M</td></tr><tr><td>Sparse + LSTM</td><td>21.7 ± 0.9</td><td>15.5M</td></tr><tr><td>Circulant + LSTM</td><td>23.9 ± 0.9</td><td>15.4M</td></tr><tr><td>Dense + LSTM</td><td>15.4 ± 0.6</td><td>15.9M</td></tr><tr><td>FFT + LSTM</td><td>15.7 ± 0.1</td><td>15.4M</td></tr><tr><td>Identity + LSTM</td><td>20.7 ± 0.3</td><td>15.4M</td></tr><tr><td>Kaleidoscope + LSTM</td><td>14.6 ± 0.3</td><td>15.4M</td></tr><tr><td>MFSC features + LSTM</td><td>14.2 ± 0.2</td><td>14.3M</td></tr><tr><td>SincNet (Ravanelli et al., 2019)</td><td>17.2</td><td>10.0M</td></tr><tr><td>LiGRU</td><td/><td/></tr></table>

ShuffleNet is a convolutional neural network with residual (skip) connections that uses a permutation matrix to shuffle the channels after each grouped 1x1 convolution, sending the i-th channel to the (i mod g)-th group, where g is the total number of groups. The architecture for each residual block in ShuffleNet is: 1x1 group conv → Batch norm, ReLU → Permutation → 3x3 depthwise conv → Batch norm → 1x1 group conv. The permutation is fixed.Zhao et al. (2019) propose to instead use the Hadamard transform before and after each grouped 1x1 convolution to mix the channels. Note that the Hadamard transforms are placed before the batch normalization and ReLU layer (unlike the permutation matrix in the original ShuffleNet design). In particular, the architecture for each block is: Hadamard → 1x1 group conv → Hadamard → Batch norm, ReLU → 3x3 depthwise conv → Batch norm → 1x1 group conv. The Hadamard transform is fixed.

Table content: <table><tr><td>B.2 REPLACING CNN CHANNEL SHUFFLE</td></tr><tr><td>B.2.1 MODEL ARCHITECTURES</td></tr></table>

Top-5 classification accuracy of ShuffleNet on ImageNet validation set. We report results for different network width multipliers (number of channels), and for different kinds of matrices used for channel mixing. Using a learnable K-matrix in place of each fixed permutation (shuffle) or Hadamard matrix improves top-5 accuracy by up to 4.8%. Parameter counts are the same as in Table2.

Table content: <table><tr><td/><td colspan="4">Shuffle Hadamard Kaleidoscope (K.) K. vs. Shuffle</td></tr><tr><td colspan="2">0.25 ShuffleNet g8 68.6%</td><td>68.4%</td><td>73.4%</td><td>+4.8%</td></tr><tr><td>0.5 ShuffleNet g8</td><td>79.9%</td><td>79.2%</td><td>81.7%</td><td>+1.8%</td></tr><tr><td>1.0 ShuffleNet g8</td><td>86.0%</td><td>85.8%</td><td>86.8%</td><td>+0.8%</td></tr></table>

Expressiveness of different classes of structured matrices: Frobenius error of representing common structured matrices (columns) of dimension 256 using three structured representations of matrices with adjustable numbers of parameters. (Left group: Target matrices in the same class as the methods. Middle group: Target matrices with fixed number of parameters. Right: Random matrix to show typical scale of error.) Each method is allotted the same number of parameters, equal to a log n factor more than that of the target matrix. Low-rank and sparse matrices are unable to capture any structure outside their own class, while the minima for kaleidoscope matrices found via optimization better capture the actual structure for out-of-class targets better than the baselines.

Table content: <table><tr><td/><td>Method</td><td colspan="7">Target Kaleidoscope Low-rank Sparse Convolution Fastfood Random</td></tr><tr><td/><td colspan="2">Kaleidoscope</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td></tr><tr><td colspan="3">Global Min. Low-rank</td><td>14.9</td><td>0.0</td><td>10.8</td><td>14.6</td><td>11.6</td><td>15.5</td></tr><tr><td/><td>Sparse</td><td/><td>11.7</td><td>12.2</td><td>0.0</td><td>13.1</td><td>7.1</td><td>14.1</td></tr><tr><td>With SGD</td><td colspan="2">Kaleidoscope</td><td>0.0</td><td>0.01</td><td>8.0</td><td>0.0</td><td>5.1</td><td>14.5</td></tr></table>

H 1 , H 2 are n 2 × n 2 horizontal step matrices and each D k is a n

Table content: <table><tr><td>the four, n 2 × n 2 corner submatrices of H by:</td><td/><td>2 × n 2 diagonal matrix. Denote</td></tr><tr><td>H =</td><td>H 11 H 12 H 21 H 22</td><td>.</td></tr></table>

Otherwise, for sake of contradiction, suppose that H 21 [i, :] = 0. Then, two of the first n 2 columns of H would have non-zero entries n 2 rows apart, contradicting H being a horizontal step matrix. Hence, H 21[i, :]

Table content: <table/>

and each D k diagonal. Since P 2 is a permutation, by Corollary K.5, we can write it as Õ3 Õ3 for some Õ3 ∈ B * and Õ3 ∈ B. Moreover, noting that O 1 and O 5 are permutations, we make use of Lemma K.6 to re-express S as:

Table content: <table/>

Peer review (iclr2020):
Summary The authors introduce kaleidoscope matrices (K-matrices) and propose to use them as a substitute for structured matrices arising in ML applications (e.g. circulant matrix used for the convolution operation). The authors prove that K-matrices are expressive enough to capture any structured matrix with near-optimal space and matvec time complexity. The authors demonstrate that learnable K-matrices achieve similar metrics compared to hand-crafted features on speech processing and computer vision tasks, can learn from permuted images, achieve performance close to a CNN trained on unpermuted images and demonstrate the improvement of inference speed of a transformer-based architecture for a machine translation task. Review The overall quality of the paper is high. The main contribution of the paper is the introduction of a family of matrices called kaleidoscope matrices (or K-matrices) which can be represented as a product of block-diagonal matrices of a special structure. Because of the special structure, the family allows near-optimal time matvec operations with near-optimal space complexity for structured matrices which are commonly used in deep architectures. The proposed approach is novel. It gives a new characterization of sparse matrices with optimal space complexity up to a logarithmic term. Moreover, the proposed characterization is able to learn any structured matrix and matvec time complexity of the K-matrix representation is near-optimal matvec time complexity of the structured matrix. Even though in the worst-case complexity is not optimal, the authors argue that for matrices that are commonly used in machine learning architectures (e.g. circulant matrix in a convolution layer) the characterization is optimal. This results in a new differentiable layer based on a K-matrix that can be trained with the rest of an architecture using standard stochastic gradient methods. However, it is worth noting that the reviewer is not an expert in the field, and it is hard for him to compare the proposed approach with previous work. The paper is generally easy to follow. Even though the introduction of K-matrices requires a lot of definitions, they are presented clearly and Figure 1 helps to understand the concept of K-matrices. The experimental pipeline is also clear. Given the special structure of the family, the reviewer might guess that having K-matrices can slow down the training, i.e. it might require more epochs to achieve the reported results compared to baselines. Providing training plots might increase the quality of the paper. The experimental results are convincing. First, the authors show that K-matrices can be used instead of a handcrafted MFSC featurization in an LSTM-based architecture on the TIMIT speech recognition benchmark with only a 0.4% loss of phoneme error rate. Then, the authors evaluate K-matrices on ImageNet dataset. In order to do so, they compare a lightweight ShuffleNet architecture which uses a handcrafted permutation layer to the same architecture but with a learnable K-matrix instead of the permutation layer. The authors demonstrate the 5% improvement of accuracy over the ShuffleNet with 0.46M parameters with only 0.05M additional parameters of the K-matrix and the 1.2% improvement of accuracy over the ShuffleNet with 2.5M parameters with only 0.2M additional parameters of the K-matrix. Next, the authors show that K-matrices can be used to train permutations in image classification domains. In order to demonstrate so, they take the Permuted CIFAR-10 dataset and ResNet-18 architecture, insert a trainable K-matrix at the beginning of the architecture and compare against ResNet-18 with an inserted FC-layer (attempting to learn the permutation as well) and ResNet-18 trained on the original, unpermuted CIFAR-10 dataset. With K-matrix, the authors achieve a 7.9% accuracy improvement over FC+ResNet-18 and only a 2.4% accuracy drop compared to ResNet-18 trained on the original CIFAR-10. Finally, the authors demonstrate that K-matrices can be used instead of the decoder’s linear layers in a Transformer-based architecture on the IWSLT-14 German-English translation benchmark which allows obtaining 30% speedup of the inference using a model with 25% fewer parameters with 1.0 drop of BLEU score. Overall, the analysis and the empirical evaluations suggest that K-matrices can be a practical tool in modern deep architectures with a variety of potential benefits and tradeoffs between a number of parameters, inference speed and accuracy, and ability to learn complex structures (e.g. permutations). Improvements 1. Even though K-matrices are aimed at structured matrices, it would be curious either to empirically compare K-matrices to linear transformations in fully-connected networks (i.e. dense matrices) or to provide some theoretical analysis. 2. Section 3.3 argues that K-matrices allow to obtain an improvement of inference speed, however, providing the results of convergence speed (e.g. training plots with a number of epochs) will allow a better understanding of the proposed approach and will improve the quality of the paper.