Qwen Team Releases FlashQLA for High-Performance Linear Attention

Once large model sequence lengths exceed one million, the O(N²) overhead of traditional softmax attention becomes a bottleneck. Models like Qwen are shifting toward hybrid architectures that combine linear and softmax attention. FlashQLA is a high-performance GPU operator library specifically designed for GDN (Gated Delta Network) linear attention, addressing challenges such as low parallelism, redundant memory access, and SFU overload.

The implementation focuses on hardware-friendly algebraic rewrites tailored to GDN structural characteristics. It proposes co-optimization through sequence-dimension chunking and warp specialization. By leveraging sliding window properties for reuse, it significantly reduces memory bandwidth pressure.

This content is suitable for large model system engineers, GPU operator developers, and inference engine architects.

As large model sequence lengths advance from tens of thousands to millions, the memory and computational overhead of O(N²) softmax attention has become a hard constraint. In recent years, mainstream models represented by Qwen and Kimi have gradually shifted toward hybrid architectures combining "linear attention + softmax attention." However, efficient GPU implementations of linear attention still face systemic challenges, including insufficient parallelism, redundant memory access, and excessive load on Special Function Units (SFUs).

This article summarizes the talk "High-Performance Optimization Practices for Linear Attention in Qwen Series Models" delivered by Cheng Zhangruidong, an algorithm engineer at Alibaba's Qwen team, at AICon 2026 Shenzhen. He detailed FlashQLA—a high-performance GPU operator library for GDN (Gated Delta Network) linear attention—covering its positioning and background, sequence parallelism strategies, utilization of sliding window properties, hardware-friendly algebraic rewrites, and warp specialization implementation. The talk also provided performance benchmarks and updates on ecosystem integration.

The following is a transcript of the speech (edited by InfoQ without altering the original meaning).

Positioning and Background: Why FlashQLA Was Created

The core motivation for linear attention stems from a simple observation regarding matrix multiplication chains. The computation flow for softmax attention involves first performing matrix multiplication between Q and K, then applying row-wise softmax normalization, and finally multiplying with V. Even though FlashAttention avoids explicitly instantiating the N×N P matrix, it still requires completing row-wise, data-dependent softmax operations on P. Consequently, the computational complexity remains O(N²d).

However, if we temporarily ignore the non-linear softmax operation and reconsider the process purely from the perspective of matrix multiplication chains: Both Q and K have shapes defined by sequence length multiplied by head dimension. Computing QKᵀ first yields an N×N matrix, which is then multiplied by V, resulting in a total complexity of 2N²d. Conversely, if we compute KᵀV first, we obtain a d×d matrix, which is then multiplied by Q. In this case, the complexity of both matrix multiplications drops to Nd². This means that as long as we can bypass the intermediate row-wise softmax, there is potential to reduce the time complexity of attention calculation from O(N²d) to O(Nd²).

Following this line of thought, the most direct approach is to retain the causal lower-triangular matrix while removing softmax. Causality arises from temporal dependencies and cannot be discarded; however, softmax can be removed. This results in a recursive form: instead of computing KᵀV for the entire sequence at once, we accumulate token by token. The outer product of K and V for each token is written into a d×d state matrix, which then propagates to the next step to ultimately compute the output. This constitutes the fundamental O(Nd²) logic of linear attention.

Structurally, this recursive formula resembles an RNN. Its advantage is that time complexity grows linearly with sequence length, but the cost is equally apparent: obtaining the state for the next token strictly depends on the state of the previous token, creating strong data dependencies. From 2022 to 2025, academia derived numerous linear attention variants based on this framework. The first variant truly capable of scaling up was the Gated Delta Network (GDN), proposed by Yang Songlin in early 2025. The state update formula for GDN evolved to be quite complex, introducing two key parameters: the update gate beta and the forget gate alpha. The forget gate alpha plays a critical role in subsequent optimizations. Currently, GDN represents a well-balanced linear attention solution in terms of expressiveness and efficiency, offering significant improvements over original linear attention. Since then, several DeltaNet variants have emerged, such as KDA in Kimi K3.

The hybrid architecture currently dominated by Qwen and Kimi shares a common design philosophy: replacing most attention layers with linear attention. A notable algorithmic flaw in linear attention is that regardless of sequence length, all information must be compressed into a d×d state space. This inherently limits information capacity, making it difficult to achieve effective long-range recall. To compensate for this, models insert original softmax attention every few layers. The Qwen series adopts a 3 structure consisting of three GDN layers followed by one softmax attention layer, while Kimi K3 uses three KDA layers followed by one softmax attention layer.

In standard prefill operations for GDN modules, calculating recursively token-by-token according to the formula leads to a specific issue: the k and v vectors in the formula rely on vector dot products and outer products, which cannot efficiently utilize GPU tensor cores. The solution is to extend the recursion step size from single tokens to the input dimensions of matrix multiplication units. This requires introducing a mathematical transformation known as the WY representation. After this transformation, the recursion shifts from token-by-token to chunk-by-chunk, with a typical chunk size of 64.

The process consists of four steps: First, compute $K K^T$ for the $K$ matrix within each chunk to obtain a chunk_size × chunk_size matrix and invert it. Second, use the inverse matrix to perform intermediate projections on $K$ and $V$, yielding $W$ and $U$. Third, execute the state recurrence step-by-step across chunks. Fourth, calculate the output based on the state $S$ derived from the recurrence. Notably, except for the third step, the first, second, and fourth steps have no data dependencies, allowing work on different chunks to be executed independently in parallel. However, the third step requires the state from the previous chunk to compute the state for the current chunk, which is the critical bottleneck for parallelism.

In the existing implementation of the FLA repository, these four steps correspond to four independent operators, with three variables—$W$, $U$, and $S$—instantiated as intermediates. These variables exist only as temporary values during the recurrence process; inference frameworks do not need to retain them. Nevertheless, they are written back to HBM and then read by the subsequent operator, resulting in pure HBM bandwidth overhead. Additionally, there is a parameter gamma, representing the cumulative product of gates within each chunk, where "gate" refers to the previously mentioned alpha.

This four-step forward pass has three main pain points. The first is insufficient parallelism. The third step requires that the state of each chunk depends on the previous chunk, meaning each attention head can only launch one thread block (CTA), corresponding to a single SM. Taking smaller Qwen models as an example, the number of linear attention heads might be as low as 32, whereas H100 or GB series cards often have over a hundred SMs, leading to severely underutilized GPU resources. This limitation on parallelism also poses challenges for Tensor Parallelism (TP) optimizations. The second issue is state transfer overhead. The three intermediate variables $W$, $U$, and $S$ must be written back to HBM and read by subsequent operators, yet they do not enter the persistent state of the inference framework, effectively wasting bandwidth. The third problem is the capability ceiling of the Triton operators used in FLA. Triton cannot finely control the internals of operators and does not provide interfaces for hand-written warp specialization. When multiple matrix multiplications are merged into a single operator, along with chunk-by-chunk recurrence, the structure resembles a traditional matrix multiplication pipeline, making it difficult for Triton to achieve ideal pipeline orchestration. Therefore, we need lower-level tools to manually orchestrate reasonable pipelines to fully utilize various hardware units on the GPU.

Core Technical Principles

To address insufficient parallelism, a natural approach is to parallelize along the sequence dimension. Prior work in the FLA repository already includes GDN sequence parallelism. The core problem can be formulated as follows: If we split the sequence in the middle and start the recurrence for the latter half from a zero-initial state, the resulting state will inevitably be incorrect. The question is: how do we correct this error?

Let $S$ be the correct state sequence, and $S^*$ be the state sequence obtained by starting the recurrence from an incorrect initial state. By subtracting the recurrence formulas for both, we derive an important property: the difference between states of adjacent chunks follows a matrix multiplication relationship. Specifically, the difference at the next step equals the difference at the previous step multiplied by a $d_K \times d_K$ matrix. Further leveraging the associativity of matrix multiplication, we can independently compute correction matrices for each chunk derived from the incorrect state. In other words, instead of multiplying the correction matrix onto $S$ immediately, we can first compute the correction matrix corresponding to each chunk and then apply a single matrix multiplication at the end to complete the correction. This observation—that the recurrence of the $d_K \times d_K$ matrix can be "associated"—is the theoretical prerequisite for sequence parallelism.

However, directly storing a $d_K \times d_K$ correction matrix for every chunk is too costly. For example, with $d_K=128$ and chunk_size=64, storing a $128 \times 128$ matrix per chunk essentially reintroduces the state transfer overhead that was previously saved through operator fusion. In the actual implementation, we adopt a "compute-for-memory" strategy: each CP rank needs to save only one correction matrix. Specifically, within each CP rank, we independently compute the incorrect states $S^*$ and correction matrices $M$ for all chunks but retain only the result for the last chunk of that rank. Then, we perform a state-dependent recurrence across the length of the CP size to sequentially obtain the precise corrections for each rank. Finally, we recalculate the true $S$ and $O$. This scheme implies that $S$ is computed twice—the first pass yields $S^*$, and the second yields the true $S$—resulting in some computational waste, but it is highly efficient in terms of memory usage. Logically, this solution does not distinguish between intra-card sequence parallelism and inter-card sequence parallelism. Currently, the FLA repository only supports inter-card sequence parallelism, but the underlying principles are identical, and there is no algorithmic necessity to explicitly differentiate between the two.

After solving the parallelism issue, can we directly implement a fully fused operator? The answer is no. This is because certain computations within a fully fused operator depend on prior data required for sequence parallelism, while others do not. The most critical dependent quantity is matrix inversion. Matrix inversion is not cheap on GPUs, and computing the sequence parallelism correction matrices requires the result $W$ after inversion. If inversion were moved to a subsequent operator, calculating the prior data dependencies $S^*$ and $M$ for sequence parallelism would require recomputing the inversion, which is highly inefficient. Therefore, we choose to split the workflow after the first step of the four-step plan: matrix inversion is isolated, the remaining three steps can be fused, and operations related to sequence parallelism are inserted between these two parts.

The key difference between intra-GPU sequence parallelism and inter-GPU sequence parallelism is that the degree of parallelism across GPUs is manually fixed, whereas for intra-GPU parallelism, we do not know how many chunks a given sequence length should be divided into. This number can theoretically be estimated. Returning to the previously mentioned three-step parallel scheme: Step 1 and Step 3 can run independently across each CP rank, with a recurrence length equal to the number of chunks $L$ within each rank; Step 2 has a recurrence length of $N/L$, which corresponds to the CP size. In the simplified "one wave" model, the time for Step 1 and Step 3 is proportional to $L$, while the time for Step 2 is proportional to $N/L$. Thus, the total time $T = \lambda_1 L + \lambda_2 N/L$. Minimizing this expression yields the optimal $L$, which follows a linear relationship with $\sqrt{N}$. The specific values of $\lambda_1$ and $\lambda_2$ depend on the actual performance of operators on specific hardware and can be obtained through tuning. We are currently optimizing this model, as the "one wave" simplification is too coarse. Future work will introduce more parameters for each operator to achieve a more precise estimation of the sequence parallelism scheme.

Although parallel schemes can be inserted, the insertion itself introduces overhead, so it is necessary to determine when to enable them. An obvious condition is to trigger when the total workload is small enough, i.e., when the GPU's SMs are clearly underutilized. However, there is another easily overlooked hidden condition: if the sequence itself is not long enough, inserting sequence parallelism offers little benefit in terms of recurrence. This is not a simple AND or OR relationship. If we plot batch size multiplied by head count on the x-axis and sequence length on the y-axis, the reasonable boundary for enabling auto-CP should be a quadratic curve. The current implementation uses a simplified trapezoidal partitioning, which will later be replaced with a more precise scheme.

GDN itself possesses a mathematical property that was previously underutilized: the sliding window property. Observing the original GDN recurrence formula: $\alpha$ is a real number between 0 and 1, $\beta$ is also a real number between 0 and 1, and $K$ is L2-normalized, meaning the magnitude of $K$ for each token is constrained. Therefore, the determinant of the matrix $I - \beta K K^T$ has a well-defined numerical range. Since each step accumulates a multiplication by $\alpha$, the influence of each state on subsequent states decays exponentially with $\alpha$. Specifically, the impact of $S_t$ on $S_{t+1}$ is attenuated once by $\alpha$, and its impact on $S_{t+8}$ is attenuated eight times. When the window is sufficiently long, the influence of early states on the current state can be considered zero. In other words, as long as $\alpha$ for a head is not constantly 1, a sliding window will inevitably exist, with the window size determined jointly by the distribution of $\alpha$ and the gate data distribution.

We performed statistics on trained models and found that approximately 60% to 80% of linear attention heads have an $\alpha$ that is not constantly 1. This proportion is quite reasonable: the gate parameter $\alpha$ added in GDN is precisely intended to enhance expressiveness. If the gate does not function effectively, the model degenerates into DeltaNet, which lacks sufficient expressiveness. In other words, as long as training causes the gate to play a role, the sliding window phenomenon will emerge. For these heads exhibiting the sliding window property, the sequence parallelism scheme can be significantly simplified: instead of starting the recurrence from the initial state at index 0, we only need to backtrack a sufficiently large window size from the start of each chunk to directly obtain a starting state with negligible error. This aligns with the sequence parallelism approach used in sliding window attention.

More advantageously, this sliding window scheme can reuse operators from the aforementioned sequence parallelism scheme. The sequence parallelism scheme necessarily requires computing $S^*$, but for sliding attention heads, $S^*$ is already the exact $S$ and can be used directly, eliminating the need to compute the $M$ matrix and perform the correction in Step 2. The operator for Step 1 can be fully reused. To decide which heads adopt the sliding window optimization, an additional lightweight kernel must be launched to scan the gate data distribution for each head. For GDN, the gate is a single real number per token, so the scanning cost is comparable to gate accumulation, making the overhead negligible.

Regarding the GDN recurrence process, we performed two types of hardware-friendly algebraic rewrites: one targeting the load on SFUs (Special Function Units), and the other targeting the time complexity of chained matrix multiplications.

First, let us examine the SFU load. In the GDN calculation formulas, lowercase $\gamma$ represents the cumulative product of $\alpha$, while uppercase $\Gamma$ is the outer product matrix of lowercase $\gamma$ within a chunk. Formally, the element at row $s$ and column $t$ of uppercase $\Gamma$ equals $\gamma_{i_s} / \gamma_{i_t}$. Directly calculating the reciprocal of lowercase $\gamma$ numerically poses a risk of overflow, especially when the chunk size is large, as $\gamma$ can easily approach zero, potentially exceeding BF16 precision limits. A safer approach is to explicitly expand the calculation, using stored $\log \alpha$ values and replacing multiplication with addition: $\gamma_{i_s} / \gamma_{i_t}$ equals $\exp(g_{i_s} - g_{i_t})$. Consequently, computing each term of this $C \times C$ matrix requires one EXP operation, totaling $C^2$ EXP operations. This is the source of the SFU load.

The key observation is that during matrix inversion, uppercase $\Gamma$ can be factored out. In the original inversion formula, uppercase $\Gamma$ appears after $K K^T$. Since the diagonal elements of uppercase $\Gamma$ are always 1, it can be merged with the identity matrix, allowing uppercase $\Gamma$ to be moved outside.

Formally, uppercase $\Gamma$ can be written as lowercase $\gamma$ multiplied by the identity matrix and then by the inverse of lowercase $\gamma$. During lower triangular matrix inversion, this diagonal scaling can be factored out entirely. Let the resulting inverted matrix be denoted as $A'$. Subsequently, when calculating $W$ and $U$, uppercase $\Gamma$ is multiplied back onto $A'$. Since $W$, $U$, and $S$ inherently require calculating uppercase $\Gamma$ in subsequent operators, this approach avoids the problem of separately computing $O(C^2)$ EXP operations in two different operators, significantly reducing the pressure on SFU units.

The second category of rewrites focuses on chained matrix multiplications. Taking the sequence parallel transition matrix $M$ as an example: originally, $M$ is the cumulative product of several $128 \times 128$ matrices. The time complexity of each matrix multiplication is $128 \times 128 \times 128$, which is computationally expensive. By leveraging the aforementioned large Gamma extraction technique, we can derive an inverse matrix $A'$ that excludes the large Gamma terms. Re-deriving the recurrence formula for $S$ allows us to omit the term in the recurrence where the last gamma of each chunk is divided by all gammas. At this point, the variables required for the recurrence can be organized into $X$, $Y$, and $Z$ forms, which also apply to $M$. Intuitively, since $M$ and $S$ are very similar, there must exist a way to compute $M$ with the same computational complexity as $S$. After restructuring, the gate cumulative product term within each chunk of the $M$ recurrence is factored out of the parentheses. Inside the parentheses, $X$ is the inverse matrix multiplied by $\beta K K^T$, ensuring bounded numerical values; $Z$ is $K$ multiplied by $M$, which is also bounded. Consequently, the numerical range of the entire $M$ recurrence is controllable. When combined with the gamma decay term, the $M$ matrix itself exhibits a decay property proportional to the cumulative gate product. This effectively validates the sliding window nature of GDN from another perspective.

Similar optimizations for chained matrix multiplications apply to both forward and backward recurrences. For the forward pass, since $W$ and $U$ contain similar terms, we can merge the alpha coefficients. By first calculating $V - S$ and then multiplying by gamma, we reduce the number of matrix multiplications. This approach already existed in earlier fully fused GDN implementations; our main contribution was porting and integrating it. Optimizations for the backward recurrence are more fragmented. One representative work involves backpropagating the gradient of $K K^T$ to $K$. This requires computing $A \times K$ plus $DA^T \times K$. Here, we can combine like terms for $DA$ and $DA^T$ before multiplying by $K$, thereby reducing the matrix multiplication load. Although performing a shared memory transpose on 32-bit floating-point matrices places higher demands on the already tight shared memory bandwidth, this merging operation benefits subsequent reduce operations for gate and beta gradients throughout the backward flow. Therefore, it represents a holistic optimization choice.

After algebraic rewriting and fusion, the computational volume of the last three operators remains smaller than the memory access volume. Using our custom units, the memory access volume for a $64 \times 64$ BF16 matrix within one chunk is defined as one unit. The total memory access volume for $Q$, $K$, $V$, $A$, and the output $O$ is approximately 18 units. On the H200, after mapping the tensor core computation bandwidth of WGMMA equivalently, the computational volume for matrix multiplications is approximately 30 units. This means the computational volume is not significantly lower than the HBM memory access volume. However, note that this calculation uses direct HBM bandwidth without considering L2 cache hits. For GDN prefill, L2 hit rates are indeed low because each block of $Q$ and $K$ is not reused across different SMs after being consumed, unlike in standard matrix multiplications. Overall, this remains a memory-bound operator.

Based on this, we designed a ping-pong style two-stage pipeline. Traditional multi-stage pipelines often aim to leverage L2 cache, but since the GDN operator has low L2 hit rates, a two-stage pipeline is sufficient. Regarding registers, 64K 128-bit registers are evenly distributed among four warp groups, leaving approximately 128 registers per thread, which is quite tight. Specifically, the $S$ matrix is $128 \times 128$; when distributed across each warp group, each thread exactly fills its 128 32-bit registers. Shared memory is also constrained. If double buffering is applied to inputs with head dimensionality such as $Q$, $K$, and $V$, the total requirement is about 9 units. Each unit corresponds to 28KB, which just fits the ping-pong structure, but intermediate variables require extreme shared memory reuse.

In terms of pipeline division of labor, we assign responsibilities based on warp groups. On Hopper architecture, WGMMA instructions are issued by the entire warp group, so orchestration must be done at the warp group granularity. The design utilizes four warp groups: one producer warp group is responsible for transferring $Q$, $K$, and $V$ from HBM via TMA, while also acting as a consumer to write the output $O$ back to HBM. The remaining three consumer warp groups handle the computation of $V'$, $O$, and $S$ respectively. Why dedicate a single warp group exclusively to $S$? Because $S$ needs to reside permanently in registers, making exclusive allocation the most logical choice. Since $V$ and $O$ involve many temporary intermediate variables, allocating them to separate warp groups ensures sufficient register space. Through this warp specialization design, we overlap data transfer, tensor core computation, CUDA core computation, and SFU computation entirely. As a memory-bound operator, the computational load is not excessively high, but NCU profiling shows that the memory access load exceeds 70%.

Performance and Benchmarks: 2-3x Forward · 2x Backward

In terms of final performance, we achieved a 2x to 3x speedup for the forward pass and approximately a 2x speedup for the backward pass. The optimization logic for the backward pass is largely consistent with that of the forward pass, including an intra-card sequence parallelism scheme based on sliding windows. Although warp specialization in backward operators is more complex, the underlying framework remains the same. This operator repository is designed for both training and inference; currently, it has been integrated into Qwen’s internal training and inference pipelines. The acceleration effect from auto CP is particularly significant in scenarios with high TP parallelism or smaller models. Compared to FLA’s Triton implementation, our solution delivers a 2x to 3x speedup in common scenarios, approaching nearly 3x in selected favorable cases. When compared to earlier versions of FlashInfer, our approach offers greater compression potential at higher TP levels, primarily due to resolving SM utilization issues. Benchmarks cover common configurations such as the Qwen3.5 / Qwen3.6 series, with GDN head counts ranging from a minimum of 8 to a maximum of 128, demonstrating reproducible speedups across different head configurations. Regarding hardware support, the initial design targeted Hopper, followed by subsequent support for SM100, SM103, SM120, and SM121, covering the entire Blackwell family.

Engineering Progress and Open Source Ecosystem

FlashQLA aligns with FLA’s interface conventions in engineering practice. It provides high-level APIs with signatures identical to those of FLA, allowing code originally using chunk_gated_delta_rule to switch directly via from flash_qla import chunk_gated_delta_rule. An auto CP toggle has been added to the interface. When enabled, we guarantee no performance degradation, but cannot ensure special requirements such as deterministic inference or consistency between training and inference; users with such needs can manually disable it. The project has been published on PyPI and can be installed directly via pip. Since version v0.1.2, FlashQLA has become the official backend for FLA, meaning that after installing FLA, users can directly invoke FlashQLA provided their hardware meets the conditions. The community has also helped us complete partial integration with vLLM.

On the project timeline, the first open-source release occurred on April 24, 2026. In June, we completed intra-card sequence parallelism for the backward pass and added SM100 support. In July, with community assistance, we completed SM120 support and officially became the FLA backend. Currently, we are advancing work in two areas: first, extending inter-card sequence parallelism. Since sequence parallelism operators do not distinguish between intra-card and inter-card scenarios at the algorithmic level, this extension is a natural progression, primarily targeting ultra-long sequence inference in PD-disaggregated setups and ultra-long data pre-training. Second, developing new auto CP algorithms. The current one-wave model is overly simplified; future work will involve more detailed modeling of operator performance to achieve more precise decisions for sequence parallelism splitting.

Application Scenarios and Outlook

Overall, FlashQLA aims to cover scenarios such as large model pre-training, edge-side inference, and long-sequence deployment with TP, thereby completing the open-source model ecosystem. All large models are supported; small models benefit from optimizations in edge-side inference and single-card training, while large models gain optimizations in sequence parallelism when TP is enabled, along with support for variable-length sequences. There are three core contributions: First, auto CP intra-card sequence parallelism, which leverages the sliding window property of GDN that was previously not systematically emphasized. Second, hardware-friendly algebraic rewrites, including extracting large Gamma values to reduce SFU computation and optimizing matrix multiplication order. Third, hand-writing multiple warp-specialized operators using TileLang. Overall performance is reproducible on both Hopper and the Blackwell family, achieving 2x to 3x speedups for both forward and backward passes, with out-of-the-box usability as an official backend. Finally, we thank the teams behind Flash Linear Attention, TileLang, and FlashInfer, as well as prior work on GDN sequence parallelism.

Conference Recommendation

The 2026 AICon Global Artificial Intelligence Development and Application Conference (Beijing) has officially launched! Featuring 12 major tracks covering model capabilities, inference infrastructure, Agent engineering, software development, product design, and commercialization practices. The conference brings together frontline technical experts and industry practitioners to share deep insights on cutting-edge trends and practical experience. Call for papers is now fully open; we sincerely invite you to take the stage and share your hands-on experiences. AICon 2026 Beijing looks forward to joining you on this journey.