BLAST: Block-Level Adaptive Structured Matrices for Efficient Deep Neural Network Inference

Authors: changwoo  lee
soo min kwon
qing  qu
hun-seok  kim
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 17:06:04

Abstract: Large-scale foundation models have demonstrated exceptional performance in language and vision tasks. However, the numerous dense matrix-vector operations involved in these large networks pose significant computational challenges during inference. To address these challenges, we introduce the Block-Level Adaptive STructured (BLAST) matrix, designed to learn and leverage efficient structures prevalent in the weight matrices of linear layers within deep learning models. Compared to existing structured matrices, the BLAST matrix offers substantial flexibility, as it can represent various types of structures that are either learned from data or computed from pre-existing weight matrices. We demonstrate the efficiency of using the BLAST matrix for compressing both language and vision tasks, showing that (i) for medium-sized models such as ViT and GPT-2, training with BLAST weights boosts performance while reducing complexity by 70\% and 40\%, respectively; and (ii) for large foundation models such as Llama-7B and DiT-XL, the BLAST matrix achieves a 2x compression while exhibiting the lowest performance degradation among all tested structured matrices. Our code is available at https://github.com/changwoolee/BLAST.

Main paper: 
 introduction: 
 Foundation models built on large deep neural networks (DNNs) have demonstrated remarkable performance in vision and language tasks. However, the size of these large networks poses both computational and storage challenges, especially in resource-constrained environments such as edge devices. The size of a single DNN often exceeds the capacity of the supporting hardware devices [1] [2] [3] [4] [5] . For example, Llama-70B [1] demands at least 140GB of memory solely for loading its weights in half-precision floating point representation, while the state-of-the-art commercial GPU only accommodates 80GB of memory. Furthermore, inference with these networks involves numerous dense matrix-vector operations, which can be limiting when computing power is constrained.Fortunately, large (overparameterized) DNNs often exhibits parameter redundancy, where the intrinsic dimension of the weights is much lower than the ambient dimension. As such, the weights should be structured, possessing hidden properties such as low-rankness [6] [7] [8] [9] or sparsity [10, 11] . Hence, it is possible to replace (or factorize) these dense existing weight matrices with structured ones without degrading performance [10] [11] [12] . However, using structured matrices that do not align with the true underlying structure of the weight matrices can result in significant performance degradation. We demonstrate this point in Figure 1 where we attempt to capture the structure of a diffusion model transformer (DiT) [13] using the low-rank structure to generate synthetic images. In Figure 1 , we compress the model's linear layers by approximately 50% of the total number of parameters using low-rank weight matrices via singular value decomposition (SVD) and generate images with the compressed model (see Section 4.2 and Appendix C.3 for details). As shown in Figure 1 (middle), simply using the low-rank structure introduces unwanted artifacts in the generated images. [13] starting from the same noise vectors and a deterministic solver. The original model is compressed by 50% through BLAST or low-rank matrices and re-trained for 10 epochs on ImageNet. The images from the model compressed via BLAST preserve the quality of the images of the original model, whereas the images generated by the low-rank model contain more undesired artifacts.To address this issue, many flexible structures for modeling DNN weights have been proposed to minimize the misalignment between imposed and true low-dimensional structures. For example, Dao et al. [14] proposed the Monarch matrix, a specific type of Block Low-Rank (BLR) structure [15] , in which all blocks share the same rank, intended for use in the linear layers of transformers [16] . Matrix multiplication with a Monarch matrix can be performed efficiently using batched matrix multiplication routines. Additionally, Chen et al. [17] investigated a block sparse plus low-rank structure. However, all of these methods still suffer from the fact that the underlying structure of each weight matrix is not known a priori. By imposing one of these structures, performance degradation may still occur due to misalignment. Recently, Lee and Kim [12] introduced a data-driven design called Generalized Block Low-Rank (GBLR). This approach employs multiple rank-1 blocks with various sizes and locations learned from data via differentiable masks. Unfortunately, the GBLR matrix is optimized for custom-designed hardware, as the learned block patterns are random. It has limited usability on general GPUs as the computation of GBLR matrices does not accelerate well on them.In this work, we introduce the Block-Level Adaptive Structured (BLAST) matrix, a versatile and efficient design tailored to uncover various low-dimensional structures in the weight matrices of DNNs for accelerated inference on GPUs. Our matrix structure leverages shared bases across block matrices with block-wise diagonal coupling factors. This structure encapsulates different structures such as low-rank, block low-rank, block-diagonal matrices, and their combinations. BLAST matrices can be applied to the training scenario from scratch or compression after training. For training from scratch, we let the linear layers of the DNN to directly adopt the BLAST structure and learn its factors from data. The factors of the BLAST matrix are constructed to have well-defined gradients, allowing them to be optimized using popular methods like stochastic gradient descent (SGD) or Adam [18] . For compressing existing weights, we propose a factorization algorithm to learn the BLAST factors from pre-trained weights. The compression performance can be further improved by updating the BLAST factors using data, a process we call "re-training".We demonstrate the efficiency of BLAST by training Vision Transformers (ViT) [19] and GPT-2 [20] from scratch on various datasets, showing that it can reduce complexity by 70% and 40%, respectively. We also compress existing ViT and Diffusion Transformer (DiT) [13] models with BLAST matrices by 70% and 50%, respectively, demonstrating that BLAST compression (and re-training) achieves higher accuracy / quality compared to existing methods for ViT and DiT (see Figure 1 ). For the language tasks, we compress Llama-7B [1] by 50% via BLAST and re-train on 0.49B tokens, showing the lowest accuracy degradation with significant inference speedup on a NVIDIA A100 GPU. Overall, our contributions can be summarized as follows: • We propose a novel block-structured matrix called BLAST that encompasses a wide range of matrix structures, allowing for faster matrix multiplication. Various existing structured matrices such as Low-Rank, Monarch [14] , and Block Diagonal matrices can be expressed using the BLAST matrix. • We provide gradient descent-based methods to find the BLAST factors for DNN weights. We empirically show that standard DNN training with the BLAST weight matrices effectively recovers the original accuracy while achieving up to a 70% reduction in computational complexity. • In cases where pre-trained dense weights are available, we propose a preconditioned gradient descent factorization algorithm to decompose the weights to BLAST factors for compression and further re-training. Our experimental results show that pre-trained foundation models for vision or language tasks can be compressed by 50% using BLAST matrices.Ui,j V T i,j Ai,j = Ui,jV T i,j Monarch [14] U V T A = U V T Low-Rank Di,i Ai,j = Di,i if i = j, O o.w. Block-Diagonal A = r ′ k=1 ukv T k GBLR [12] r Ui V T j Si,j n b r Ai,j = UiSi,jV T j BLAST (Proposed)Notation and Organization. We use σ 1 (X) to denote the largest singular value of the matrix X.The notation ⊙ indicates Hadamard product.The rest of the paper is organized as follows. In Section 2, we introduce the BLAST matrix and discuss its properties. In Section 3, we propose a methodology to train/compress DNNs with BLAST weight matrices. In Section 4, we demonstrate the effectiveness of the BLAST weights in improving efficiency without noticeable accuracy degradation. We discuss related works in Section 5, and conclude in Section 6. 
block-level adaptive structured (blast) matrix: 
 Consider a square matrix 1 A ∈ R n×n for some n ∈ N, which has an unknown intrinsic lowdimensional structure. We first equally partition the matrix A into b × b blocks of size p × p where b, p ∈ N are constants such that n = bp:EQUATIONThen, the BLAST matrix parameterizes each block matrix A i,j using three factors:EQUATIONwhere U i , V j ∈ R p×r are the left and the right factors, respectively, and S i,j = diag(s i,j ) is an r × r diagonal matrix whose diagonal entries are s i,j ∈ R r . We provide a visual representation on the rightmost side of Figure 2 , and illustrate how this structure differs from other types of matrices. While the BLAST structure may appear similar to SVD, there are two notable differences: (i) the left and right factors do not need to be orthonormal, and (ii) the diagonal entries do not need to be positive. These distinctions make it more flexible in capturing different types of low-rank structures.As illustrated in Figure 2 , the BLAST matrix also comes with two unique properties:• Factor Sharing: The left factor matrix U i of size rp is shared across b blocks at the i th row, i.e., A i,1 , . . . , A i,b . Likewise, the right factor V j is shared across the blocks at the j th column. On the other hand, the diagonal factor s i,j of size r is specific to each block A i,j . Hence the total number of parameters of an n × n BLAST matrix with b × b number of blocks of rank r is 2rpb + rb 2 = 2nr + rb 2 . This reduces the number of parameters b times by enforcing the blocks at the same row or column share the same bases. • Individual Diagonal Factors: The individual diagonal factors of each block matrix are the source of the adaptivity and flexibility of the BLAST matrix. By changing the values of the diagonal factors, the BLAST matrix can encompass a wide variety of matrix structures. These factors can be estimated using gradient descent, since s i,j is a real-valued vector and A i,j = U i diag(s i,j )V T j is linear to s i,j .Low-Rank Matrices as Special Cases of BLAST To demonstrate how the BLAST matrix can capture different types of structures, we present an example showing how the BLAST matrix can encompass a low-rank matrix. Consider the case where all the diagonal factors are ones, i.e., s i,j = 1 r for all i, j = 1, 2, . . . , b. Then, we can write the block matrix as follows:U V T =     U 1 U 2 . . . U b     [V 1 V 2 • • • V b ] =      U 1 V T 1 U 1 V T 2 • • • U 1 V T b U 2 V T 1 U 2 V T 2 • • • U 2 V T b . . . . . . . . . . . . U b V T 1 U b V T 2 • • • U b V T b      .Hence, if the true underlying structure is low-rank, we can expect the BLAST matrix to learn this specific structure. Similarly, we show in Section A.1 that the BLAST matrix can construct low-rank, block-diagonal, and block low-rank matrices through different diagonal parameters. A combination of these canonical structured matrices, such as a low-rank with block-diagonal matrix, can also be achieved by simply concatenating the factors of each matrix.Algorithm 1 BLAST Matrix-Vector Product Require: U , V , s, x 1: [x T 1 , x T 2 , • • • , x T b ] T ← x 2: for j = 1, 2, . . . , b do ▷ #Parallel 3:z j ← V T j x j 4: end for 5: for i = 1, 2, . . . , b do ▷ #Parallel 6:y i ← U i b j=1 s i,j ⊙ z j 7: end for 8: return y ← [y T 1 , . . . , y T b ] TMatrix Multiplication DNNs involve numerous matrix-vector (matrix-matrix) multiplications in the form of y = Ax (Y = AX). Algorithm 1 depicts the BLAST matrix-vector multiplication procedure. Consider the partitioned input vectorx = [x T 1 , x T 2 , • • • , x T b ] T and the partitioned output vector y = [y T 1 , y T 2 , • • • , y T b ]T . The i th partitioned output vector y i is then computed by the sum of the b block-wise matrix-vector multiplications along j = 1, . . . , b:EQUATIONThe number of multiplications required to perform the matrix-vector multiplication y = Ax is (2n + b 2 )r. The matrix multiplication z j = V T j x j , j = 1, . . . , b is computed once and shared across i = 1, . . . , b, whereas the matrix multiplications in Line 3 and Line 6 of Algorithm 1 can be executed in parallel, e.g., by torch.bmm in PyTorch [21] . An implementation of Algorithm 1 for general matrix or tensor inputs can be found in Appendix A. 
applications of blast matrices: 
 There are two main applications of BLAST matrices: (i) training from scratch with the BLAST structure and (ii) compression of pre-trained weights using BLAST factorization.To train a DNN on a dataset, parameters are typically initialized randomly and updated through stochastic gradient descent. In this setting, BLAST can replace dense weights to learn structures from the training data. Instead of using random dense weight matrices, the model is initialized with random BLAST factors U i , V j , s i,j . Since the forward and the backward path of the linear layer involving the weight matrix is composed of three linear operations as in Equation (3), the derivatives of the minibatch loss can be back-propagated by automatic differentiation frameworks [21] . Hence, all of the trainable parameters of BLAST can be updated using conventional optimizers (e.g., Adam [18] or AdamW [22] ) without additional treatment.BLAST Factorization via Gradient Descent Given pre-trained dense weights of a DNN, we can compress the weights using BLAST matrices. Let A denote the weight matrix and A i,j denote its blocks. We estimate the BLAST factors of A i,j by finding the factors of the BLAST matrix that minimize the Frobenius norm error between the original weight matrix and the BLAST structure:EQUATIONwhere * denotes the collection of all b components along the axis. This problem shares many characteristics with the classical matrix factorization problem [23] [24] [25] [26] , and hence we can solve for the factors using alternating gradient descent starting from small random initialization (e.g., Line 1 of Algorithm 2) [27, 8] . That is, the k th gradient descent step is composed of three alternating updates with a step size η > 0:EQUATIONEQUATIONEQUATION, sEQUATIONWith properly chosen step sizes, Equations ( 5) to (7) always decrease the loss value whenever the current variables do not have any infinite entries and the gradient is non-zero. Using notationsV (k) i = S (k) i,1 V (k)T 1 • • • S (k) i,b V (k)T b T and Ū (k) j = (U (k+1) 1 S (k) 1,j ) T • • • (U (k+1) b S (k) b,j ) T Tto indicate the concatenation of the right and left factors scaled by the diagonal components, the loss is monotonically non-increasing as in the following theorem. Theorem 1. Let A i,j ∈ R p×p be a target block andU (k) i , V (k) j∈ R p×r , and s (k) i,j ∈ R r be factors of a block in the BLAST matrix to be optimized. With the step sizes 0< η U (k) i ≤ 1/σ 1 V (k)T i V (k) i , 0 < η V (k) j ≤ 1/σ 1 Ū (k)T j Ū (k) j , 0 < η s (k) i,j ≤ 1/σ 1 (U (k+1)T i U (k+1) i ) ⊙ (V (k+1)T j V (k+1) j) , the gradient descent updates in Equations ( 5) to (7) monotonically non-increase the loss:ℓ(U (k+1) * , V (k+1) * , s (k+1) * , * ) ≤ ℓ(U (k) * , V (k) * , s (k) * , * ).The proof of Theorem 1 is in Section B, which is an application of the descent lemma from classical optimization theory.Blast Factorization via Preconditioned Gradient Descent Recall that in order to estimate the BLAST factors given a pre-trained weight A, we need to choose a rank r. Since we do not know the rank of A a priori, we may have to overestimate the rank. However, overestimating the rank may slow down the convergence rate for solving Equation (4). To illustrate this, we performed an empirical analysis of the convergence behavior on a synthetically generated target low-rank matrix A, whose dimension is 256 × 256 with a true rank of r * = 8. For the analysis, we computed the factors of a BLAST matrix with b = 16 for various values of r. We used linearly decaying step sizesη (k) = η U (k) i = η V (k) j = η s (k) i,j. When r = r * = 8 (the ranks of the left and right factors U * , V * match the actual rank of A), gradient descent finds a low-rank solution with minimal error within 30 iterations, as shown by the blue curve in Figure 3 -left. However, in the case of r = 32 > r * where the BLAST factors are overparameterized, we observed slower convergence and a substantial residual Algorithm 2 Preconditioned BLAST Factorization (see Equations ( 8) and ( 9))EQUATIONi,j -η (k) • P (k) si,j W (k) i,j s i,j -diag U (k+1)T i A i,j V (k+1) j , ∀i, j = 1, . . . , b 6: end for 7: return U (K) * , V (K) * , s (K) * , *error after 100 steps as shown by the blue curve in Figure 3 -right. This behavior is consistent with previous observations of slow convergence in ill-conditioned matrix factorization problems [23, 24] .The convergence rate of solving the overparameterized low-rank factorization by gradient descent can be improved via inexpensive preconditioners [23, 24] which effectively decrease the condition number at each iteration. Inspired by the preconditioned gradient descent for low-rank factorization, we generalize the idea to solve our problem by multiplying preconditioning matrices to the gradients in Equations ( 5) to (7) . We summarize the preconditioned gradient descent method for the BLAST factorization in Algorithm 2, where the following preconditioning matrices are used:EQUATIONEQUATIONδ is proportional to the square root of the error in Equation ( 4). The derivations are presented in Appendix A.2. Figure 3 shows that preconditioning improves the convergence of the overparameterized BLAST factorization. The preconditioned gradient descent (yellow curve) finds the points with low error after 100 steps, whereas the gradient descent without preconditioning fails to achieve a small error. More empirical studies on the BLAST factorization with or without preconditioning can be found in Appendix D.1. We summarize the compression procedure in Algorithm 2. The computational complexity of this compression algorithm is O(nr 2 + r 3 ) where the cubic dependency on r is from the matrix inversion steps. We emphasize that these matrix inversion steps are not substantial computational bottlenecks because r is smaller than n as in prior work [23, 24] .After BLAST compression outlined in Algorithm 2, the BLAST factors can be used directly for inference to save both storage and computational costs. However, we observe that, instead of using the estimated factors directly, using them as initial points and refining the estimates by re-training the model with BLAST factors can further improve the performance of the compressed model. We refer to this process as "re-training" after BLAST compression. Accuracy (%) Relative FLOPs (%) Dense ViT-Base 78.7 100 Low-Rank 78.9 33.5 Monarch [14] 78.9 33.5 Gaudi-GBLR [12] 78.  
experimental results: 
 We evaluate the BLAST matrix under two settings: (i) training from scratch with random initialization in the BLAST format, and (ii) re-training after compressing the dense weights to BLAST matrices via Algorithm 2. We compare the performance of BLAST with both non-adaptive and adaptive structured matrices. Among the non-adaptive approaches, we include low-rank (LR) matrices, Monarch for block low-rank (BLR) [14] , and Pixelfly [17] or Block-Diagonal for block sparse matrices. For the adaptive and learnable structured matrix category, we evaluate Gaudi-GBLR [12] . We report the number of floating point operations (FLOPs) by counting the number of multiplications. The BLAST matrix with b × b number of blocks is denoted by BLAST b . We used the same hyperparameter r for every target weight matrix by setting it to meet the computational budget of the DNN. All experimental details can be found in Appendix C. Image Classification We train the reduced-size Vision Transformers (ViT) [19] with BLAST 3 (BLAST with b = 3) weight matrices on CIFAR-10, 100 [28] , and ImageNet-1k [29] for 310 epochs from random initialization, and compare with other structured matrices. In the CIFAR-10 and CIFAR-100 benchmarks, BLAST outperforms several non-adaptive baselines, such as Low-Rank, Pixelfly, and Monarch with higher accuracy at the same FLOPs complexity (Figure 4 ). Gaudi-GBLR presents the most favorable accuracy-to-FLOPs tradeoff due to its capability of learning the adaptive resource/budget allocation for each weight matrix, which is a feature that our BLAST setting lacks in this particular evaluation (as we force it to use the same r for all matrices). However, in the context of ImageNet-1k in Table 1 , weight matrices trained using BLAST with b = 3 attain the highest levels of accuracy with the least FLOPs. This superior performance of BLAST (despite the common r for all matrices) over Gaudi-GBLR can be attributed to its simpler training process with fewer hyperparameters. In contrast, the more complex training requirements of Gaudi-GBLR, which involve smoothness annealing and proximal gradient descent, may lead to suboptimal results for a large model such as ViT-Base in Table 1 .We validate the training performance of BLAST weights on language models. We replace the weights of GPT-2 [20] with random BLAST 6 matrices and trained the network from scratch on the WikiText 103 [30] dataset for 100 epochs. In Figure 5 , we compare the test set perplexity of BLAST with the perplexity from low-rank, block-diagonal, Monarch, and Gaudi-GBLR matrices. Similar to the ImageNet training, we found that BLAST achieves the best perplexity-FLOPs trade-off. Compared to Gaudi-GBLR, BLAST obtains a significant perplexity gain. We attribute this improvement to the simple training process of BLAST which requires less hyperparameter tuning than that of Gaudi-GBLR.In this section, we discuss the performance of BLAST weights when pre-trained dense weights are available. We first compress the dense weights using Algorithm 2 and re-train the model on the training data with the cross-entropy loss. Diffusion Models We compress the weights of a Diffusion Transformer (DiT) [13] pre-trained on ImageNet using BLAST matrices and compare its performance to SVD-based low-rank approximation. For both techniques, we match the compression ratio such that both decrease the total number of model parameters by 50%, and re-train each model for 10 epochs on the ImageNet training set. We evaluate the models by generating a total of 50,000 images using the original, low-rank, and BLAST compressed models, and compute the FID [31] , sFID [32] , and IS [33] metrics with respect to the ImageNet validation set.The objective is to observe if the compressed model can generate images as realistic as the original uncompressed model.In Table 2 , we show quantitatively that the model compressed via BLAST significantly outperforms the model compressed via SVD. The low-rank compressed model often generates unrealistic images, leading to poor metrics such as the inception score. Figure 1 also contrasts how the BLAST matrices contribute to maintaining high perceptual quality as well as a close instance-wise resemblance with the uncompressed model outputs. Due to space limitations, we defer additional qualitative results and experimental setup to Appendix D.2.Large Language Models (LLMs) We compress the weights of Llama-7B [1] with BLAST matrices using Algorithm 2 by 20% and 50%, and re-train the models for 400 steps on a subset of SlimPajama [34] dataset using 0.49B tokens. The number of blocks b in the BLAST matrices is fixed at 16, and we use r = 1024 for the attention modules and r = 1488 for the MLP modules to achieve a 50% compression ratio. We test the WikiText-2 perplexity and the zero-shot task classification accuracy on common sense reasoning datasets including PIQA [35] , HellaSwag [36] , WinoGrande [37] , BoolQ [38] , OpenBookQA [39] , ARC-easy and challenge [40] . We report the performance of Low-Rank, Monarch, and Block-Diagonal weight matrices after compression at the same rate and re-training. In Table 3 , the first row presents the performance of the original Llama-7B model. On 50% compression ratio in the last five rows, the Monarch and Block-Diagonal matrices fail to recover the acceptable performance. Compared to Low-Rank weights, BLAST weights achieve the lowest performance degradation in WikiText-2 perplexity and zero-shot classification accuracy. The accuracy of each common sense reasoning benchmark and extended results can be found in Appendix D.3.We provide an analysis to quantify the performance impact of compression and re-training. We first quantify the weight compression performance at 20% compression ratio in significant performance loss, at most 4x higher perplexity and 25% accuracy degradation. On the other hand, BLAST 16 compression without re-training maintains reasonable accuracy and perplexity. This shows that the flexible and adaptive structure of BLAST matrices captures more information than other types of structured matrices. Yet, BLAST compression also exhibits noticeable performance degradation on a more intensive compression ratio (see the yellow curve in Figure 7 ). We provide more compression-only results on Diffusion Models and LLMs in Appendix D.We find that the re-training stage is crucial for converting a pre-trained model into an efficient version without losing significant accuracy when the compression ratio is high. In Figure 7 , we show the average zero-shot classification accuracy of the compressed models with 50% compression. The models with BLAST weights before re-training (yellow curve) exhibit substantial accuracy degradation at higher compression ratios. However, re-training (blue curve) effectively recovers performance using only 0.49B tokens and 400 steps.We evaluate the runtime of the Llama-7B compressed by the BLAST matrices on the text generation task. For evaluation, we let the model generate the sequences of length L =10, 100, and 1000 ten times and report the average and the standard deviation of model runtime in Table 4 . The instruction we use to generate the desired sequence length is "Increasing sequence: one," and the model generates the text sequence such as "two, three, " and so on, with a batch size of 1. All runtime evaluation tests were conducted on a single 40GB NVIDIA A100 GPU after compiling the script using torch.compile(). The 20% compressed model shows a 12%∼15% runtime reduction without any library function customization for BLAST matrix multiplication. The speedup when b = 2 is higher than when b = 16 because a larger number of blocks increases the computation overhead to perform Equation (3). Notably, the 50% compression ratio provides 32%∼35% runtime reduction when b = 16. We note that the test is highly memory-bandwidthbounded. Thus the speedup reported in Table 4 can be mostly attributed to the reduction of parameters (i.e., memory accesses) rather than the reduction in FLOPs due to BLAST compression. 
related works: 
 Structured Matrix with Shared Bases Sharing the bases of block-structured matrices has recently drawn interest due to its considerable memory savings. BLAST matrices exemplify this approach.Ashcraft et al. [41] extend the BLR [15] format to BLR 2 , incorporating shared bases and block-wise low-rank coupling factors determined through LQ or QR factorization. Similarly, Yen et al. [42] apply the concept of shared bases in the Block-Diagonal preconditioner for DNN weights. While BLAST also shares the bases of blocks, it is distinct in having a diagonal coupling matrix, as shown in Equation ( 2). The design of BLAST matrices aims to enhance efficiency and learn a variety of structures, from low-rank to high-rank block matrices. More importantly, identifying the diagonal coupling factors in BLAST matrices does not necessitate QR decomposition. Instead, they can be updated via gradient descent, making this approach well-suited for modeling the weight matrices in deep learning models.DNN Weight Pruning and Decomposition Earlier work on DNN pruning [43] [44] [45] identifies less important parameters from the Hessian or magnitude to sparsify the weight matrix. Unlike general sparse matrices, a group sparse matrix skips computation in a group-wise manner by pruning channels [46, 47] . Sparse GPT [11] successfully prunes large language models with 50% sparsity without significantly degrading the performance of the original model. The model can achieve actual speedup utilizing 2:4 sparsity [48] on specific GPUs. However, 2:4 sparsity requires accelerators with specific architectures (e.g., NVIDIA A100 GPUs) and supports only the 50% compression ratio. On the other hand, BLAST is device-agnostic since it can be implemented with basic matrix arithmetic operations and offers diverse compression ratios.Low-rank matrices have been widely adopted for CNN compression [49, 50] and Transformer compression [51, 52] . Additionally, Butterfly [53] and Monarch [14] Learning Low-dimensional Structures of DNN Weights Similar to the BLAST matrix, a Gaudi-GBLR matrix [12] enables learning low-dimensional structures by gradient descent in the generalized structured matrix space. Gaudi-GBLR overlaps a variable number of zero-padded rank-1 blocks to model high-rank submatrices. Although Gaudi-GBLR can express a wider spectrum of matrices than BLAST, the matrix-vector multiplication for Gaudi-GBLR is less efficient because GPUs and typical neural processors cannot handle zero-padded vectors efficiently. In contrast, the BLAST matrix-vector operation does not involve zero padding, allowing for more efficient execution in hardware for the same FLOPs (as shown in Figure 4 , Figure 6 , and Table 1 ). 
conclusion and future work: 
 In this work, we introduced the BLAST matrix designed to improve the inference efficiency of large DNNs. The BLAST matrix represents various low-dimensional structures of the weight matrices with fewer parameters, while enabling efficient matrix-vector products. The BLAST factors are either learnable from data or estimated from existing weights using our preconditioned factorization algorithm. Our results on both language and vision tasks highlight the effectiveness of BLAST.Limitations and Future Work The BLAST matrix-vector product consists of three steps, as detailed in Equation ( 3), which may degrade hardware-execution parallelism. In our evaluation, we used the same computational budget r for all matrices. Learning an adaptive budget per layer or matrix (e.g., via overparameterization [54, 7] ) could further improve BLAST performance, which is left for future work. The proposed method has not been evaluated on tiny (<100M parameters) or extremely large (>10B parameters) DNNs. Additionally, optimizing runtime and power consumption via BLAST matrices with customized library functions and/or hardware accelerators also remains as future work. Furthermore, a deeper theoretical investigation into the behaviors of BLAST matrices would provide a more comprehensive understanding of their capabilities and limitations. Applying advanced re-training techniques, such as knowledge distillation [55] or iterative compression and distillation [56] , to the BLAST compression pipeline is also left for future work. Finally, beyond the weight structures, we expect BLAST can also help understand and exploit low-dimensional data manifolds [57] [58] [59] in future work. 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: Instead of providing a decomposed Llama-7B model, we provide the code file to reproduce the result if one can access the original safeguarded model.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [Yes] Justification: Appendix C.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. 
 
 Appendix: Matrix Multiplication Code Implementation Following the conventional setting in Transformers [16] , the input tensor is assumed to have batch, sequence, and channel dimensions. The left and right factors are multiplied using the batched matrix multiplication routine, whereas the diagonal factors are multiplied via broadcasting and summation. Block diagonal matrix A block-diagonal matrix is a BLAST matrix when r = p and s i,j = 1 r if i = j 0 r otherwise since    A 1,1 A 2,2 . . . A b,b     =      U 1 diag(s 1,1 )V T 1 U 2 diag(s 2,2 )V T 2 . . . U b diag(s b,b )V T b      .When r < p, BLAST matrices model block-diagonal matrices which have low-rank diagonal blocks.Block low-rank (BLR) matrix For ease of understanding, let us consider a BLR matrix of 9 n 3 × n 3 rank-1 blocks. Each block is composed of the unique bases A i,j = u i,j v T i,j . Now consider U i = [u i,1 , u i,2 , u i,3 ] and V j = [v 1,j , v 2,j , v 3,j ]. Then, the BLR matrix is a BLAST matrix with r = b = 3:   u 1,1 v T 1,1 u 1,2 v T 1,2 u 1,3 v T 1,3 u 2,1 v T 2,1 u 2,2 v T 2,2 u 2,3 v T 2,3 u 3,1 v T 3,1 u 3,2 v T 3,2 u 3,3 v T 3,3   =   U 1 S 1 V T 1 U 1 S 2 V T 2 U 1 S 3 V T b U 2 S 1 V T 1 U 2 S 2 V T 2 U 2 S 3 V T 3 U 3 S 1 V T 1 U 3 S 2 V T 2 U 3 S 3 V TU i , V j ∈ R p×(bt) , s i,j ∈ R bt .By letting s i,j,k = 1 if t(j -1) + 1 ≤ k < tj + 1 0 otherwise , the BLAST matrix can model the BLR matrix.Note that the number of parameters of the BLAST matrix is 2nr + rb 2 , whereas that of the BLR matrix in this case is b 2 • (p + p)t = 2(pb)(bt) = 2nr. In other words, the BLAST matrix models various matrices with the cost of rb 2 .We rewrite the loss function in Equation ( 4) below:EQUATIONWe first derive the gradients of U i , V j , and s i,j , then discuss the preconditioning matrix for each factor.Here we derive the gradients of Equation ( 4) with respect to the BLAST factors. We begin with introducing the short-handed notation for the concatenated factors:V T i = [S i,1 V T 1 • • • S i,b V T b ], Ūj = [(U 1 S 1,j ) T • • • (U b S b,j ) T ] T .That is, the matrix Vi is composed by concatenating V T j s horizontally along j = 1, 2, . . . , b after scaling them with S i,j . Ūj is defined similarly by concatenating the scaled U i s vertically.Gradient of U i We only have to consider the loss term related to U i . Therefore, we have the following gradient expression:EQUATIONwhere for the second equality we used the concatenated version of the first line.Gradient of V j follows the similar derivation as Equation (10) :EQUATIONGradient of s i,j We consider the block-wise loss for the gradient:EQUATIONSince the Frobenius norm can be expressed by a matrix trace, the loss is written as follows:A i,j -U i diag(s i,j )V T j 2 F = Tr A i,j -U i S i,j V T j T A i,j -U i S i,j V T j = Tr V j S i,j U T i U i S i,j V T j -2A T i,j U i S i,j V T j + A T i,j A i,j = Tr S i,j V T j V j S i,j U T i U i -2S i,j V T j A T i,j U i + A T i,j A i,j, where Tr(X) is the trace of X. Note that the derivative of product in trace is given by ∇ X Tr(XY ) = Y T for any two conformal matrices X and Y . Therefore, we have∇ Si,j 1 2 A i,j -U i diag(s i,j )V T j 2 F = U T i U i S i,j V T j V j -U T i A i,j V j .Now Equation ( 12) becomes as follows:EQUATIONThe first term on the right hand side is further arranged by using the fact that diag XY T = (X ⊙ Y ) 1 for any two matrices X, Y of the same size.EQUATIONHence, the gradient of s i,j is now expressed as follows:EQUATIONNow let us derive the preconditioning matrices used in Algorithm 2.Preconditioning matrix P Ui for U i Let V j and s i,j are given. Let us consider the case when U i is at the stationary point Û which satisfies∇ Û 1 2 ∥A i, * -Û V T i ∥ 2 F = Û V T i Vi -A i, * Vi = O,where O is the zero matrix. This gives us the normal equationEQUATIONNow consider a preconditioned gradient descent with P Ui ∈ R r×r :U ′ i = U i -U i V T i Vi -A i, * Vi P Ui = U i -U i V T i Vi -Û V T i Vi P Ui = U i -U i -Û V T i Vi P Ui =⇒ U ′ i -Ûi = U i -Û I -V T i Vi P Ui Suppose V T i Vi is invertible.Then the ideal preconditioner isEQUATIONsince it brings U ′ i to the stationary point Û . However, directly use the inverse of V T i V might result in numerical instability or complete breakdown of the algorithm when the matrix is singular. By following [24] , we use the regularized versionEQUATIONwhere δ is chosen byEQUATIONPreconditioning matrix P Vj for V j The preconditioner for V j can be derived by following the similar steps for P Ui . Here we present the result:EQUATIONwhere δ is chosen by Equation (19) .Preconditioning matrix P si,j for s i,j We again consider the stationary point ŝ or the diagonal version Ŝ = diag(ŝ) which has a zero gradient:∇ ŝ 1 2 ∥A i,j -U i ŜV T j ∥ 2 F = (U T i U i ) ⊙ (V T j V j ) ŝ -diag U T i A i,j V j = 0 =⇒ (U T i U i ) ⊙ (V T j V j ) ŝ = diag U T i A i,j V j .Therefore, Equation ( 15) can be written as follows:∇ si,j ℓ(U * , V * , s * , * ) = (U T i U i ) ⊙ (V T j V j ) (s i,j -ŝ). Now consider the preconditioned gradient descent:s ′ i,j = s i,j -P si,j (U T i U i ) ⊙ (V T j V j ) (s i,j -ŝ) =⇒ s ′ i,j -ŝ = s i,j -ŝ -P si,j (U T i U i ) ⊙ (V T j V j ) (s i,j -ŝ) = I -P si,j (U T i U i ) ⊙ (V T j V j ) (s i,j -ŝ)The Hadamard product of two positive definite matrices are still positive definite (see [60, Theorem 7.5.3], also known as Schur's Product Theorem). Hence, if both U T i U i and V T j V j are positive definite so that invertible,(U T i U i ) ⊙ (V T j V j) is also invertible. The ideal preconditioner when both matrices are invertible is thereforeEQUATIONSame as for P Ui and P Vj , we also consider the regularized version:EQUATIONHere, δ is chosen by Equation (19) .We first introduce the following properties: Lemma 2. [61, Lemma 3.4 ] Assume f : R n → R is convex and continuously differentiable, and its gradient is L-Lipschitz continuous. Then for any x, y ∈ R n , one hasf (y) -f (x) -⟨∇f (x), y -x⟩ ≤ L 2 ∥y -x∥ 2 2 .Proof. See [61, Lemma 3.4 ].Lemma 3. Assume f : R n → R is convex and continuously differentiable, and its gradient is L-Lipschitz continuous. Consider a gradient descent updatex (k+1) = x (k) -η • ∇f (x (k) ).Then, with the step size 0 < η ≤ 1 L , the following holds:f (x (k+1) ) ≤ f (x (k) ) - 1 2L ∥∇f (x (k) )∥ 2 2 .That is, the gradient descent update does not increase the function value.Proof. By Lemma 2,f (x (k+1) ) ≤ f (x (k) ) + ⟨∇f (x (k) ), x (k+1) -x (k) ⟩ + L 2 ∥x (k+1) -x (k) ∥ 2 2 = f (x (k) ) -η • ∥∇f (x (k) )∥ 2 2 + η 2 L 2 ∥∇f (x (k) )∥ 2 2 = f (x (k) ) -η • 1 - ηL 2 ∥∇f (x (k) )∥ 2 2 ≤ f (x (k) ) - η 2 ∥∇f (x (k) )∥ 2 2 (since ηL ≤ 1) ≤ f (x (k) ) - 1 2L ∥∇f (x (k) )∥ 2 2 .Now we prove Theorem 1, which we restate below.Theorem 1. Let A i,j ∈ R p×p be a target block andEQUATION∈ R p×r , and s (k) i,j ∈ R r be factors of a block in the BLAST matrix to be optimized. With the step sizes 0< η U (k) i ≤ 1/σ 1 V (k)T i V (k) i , 0 < η V (k) j ≤ 1/σ 1 Ū (k)T j Ū (k) j , 0 < η s (k) i,j ≤ 1/σ 1 (U (k+1)T i U (k+1) i ) ⊙ (V (k+1)T j V (k+1) j) , the gradient descent updates in Equations ( 5) to (7) monotonically non-increase the loss:ℓ(U (k+1) * , V (k+1) * , s (k+1) * , * ) ≤ ℓ(U (k) * , V (k) * , s (k) * , * ).Proof. To prove Theorem 1, we first show that each step of Equations ( 5) to (7) satisfies Lemma 3 under the given conditions. Then we resemble the results to construct the bound.Gradient Descent Update on U i Let us denote the loss term regarding U i byℓ(U i ) = 1 2 A i, * -U i V T i 2 F .Since (i) the Frobenius norm is convex, (ii) all linear mappings are convex, and (iii) a composition of two convex functions are convex,A i, * -U i V T i 2F is a convex function of U i . Also, the gradient we derived in Equation ( 10) always exists and has the following vectorized form:vec(∇ℓ(U i )) = vec(U i V T i Vi ) -vec(A i, * Vi ) = (( V T i Vi ) T ⊗ I)u i -vec(A i, * Vi ), where ⊗ denotes the Kronecker product and u i = vec(U i ). The Lipschitz constant of the gradient is the largest singular value of the matrix ( V T i Vi ) T ⊗ I, which is the largest singular value of V T i Vi since the Kronecker product of two matrices of singular values Σ 1 and Σ 2 has the singular values of Σ 1 ⊗ Σ 2 (see [62, Theorem 4.2.15] ). Therefore, from Lemma 3, we obtain the following bound:EQUATIONGradient Descent Update on V j The loss function ℓ(V j ) = 1 2 ∥A * ,j -Ūj Ū T j ∥ 2 F with respect to V j is convex to V j and the gradient of V j in Equation ( 11) also always exists and can be rewritten as follows:vec(∇ℓ(V j )) = vec(V j Ū T j Ūj ) -vec(A T * ,j Ūj ) = (( Ū T j Ūj ) ⊗ I)v j -vec(A T * ,j Ūj ).The Lipschitz constant of the gradient is again the largest singular value of Ū T j Ūj . We have the bound from Lemma 3 similar to Equation ( 23):EQUATIONGradient Descent Update on s i,j The loss functionℓ(s i,j ) = 1 2 A i,j -U i diag(s i,j )V T j 2 Fis also convex in s i,j since diag(•) is a convex mapping. We know the gradient exists from Equation (15): 8 : Hyperparameters used for the BLAST 9 -compressed DiT-XL/2 with 20% compression ratio. m, n: size of the original matrix, b: number of row/column partitions, r: BLAST rank parameter, Layer Indices: indices of layers that the BLAST matrix replaces the weight. δ 0 = 0.1. We linearly decayed the step size from 1.0 to 0.0. The factorization process takes 3.38 GPU hours for the BLAST weights of b = 16 on NVIDIA A40 GPUs.∇ℓ(s i,j ) = diag U T i U i diag(s i,j )V T j V j -diag U T i A i,j V j = U T i U i ⊙ V T j V j s i,j -diag U T i A i,j V j , m n bTo re-train the compressed Llama models, we used a subset2 of the SlimPajama dataset [34] for 400 steps using 0.47B tokens. The global batch size was set to 576, and the models were trained on 4 NVIDIA L40S GPUs. See Table 6 for details.We used Language Model Evaluation Harness3 [67] for the zero-shot classification accuracy evaluation. In this experiment, we test the factorization algorithms discussed in Section 3, similar to Figure 3 but with a different target matrix. To be specific, we synthesize a 256 × 256-sized BLAST 16 (i.e., b = 16) target matrix with r * = 8. Then, we compare the error curve along the iterates of the gradient descent without preconditioning (GD) in Equations ( 5) to (7) , and the preconditioned gradient descent (PrecGD) in Algorithm 2. Here, we consider the exact parameterization setting when r = r * = 8 and the over-parameterized setting r = 32 > r * . In Figure 9 , unlike the low-rank target matrix, GD does not converge in both cases. However, the preconditioned version easily finds the low-error solution for the exact parameterization. For the overparameterized case, the preconditioned gradient descent method in Algorithm 2 achieved in two orders of magnitude improvement from the simple GD. We include extended experimental results of the diffusion model compression in Section 4.2Compression-only In Figures 10 and 11 , the additional image samples of original uncompressed, low-rank-compressed, and BLAST 9 -compressed DiT [13] models are presented. The images in the same column were sampled using 250 DDIM steps, starting from the same noise vector. The compression ratio was set to 20% for both models. The figures show that the outputs of the model compressed by BLAST maintain similar features and perceptual quality to the outputs of the original DiT.We present additional samples from the low-rank and BLAST DiT models at the 50% compression ratio after re-training in Figure 13 . Similar to Figure 1 , the images generated by the low-rank DiT lose significant image quality, whereas the images from the BLAST DiT preserve the quality and semantics of the original samples. Overall, the low-rank approximated model often generates unrealistic images, which contributes to low scores evident in Table 2 .Evidence of low performance of low-rank-compressed model Some images generated by the 20% low-rank-compressed model (Figure 12 -left) are highly unrealistic and have inferior quality compared to the images generated by original and BLAST-compressed models (Figure 12 -right). We observe that these samples contribute to the low scores in Table 2 . These samples were computed using 250 DDPM steps, as done by the original DiT work [13] .Compression-only We report the performance of LLM-Pruner [47] and Joint Rank-k [66] for the same compression task. LLM-Pruner [47] identifies sparse weights with data to pinpoint unimportant neurons. Joint Rank-k [66] performs a low-rank approximation jointly on weight matrices with similar column spaces by stacking them and applying a truncated SVD.In Table 12 , we compare the performance degradation of LLM-Pruner, Joint Rank-k, Low-Rank, BLAST 2 , and BLAST 16 , as well as their absolute performance. The first four rows represent the zero-shot performance of Llama-7B [1] from the literature, while the row marked with an asterisk ( * ) indicates our results.For 10% compression, Joint Rank-k achieved the lowest performance degradation, although BLAST 16 also exhibited a similar performance drop. When the model is compressed by 20%, BLAST 16 surpasses Joint Rank-k [66] , LLM-Pruner [47] , and low-rank schemes. The zero-shot accuracy versus compression ratio curve in Figure 7 shows that BLAST compression results in less performance drop for the same compression ratio. Our proposed method targets improving the efficiency of the DNN inference. This might have a negative social impact by promoting the accessibility and usability of malicious DNNs such as DeepFake. However, at the same time, we expect the BLAST matrix will bring a tremendous positive social impact. First, it contributes to sustainability by cutting down the energy consumption for the DNN inference. Moreover, the BLAST matrix can improve the accessibility of AI-based medical, educational, and social services by providing a foundation for running those models on mobile devices. Therefore, we believe BLAST will give tangible benefits to our society.• The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML).For an m × n rectangular matrix, we partition m rows into b chunks assuming that b divides m as well.https://huggingface.co/datasets/DKYoon/SlimPajama-6Bhttps://github.com/EleutherAI/lm-evaluation-harness 


Decision: Overall: this paper describes an interesting variation on low-rank/structured weight matrices with enough experimentation to suggest that the structure might be worth further investigation, and hence sharing with the NeurIPS audience.  The reviewers, particularly xrxa, draw attention to many deficiencies of experimentation and presentation, which it is strongly expected the authors will take into consideration when revising the paper.

On balance, however, I feel that the objections of reviewers V5gm and 2dKz are adequately rebutted, which puts this paper into the "borderline accept" category.  This is largely due to xrxa's observation that the "proposed approach seems to be entirely novel, building on existing low rank and block sparse techniques."  It feels as if it is valuable to have this approach described, even without the custom kernels, and  that would make its impact more widely 

As AC, I agree with xrxa: "the speedups should be rewarded, whether they are due to FLOPs or memory savings."  
However, this discussion does become relevant when rebutting the observation that a very short sequence length used, which artificially amplifies the performance gains by loading much more work onto the linear layers.  The rebuttal claims that "we observed the FLOPs saving is relevant to the runtime reduction", but this is more contentious given xrxa's observations.

Reviews: review: paper_summary: paper_summary: This paper attempts to find efficient structures in weight matrices of deep learning models. The basic idea is to learn a group of block-wise low-rank matrices via gradient descent. The proposed method replaces original dense weight matrices and hence needs retraining.
strength_weakness: strength_weakness: * The significance of the proposed method remains unclear. Compared with prior methods such as Gaudi-GBLR [14], the accuracy results of this work do not show consistent superior performance nor with a convincing explanation.
* The method does not seem to improve efficiency much, showing limited significance.
questions: questions: * How different is this work compared with block sparse + low-rank [23] and Gaudi-GBLR [24]? Why is adaptive and learnable a new contribution when these ideas have already been proposed? The paper [23] is not introduced before use in experimental results.
* In Table 4, when L is 100, the execution time of 20% compression ratio increases compared with no compression, while other two columns show slight reduction of execution time. What is the reason? 
* Please provide more details, even if it is obvious, that the learned BLAST matrix structures do not need any library function customization. What is specifically used as the backend compute kernel?
limitations: limitations: Yes.
paper_summary: paper_summary: The authors propose a learnable compressed representation (BLAST) for weight matrices used in deep learning which enables lower complexity matrix multiplications which approximate the full, uncompressed operation.  A BLAST matrix decomposes the original matrix into a grid of blocks of diagonal matrices with shared low-rank factors along each (block) row and column, and these parameters can be learned directly from a random initialization.  Converting a dense matrix to a BLAST matrix by minimizing the Frobenius norm error between original and compressed matrices can be done through gradient descent factorization, and preconditioning can accelerate convergence.  The proposed technique is evaluated in pre-training, fine-tuning, and data-free compression scenarios, using a variety of networks and tasks (vision transformers and LLMs, image classification, image generation through diffusion, and zero-shot language modeling).
strength_weakness: strength_weakness: **Originality**: While the authors included a sparse-weight baseline (LLM-Pruner), [SparseGPT](https://arxiv.org/abs/2301.00774) deserves mention (and potentially inclusion as a baseline method to compare against).

**Quality**

The authors claim in line 310 that BLAST matrices allow "more efficient execution in hardware for the same FLOPs (in Figure 5, Figure 6, and Table 1)," but these figures and table don't tell the reader anything about hardware efficiency.  While the authors have presented an algorithm for matrix-vector products, I'm not convinced that matrix-matrix multiplication is as efficient as optimized dense routines tensor execution units in any recent hardware.  (The runtime speedups in Table 4 may simply be due to the memory and bandwidth savings, not an computational gains.)

The authors stress that U and V are shared along rows and columns, but it's unclear why this is so important.  What if each block has its own U and V -- other than increasing the number of parameters, would the factorization process, quality of the results, multiplication algorithms, etc. change?

In line 245, the authors suggest that certain techniques "have not been studied as data-free compression techniques" and are therefore "excluded from the comparison" for this task.  I'm not sure if this is fair or not: were these techniques simply not designed for this task, or are they harder to adapt to this task, or are they entirely incompatible?  

**Clarity**: there are some configuration details missing:

- The rank *r* used for the experimental results.  (I see this detail is in the appendix, but I think it's useful to give at least a range of the values used in practice, as was done with *b*, in the main content.)

- The batch size and input sequence length used in the LLM runtime analysis.


Figure 4 is lacking vertical axis labels.

I'm not sure of the difference between, in line 100, "learning" and "finding" the structure of a weight matrix.


**Significance**: the quality of the results lets the submission down, in a way.  Compression rates of 10%-20% are fairly limited, and I'd be curious to know what happens, even if it's catastrophic, at higher compression rates -- even 50% would be informative.  Further, the quality of the BLAST diffusion model's output, while better than SVD's, still fall short of the baseline.  Fine detail is lost (in e.g. the water and rocks of the fourth column, the feathers in the fifth column, and the jersey in the sixth), and larger distortions appear (the hippo's snout no longer looks accurate, the bird's eye has disappeared, and the hot air balloon's shape is extremely lopsided).  Again, these results are better than the other Low-Rank technique, but practitioners may be hesitant to use the technique for this task without significant gains.
questions: questions: 1. Is matrix-matrix multiplication built by performing multiple matrix-vector products (Algorithm 1), or is there something more clever that can be done?

2. What is the difference, in line 100, between "learning" and "finding" the structure of a weight matrix?

3. What is the benefit (other than parameter reduction) in sharing the U and V factors along rows and columns?  Can it be quantified?

4. What rank *r* is used for the experimental results?

5. Regarding skipping some baselines for comparisons, were these techniques simply not designed for this task, or are they harder to adapt to this task, or are they entirely incompatible?  

6. What batch size and input sequence length was used for the results in Table 4?
limitations: limitations: The authors' discussion is excellent, only missing one potential further limitation: this technique seems like it'd need extra work to apply to higher-dimensional tensors.  I could imagine applying it on matrix-slices of a tensor, but this may not work well.
paper_summary: paper_summary: This paper proposes to compress linear layers by using a block low-rank structure with shared basis -- also known as the BLR^2 matrix structure. The basis U, V, and diagonal weights s of the low-rank block are computed through a gradient-based optimization that minimizes the Frobenius norm between the dense and low-rank matrix. They also use preconditioning to accelerate the convergence.
strength_weakness: strength_weakness: What is proposed as the BLAST matrix in this paper is an existing matrix structure known as BLR^2 [https://hal.science/hal-03070416]. The original BLR^2 structure uses LQ decomposition to compress the low-rank blocks, but the name of the structure is independent of the method used to compress the low-rank blocks. Therefore, the proposed method which uses a gradient based iterative method to compress the low-tank blocks is still a BLR^2 structure. The BLR^2 paper should be cited, and all mention of BLAST should be replaced with BLR^2. 

For CIFAR-10 and CIFAR-100, Gaudi-GBLR seems to perform better than the proposed method. The authors mention that Gaudi-GBLR's "capability of learning the adaptive resource/budget allocation for each weight matrix" is the cause, but it is unclear why this is not the case for ImageNet-1k.
questions: questions: Why does the paper repeatedly mention matrix-vector multiplications? The input to DNNs are not vectors.

Why do the authors mention that "A block low-rank (BLR) [12], also known as Monarch [13] matrix"? The Monarch matrix is a very different structure from BLR. 

Why are Gaudi-GBLR, BSP+LR, and BLR excluded from the compression benchmark in Table 2? Just because they haven't been studied as data-free compression techniques doesn't mean they are not suitable for this task.
limitations: limitations: The limitations are stated clearly in a separate section.

Comments: title: Thank you for the excellent discussion
BLAST GEMM Runtime Evaluation
SparseGPT can target structure, Batched GEMMs are inefficient, image generation results are still unsatisfying
comment: Q1: We appreciate the reviewer's feedback regarding potential confusion in naming. However, it's important to note that BLAST matrices are significantly distinct from BLR and BLR^2 matrices regarding the properties of the factored matrices. Moreover, the set of BLAST matrices are not a subset of the set of BLR^2 matrices, vice versa. We believe that using the same name as BLR^2 would lead to greater confusion as these are two distinct matrix structures. 

Q5: Despite our best efforts, we were unable to identify candidate algorithms capable of compressing the dense matrix using BSP+LR or GBLR matrices. This limitation prevented us from including these structures in the data-free compression experiment baselines. We would be extremely grateful if the reviewer could suggest any potential compression methods that we might have overlooked.

We appreciate Reviewer 2dKz's effort and interest in our work.
Thank you for the insightful and constructive feedback! We acknowledge the reviewer's concern regarding the mechanism behind the speed-ups and appreciate the valuable suggestions. We also believe the source of latency improvement should not be over-specified. We will update the paper accordingly, along with some information on the achievable speed-ups for the hyperparameters. Once again, we appreciate Reviewer xrxa's effort and interest in our work.
Q1. I acknowledge the difference between BLR^2 and BLAST. I still believe that the terminology to describe the matrix structure should be decoupled from the methods used to compress them. For example, the blocks of the BLR and BLR^2 matrices can be compressed by using methods such as SVD, randomized SVD, rank revealing QR, adaptive cross approximation, interpolative decomposition, or learning. Some of these compression schemes do not result in orthogonal bases or diagonal S blocks, but the matrix structure is still called BLR and BLR^2 in previous studies. If everyone renamed the matrix structure every time they changed the compression method, it would be very confusing.

Q5. I don't understand the logic of excluding Gaudi-GBLR and BSP+LR from the comparisons, just because they don't have dedicated compression algorithms. It should be trivial to simply use the best compression algorithm for those methods and compare with them.
**Speedups**

I think this discussion is incredibly useful, and I thank the authors for continuing to explore this area to make sure the benefits are clear.  However, I do not believe their claims have been proven.

Let me take this argument in particular: "if both 2:4 sparsity and BLAST have the same compression rate, then any speedup of BLAST is due to FLOPs savings."  This isn't the case, simply because the relative efficiencies of 2:4 sparse and dense GEMMs are not the same, much less 2:4 sparse and BLAST GEMMs.  Due to the complex nature of parallelizing a GEMM among O(100) processing elements, achieved efficiency (math utilization) can change drastically with small changes in workload size.  For example, while the dense GEMM size in your new table scales roughly linearly with N^3, the same is not true of the 2:4 sparse row; its speedups fluctuate non-monotonically from 1.35x to 1.54x.  From this, I conclude that the speedup of a 2:4 sparse GEMM is a function of more than just FLOPs performed and memory required for the compressed representation: some other factor (or factors) affects the speedup.  The same is clearly true of BLAST matrices: your BLAST4 row's speedups change (monotonically, in this case) from 1.09x to 1.57x.  So, with 2:4's efficiency changing, one can't compare to it directly and suggest that any speedup over *its* results are due to a savings in FLOPs and not the other complicating factors shown to be at play in the 2:4 results.

Similarly, for Q.10, the very short input length means most of the time is spent in the autoregressive generative phase, which uses an effective sequence length of one when generating new tokens.  With a batch size of one (from your original response's Q.10), this makes the effective GEMM size not, for example, 4096x4096 * 4096x4096, but rather 4096x4096 * 4096x1, which will have very different characteristics than the workloads in your most recent table, so it's not reasonable to equate the two.  I used your pseudocode to implement your benchmark and came up with similar results to those shown in your table, but things looked very different for BS=1, seqlen=1.

… all this is to say that understanding speedups is nontrivial, and it is not simply a function of FLOPs and memory.  So, what does that mean for this submission?  I think there are two things to keep in mind:
1) provide the readers with some idea of the achievable speedups for the range of BLAST hyperparameters used in practice (b=(2..16), r=chosen to target compression rates of 10% to 75% (I think, based on my understanding of Table 8 in Appendix C)).  This could be as simple as extending your table to include b=8 and 16, or as complete as also providing another similar table with sequence length = 1 (for autoregressive generation) and batch size sweeping from 1 to N (even for a single N).
2) avoid attributing some throughput or latency behavior to FLOPs or memory without being very confident in the source of the benefit.  It's safer and less potentially misleading to simply share the throughput and latency behavior in these cases.  Concretely, for Line 296, I might suggest "… a larger number of blocks reduces the efficiency of the batched GEMM operations" - it's accurate, easy to show, and doesn't over-specify the source of the inefficiency.

I am very grateful to the authors for sharing their findings here, and I want to make it clear that I am not weighing my disagreement with their analysis too heavily!  After all, as they say, the speedups should be rewarded, whether they are due to FLOPs or memory savings.

For completeness and inspection, the runnable PyTorch code I used when benchmarking:
```
def blast_matmul(
        X,  # input,        shape=(B, n, q*b), B=batch_size, n=num_seq
        U,  # left factor,  shape=(b, p, r), b=num_blocks, r=rank
        S,  # diag factor,  shape=(b, b, r)
        Vt, # right factor, shape=(b, r, q)
    ): # output = X @ A.T where A is a BLAST matrix of shape (b*p, b*q).
    B = X.shape[0]
    n = X.shape[1]
    b = U.shape[0]
    r = U.shape[2]
    q = Vt.shape[2]
    p = U.shape[1]

    X = X.view(b, B*n, q) #rearrange(b, B*n, q)
    Y = torch.bmm(X, Vt.mT) # multiply right factor
    Z = Y.unsqueeze(0) * S.unsqueeze(2) # multiply diagonal factor
    Z = Z.sum(1) # aggregate, shape=(b, B*n, r)
    Out = torch.bmm(Z, U.mT) # multiply left factor
    return Out.view(B, n, b*p)
```

And my calculation to choose r to hit a target compression rate:
```
r = comp_rate*n / (2 + (b*b/n))  # find the exact r to achieve comp_rate (0.0 .. 1.0) memory savings
r = int((r//16) * 16)  # set to the next-lowest multiple of 16 for hw-friendliness (next-lowest results in more compression than strictly necessary)
```
Thank you for your explanation.
**Q1. SparseGPT can target structure.**

Thank you for the comments. We completely agree that SparseGPT can be accelerated on the Ampere architecture (or higher version of) NVIDIA GPUs via 2:4 sparsity, although it is our understanding that the 50% compression ratio on LLMs is very challenging to achieve with the structured methods in general (e.g., Structured Pruning \[30\]).   
Nonetheless, we would like to clarify that BLAST is less device-specific – one can use BLAST to accelerate the DNN on the legacy NVIDIA GPUs as well as the devices from other vendors. We believe this property will play a crucial role for accessibility of the proposed method.   
   
Moreover, 2:4 sparsity can be adopted to further compress the BLAST factors or vice versa, which is an interesting combination of two different worlds. We are pleased to work on fusing two methods as future work.

**Q2. Speedup of BLAST is due to only the memory reduction, not the FLOPs savings, and custom kernels may not be able to change this behavior.** 

Thank you for the insightful feedback and experiments\! We would like to provide additional results on the BLAST GEMM, inspired by your feedback. The results show that **BLAST GEMM can benefit from FLOPs savings.** 

We benchmarked (N\*N x N\*N) dense, 2:4 sparsity, and 50%-compressed BLAST matrix multiplications. The structure of another matrix is fixed to dense. In the table below, we report the cuda time of 1000 iterations using torch.profile on an NVIDIA A100 GPU with 40GB memory.

To be completely fair with the memory bandwidth, we matched the number of parameters of BLAST to that of 2:4 sparsity, i.e., 17/32 of N^2. In this manner, both 2:4 sparsity and BLAST have the same number of values to be transferred from VRAM to CUDA blocks. In other words, the memory-related latency is similar in both cases. Hence, *if the computation reduction does not contribute to the inference speed-up of BLAST, the runtime of BLAST GEMM will be similar to or slower than the runtime of 2:4 sparsity GEMM.*

However, BLAST GEMM is **faster** than Semi-structured Matmul when N\>2048 and b=2 (\# blocks=2x2) or N\>4096 and b=4. Based on the fact that their latency related to memory bandwidth is similar, the results indicate that the extra speedup compared to 2:4 sparsity comes from the *FLOPs saving*. 

For reference, we included the runtime of dense GEMM with 50% FLOPs (N\*N x N\*0.5N). The majority of BLAST GEMM runtimes are in between the runtime of 50% dense and 2:4 sparsity, which also shows that the latency improvement of BLAST GEMM is relevant to the FLOPs saving.

   
**Inefficiency of Batched Matmul.** We agree that BMM is less efficient than the dense GEMM of the same FLOPs when the batch size is large. This hinders the efficient CUDA core utilization when the number of blocks is large and each block is small. We will work on circumventing this issue for future work. 

Again, thank you so much for the constructive feedback and discussion. We will discuss this in the revised version. In particular, we will modify Line 296 to “...a larger number of blocks increases the computation and the **memory bandwidth** overhead to perform Equation (5).”  
   
Table.

| Models | N=2048 | N=4096 | N=6144 | N=8192 |
| :---- | ----- | ----- | ----- | ----- |
| Dense | 0.339s | 2.556s | 8.198s | 19.325s |
| Dense (50% FLOPs) | 0.209s | 1.313s | 5.098s | 9.861s |
| 2:4 Sparsity | 0.252s | 1.660s | 6.424s | 14.744s |
| BLAST (b=4) | 0.311s | 1.938s | 5.638s | 12.350s |
| BLAST (b=2) | 0.282s | 1.581s | 4.680s | 11.602s |

**Q10. Very short input length.** We agree that BLAST takes advantage of the short input sequence in Table 4\. However, based on our runtime analysis in Q2, we observed the FLOPs saving is relevant to the runtime reduction.

**Q6. Image quality degradation.** Thank you for the clarification. We believe there is room to improve the image quality by carefully choosing the layers for compression. But we leave this for future work.
I appreciate the authors' thoughtful responses, but I will keep my rating. 

Q1. SparseGPT also generalizes easily (as evidenced by the results in the paper and public code) to semi-structured (or N:M) sparsity, which is supported in the hardware and software you used for evaluation (NVIDIA A100 and PyTorch [1]).  They did not experiment on Llama-7B, as you did for the "aggressive compression" results in the new Figure B (thank you for including this!), but for OPT-175B (admittedly, a much larger model), 2:4 sparsity (50% sparsity with a hardware-compatible format) sees only a 1.2 percentage point drop in average zero-shot accuracy, as compared to BLAST's ~15-20 percentage point drop (it's hard to tell precisely what the drop is from the figure) at the same compression rate.  The different model sizes make it hard to compare, but it's at least been shown that SparseGPT can maintain compelling accuracy even with significant compression that is easy for hardware to exploit.

Thank you for the reminder that there is a formulation for data-free BLAST compression, which is indeed a difference.  I'd only mention here that SparseGPT uses only 128 calibration samples, and roughly the same amount of computation as BLAST (O(hours) for a single model on a single GPU), so it's otherwise a pretty fair comparison.

Further, I'm hesitant to agree that the speedup from BLAST is from reduced FLOPs and not from memory savings, which is also available via unstructured sparsity.  Please read on for more discussion of this point.

[1] https://pytorch.org/tutorials/advanced/semi_structured_sparse.html

Q2. The efficiency of a batched-GEMM compared to a standard GEMM with the same FLOPs is lower, simply because the per-GEMM arithmetic intensity is lower.  I put together a quick PyTorch script to benchmark a moderately-sized (4Kx4K * 4Kx4K) GEMM, as well as a series of batched GEMMs with increasing batch dimensions, and correspondingly decreasing GEMM-K dimensions, resulting in an overall consistent total FLOP count.  I also added a 2:4 semi-structured sparse GEMM to the mix, since I pointed them out in Q1, above.  Results were measured on a single NVIDIA A100.

The latency, in seconds, of 1000 iterations of these workloads are:
GEMM: 0.5415s
2:4 GEMM: 0.3463s
Batched GEMMs:
- b=1: 0.5415s
- b=2: 0.5506s
- b=4: 0.5723s
- b=8: 0.6218s
- b=16: 0.9853s
- b=32: 1.3376s

As you can see, while a small batch dimension is no big difference, as soon as the batch dimension is large enough to start pushing the individual GEMMs into being bandwidth-limited, the overall workload becomes slower.  I also ran a batched gemm with a batch dimension of 16, but with GEMM-K reduced by 32: in this experiment, the total number of FLOPs is halved, but the latency compared to the dense GEMM still increased to 0.6911s.

So, while you observe speedups for the workloads you tested, I believe this is due to only the memory reduction, not the FLOPs savings, and custom kernels may not be able to change this behavior.  Unless you can conclusively show that the FLOPs savings are the reason for a speedup in practice, I might consider focusing more on the compression benefits than the reduced FLOPs.  I fully agree that memory and bandwidth savings are critical components of hardware efficiency, but I still disagree that the results presented support the argument that the reduction in FLOPs (as shown in Table 1 and Figure 6) are of major significance.

Q3, Q5. Thank you for your responses.

Q6. See my point above about an alternate baseline's quality at 50% compression for LLM tasks.  I understand that the goal of BLAST in image diffusion settings is not to match the baseline. Let me clarify my concern: the images generated by BLAST are qualitatively worse than the baseline in a consistent manner.  There is a lack of fine detail (flower, water surface, feathers, jersey), and broad shapes are distorted (hippo, balloon).  If I'm using a model to generate pictures of hot air balloons, I would be unsatisfied with either BLAST or Low-Rank's results in Figure 1, but the baseline model is perfectly acceptable.  I would not judge the compressed model to be useful for the task, but I acknowledge that its results are superior to Low-Rank's.

Q7, Q8, Q9. Thank you for your responses.

Q10. Thank you for these details, they are helpful to understand the workload under test.  In particular, this is a very short input sequence length, which will reduce the time spent in processing the attention in the model, and will focus the time in the linear layers (which are those that BLAST accelerates).  Further, the relatively larger output sequence lengths will require more processing time in the generative phase, which is extremely bandwidth-limited, as opposed to the context phase, which can be math-limited.  This furthers my suspicion that the speedup you see in practice is entirely unrelated to any reduction of FLOPs.


Q11, Q12. Thank you for your responses.