DeltaProduct: Improving State-Tracking in Linear RNNs via Householder Products

We greatly thank the reviewer for the thoughtful and thorough response and apologize for the lack of clarity, typos and small mistakes in the theoretical sections. We modified our draft incorporating the reviewer feedback, which greatly improved the clarity of the paper. After careful checking and as the reviewer points out, none of our theoretical results will change. Due to space constraints, we have to be concise in our response, but will be happy to provide any clarification during the discussion period.

Direct Parametrization Questions

Lines 184-188: The issue discussed in that paragraph concerns bounding the norm of the recurrence. If we choose the second key to be different from the first we cannot in general guarantee the recurrence stability even when $n_h = 1$, as explained in Appendix B.5 for RWKV-7 which also has two different keys. Please also see our first response to Reviewer HNXg, where we expand on the expressivity of RWKV-7 products.

Appendix Issues

Proposition 1

Thank you for identifying this error in Proposition 1. Our corrected proof shows that real eigenvalues are guaranteed if at least one $\beta_i$ is in $[0,1]$ by proving the discriminant $D = [\text{tr}(\mathbf{A}_S)]^2 - 4\det(\mathbf{A}_S) \ge 0$ where $\mathbf{A} = \mathbf{H}_2\mathbf{H}_1$ and $\mathbf{A}_S$ is the subspace spanned by the vectors $\mathbf{k}_1$, $\mathbf{k_2}$.

• If one $\beta_i \in [0,1]$ and the other is in $(1,2]$, then $\det(\mathbf{A}_S) \le 0$, which ensures $D \ge 0$.
• If both $\beta_1, \beta_2 \in [0,1]$, the AM-GM inequality shows $\text{tr}(\mathbf{A}_S) \ge 2\sqrt{\det(\mathbf{A}_S)}$, which implies $[\text{tr}(\mathbf{A}_S)]^2 \ge 4\det(\mathbf{A}_S)$ and thus $D \ge 0$. Hence complex eigenvalues only arise when both $\beta_1, \beta_2 > 1$.

Following your suggestion, we computed explicit bounds on $\cos(\theta)$ for which complex eigenvalues occur:

For the special case of two reflections ($\beta_1 =\beta_2 =2$), this simplifies to $0 < \cos^2\theta < 1$, which holds whenever the keys are linearly independent but not orthogonal.

Unclear Statements and Missing Information

Line 1078-1082: We agree that that section was difficult to follow. The general idea is that for the decoder to support our theoretical constructions, we need the map from all possible states of the RNN and the inputs to the decoder to be injective so that these can be distinguished by the MLP. There are two challenges to make this happen. The first is the dimensionality reduction via the scalar product with the query vector: the vector-valued state that we generally use would be compressed to a scalar after this operation. The second is the RMSNorm operation, which, in the case of a scalar input and $\epsilon = 0$ (in practice also if too small) gives a binary output, so no matter how many states our RNN has, we only have two different inputs to the MLP input. Since for our construction the number of states that we use is always finite, by setting the query as $\**b**/||\**b**||$ and $\epsilon > 0$ (In practice not too small) we can still guarantee injectivity despite these drawbacks. We’ll include a proof of this fact in the revised paper and we’ll be happy to clarify this further during the discussion period if needed. We also noticed a small mistake: $\**b**$ should be defined as \**b** = (b, b^2, \\dots, b^d)^\\top with $b \geq \\delta_{\\max}/\\delta_{\\max} + 1$, $\\delta_{\\max} = \\max_{x,y \\in \\mathcal{S}} |x - y|$, $\\delta_{\\min} = \\min_{x,y \\in \\mathcal{S}, x \\neq y} |x - y|$ and $\\mathcal{S} \\subset \\mathbb{R}$ is the finite set of values of each element of the RNN state matrix.

Line 1064: You are correct. The restriction on operator norm applies here. We’ll rephrase this in the updated version.

Line 1069: We will explicitly state which assumption: “Alternatively, the assumption on $\**H**\_0$ (task-dependent with rank at most $n\_h$)”.

Line 1084: j is the index of the head, which is not related to the content of the hidden state which is $e_i$. We will write that $1 \\leq i \\leq d$ to make this clearer.

Line 1073, Line 1128: Thank you, we corrected it.

Theorem 3

We agree with the reviewer that this section was not clear. We’ll expand and clarify it in the revised version. The main idea of the proof is to divide the tokens into blocks of $m$ matrices. The recurrence of the last layer computes the decomposition for the tokens in the previous blocks, while the current block is computed via the last layer decoder. We’ll be happy to provide further clarification and possibly the complete proof during the discussion period if needed.

Line 1117-1119: We’ll rewrite as “We divide the input sequence into blocks of $m$ elements: we factorize the position $i \in \\{1,\dots t\\}$ into $i = lm + \\tilde i$ where $l \\geq 0$ is the index of the previous block and $\\tilde{i} \in \\{1,\\dots, m\\}$ is the position within the current block (index $l+1$).”

Line 1120: We explain what we mean by “with enough information in xi about previous tokens” later in the text (we’ll clarify this): each token of the last layer should contain the position modulo 2m and the last 2$m$ elements of the input sequence: including both the (partial) current and previous blocks. $x$ is now bold because it is the input of the last layer, which can be a vector, as opposed to the input to the model. It should be $P_{x_{lm+i}}$, not $P_{lm+i}$, we’ll fix it. The matrices $G_{l,m}...G_{l,i+1}$ are not functions of the tokens coming after the current ones because they are part of the “previous block” of tokens. The current block of permutations is computed entirely in the decoder via lookup-table.

Notation and Consistency Issues

Line 97: Thanks for pointing this out, we will make the main paper consistent with the appendix by reverting the order, i.e. using $y_i = x_i \\cdot x_{i-1} \cdots x_1$.

Lemma 2

Line 1146: You are right to highlight the role of the decoder. Our theoretical constructions (as RWKV-7’s) rely on the universality of the MLP, whose required size grows with the complexity of the task (e.g., the group size). We’ll clarify this. The MLP complexity for complex groups with small values of $n_h$ is quite high since to learn $S_n$, the lookup table that it has to learn grows exponentially with $m=(n-1)/n_h$. Indeed, in our experiments (Figure 5 bottom row) the $S_5$ world problem cannot be learned when $n_h=1$ even with 10 layers. However it might be possible to learn the construction if $n_h$ is closer in value to $n$.

Line 1141, 1142, 1145, 1153, 1198, 1247-1248, 1252, 1268-1269: Thank you, we agree and we fixed the typos in our manuscript.

Line 1265: The direction of the angle does not matter for the proof. We will change “vector” for “point” as suggested.

Theorem 4

Line 1152: We will change the word isomorphism to the word assumption. From the context, it will be clear that we are talking about the assumption of the group being isomorphic to the orthogonal or rotation groups.

Line 1159-1161: The fact that we only need products of $n_h$ Householder matrices to represent all elements in $SO(n_h + 1)$ when $n_h$ is even, is a consequence of the Cartan-Dieudonné theorem, which we already mention since it explains DeltaProduct’s expressivity, and the fact that now we are focusing on orthogonal transformations with determinant +1. The determinant of the product of matrices equals the product of the determinants and each Householder has determinant -1. Hence, an even number of Householders cannot be in $SO(n_h + 1)$. We’ll add this clarification.

Lines 1162-1164: The same reasoning can be applied at any step (not only the last one), we will change the index from t to i to highlight this.

Lemma 3

Line 1215-1229: The proof for this will follow the same structure as our clarified proof for Theorem 3, and we will refer to it to maintain conciseness.

Theorem 6

Line 1270: We have added an introductory part to the proof explaining the geometric intuition and defining the matrix representations of reflections and rotations upfront.

Line 1271-1272: We’ll ensure that $d_0$ is always bold. The unproven equation can be directly derived from the definition $d_0 = (1,0)^T$ and the definitions of rotation and reflection matrices, given that the first column of a rotation matrix is identical to that of a reflection matrix (based on our notation).

Lines 1273-1274: We do not use that equality but just use the definition $c_0 =**H**(\\pi/m)d_0$ here. See the above point for the proof of the second equation.

Lines 1274-1275 (Critical Error): You’re right, that was a typo, your expression is correct. The $-2m$ term can be removed due to the periodicity of the rotation matrix, yielding the desired result $**R**(2(i+j))\pi/m)$. We’ll fix this in the revised version.

Additional Issues

Line 1292: Correct. We will revise to remove "exponentially fast".

Lines 1375-1377 (Chomsky Hierarchy): In our experiments, we compared against the results by Grazzi et al. 2025 for DeltaNet, Mamba and mLSTM / sLSTM. Our hyperparameters were chosen based on the ones found to produce the best results on the state-tracking experiments.

Lines 1407-1411 (Language Modeling): Due to compute budget constraints, we decided to use the hyperparameters used by Yang et al. 2024 for similar experiments.

References:

Grazzi, R., Siems, J., Zela, A., Franke, J. K., Hutter, F., & Pontil, M. Unlocking State-Tracking in Linear RNNs Through Negative Eigenvalues. ICLR. 2025

Yang, S., Wang, B., Zhang, Y., Shen, Y., & Kim, Y. (2024). Parallelizing linear transformers with the delta rule over sequence length. Advances in neural information processing systems, 37, 115491-115522.

Given the constraints of the rebuttal, we think the only way to provide the proofs is in the comments. We would like to note that this requires re-formatting due to the absence of macros and the different syntax, which might produce typos. In the following comments we provide the proof and statement of Theorem 3, the proof of Proposition 1 case 3 and the first part of the Decoder implementation section. Please let us know if other parts would be helpful.

Theorem 3.

For any $n \in \mathbb{N}$ there exists a DeltaProduct model with one of the following configurations that can solve the word problem of the symmetric group $S_n$: (i) one layer with $n_h =n-1$ (ii) 3 layers with $n_h>1$ (iii) 4 layers with $n_h=1$.

The construction for (ii) and (iii) requires that the MLP at the second last layer computes a lookup-table of size $2m \times (n!)^{2m}$, function of the last $2m$ input tokens and the position modulo $2m$ with $m = \lceil(n-1)/n_h\rceil$.

Proof. One way to solve the group word problem for the symmetric group $S_n$ is to map each element of the group $g \in S_n$ to the corresponding permutation matrix $\mathbf{P}\_g \in \\{0,1\\}^n$ and then for each input sequence $x_1,\dots, x_t$ with $x_i \in S_n$ compute each element of the output sequence $y_1,\dots, y_t$ as

where $\phi$ is a surjective map from vectors in $\\{1,\dots, n\\}^n$ to the $n!$ elements of $S_n$, which we consider integers for simplicity, i.e., $x_i, y_i \in \\{1,\dots, n!\\}$.

(i). Since a permutation of $n$ elements is a series of at most $n-1$ swaps of $2$ elements, if $n_h = n-1$, then we can solve the problem with a 1-layer DeltaProduct by setting $\mathbf{H}\_0 = \mathbf{u}\_0$, $\mathrm{dec}(\mathbf{H}\_i, x_i) =\phi(\mathbf{H}\_i)$, $\mathbf{B}(x_i) = 0$ ($\mathbf{v}\_{i,j} = 0$), $\mathbf{A}(x_i) = \mathbf{P}\_{x_i}$. The latter is achieved by setting for the $j$-th element in the product $\prod_{j=1}^{n_h}(I - \beta_{i,j} \mathbf{k}\_{i,j}\mathbf{k}\_{i,j}^\top$), either $\beta_{i,j} = 2$ and $\mathbf{k}\_{i,j} = (\mathbf{e}\_k -\mathbf{e}\_p)/\sqrt{2}$ with $\mathbf{e}\_i$ being the $i$-th element of the canonical basis of $\mathbb{R}^n$ (swap element at index $k$ with the one at index $p$), or $\beta_{i,j} = 0$ (identity).

(ii) and (iii). If $n_h < n-1$, then the state transition matrix is not sufficiently expressive to represent all permutations of $n$ elements. However, we can use additional layers to overcome this issue as follows. We divide the input sequence into blocks of $m$ elements: we factorize the position $i \in \\{1,\dots t\\}$ into $i = lm + \tilde i$ where $l \geq 0$ is the index of the previous block and $\tilde{i} \in \\{1,\dots, m\\}$ is the position within the current block (index $l+1$). First, consider the case when $l \geq 1$. Let $\tilde{\mathbf{P}}\_l = \mathbf{P}\_{x_{(l-1)m+m}} \dots \mathbf{P}\_{x_{(l-1)m + 1}}$ be the product of the permutations of the previous block. Since $\tilde{\mathbf{P}}\_l$ is a permutation matrix of $n$ elements, we can factor it into $\tilde{\mathbf{P}}\_l = \mathbf{G}\_{l,m} \cdots \mathbf{G}\_{l,1}$ where we recall that $m =\lceil(n-1)/n_h\rceil$ and each of $\mathbf{G}\_{l,1}, \dots, \mathbf{G}\_{l,m}$ is a product of $n_h$ generalized householder matrices and thus can be a state-transition matrix of our model. We fix one factorization for each possible permutation matrix and we set $\tilde{\mathbf{P}}\_0 = \mathbf{G}\_{0,m}\cdots \mathbf{G}\_{0,1}$, with $\mathbf{G}\_{0,i} = I$ to handle the case when $l=0$.

Now let $\mathbf{x}\_i$ be the input of the last layer. if $\mathbf{x}\_i$ contains enough information about previous tokens (as we specify later), we can set the recurrence and decoder of the last layer as

where $\mathbf{H}\_0 = \mathbf{u}\_0$, $\mathbf{B}(\mathbf{x}\_i) = 0$, $\mathbf{A}(\mathbf{x}\_i) = \mathbf{G}\_{l,\tilde i}$, using the construction at point (i) since $\mathbf{G}\_{l,\tilde i}$ is a product of at most $n_h$ Householers. Note that $\mathbf{H}\_i$ contains the product of the input permutations only up to token $x_{(l-1)m}$ and a partial computation of previous block of permutations $\tilde{\mathbf{P}}\_l$. Hence, the decoder completes the computation by applying two additional components: (1) the remaining transformations $\mathbf{G}\_{l,\tilde{i}+1}$ through $\mathbf{G}\_{l,m}$ needed to complete $\tilde{\mathbf{P}}\_l$, and (2) the actual permutations from the current partial block $\mathbf{P}\_{x_{lm + 1}}$ through $\mathbf{P}\_{x_i}$.

We thank the reviewer for the positive feedback.

I'm curious about a FLOPs comparison in the scaling analysis in Figure 10. To match parameters, the number of heads was reduced, but how do these variants compare w.r.t. FLOPs?

This is a great question. The FLOPs of the recurrence scale linearly in $n_h$ such that the complexity is $O(n_h \times (LCd + Ld^2))$ (Yang et al., 2024), where $L$ is the sequence length, $C$ is the chunk size and $d$ is the dimension of keys and values. The recurrence computation also requires $n_h \times L/C$ sequential steps. However, in our setup the total FLOPs of the model are dominated by the large matrix multiplications in the input/output projections and the MLP, rather than the recurrence itself, and so, the contribution of the recurrence to total model FLOPs is minor.

In our parameter-matched scaling experiments (Figure 10 and Table 4), we adjusted the number of heads to equalize parameters. For instance, at the ~800M scale, we compare DeltaNet (16 heads) with DeltaProduct$_2$ (12 heads). The increase in FLOPs from the additional key/value projections and the doubled recurrence in DeltaProduct$_2$ is largely offset by the reduction in FLOPs from having smaller query, key, and value projections due to the fewer number of heads. Our results in Figure 13 in Appendix C.4 also support this. There we can see that for non-parameter matched models, the training throughput scales sub-linearly in $n_h$. We conducted additional throughput experiments with an improved model implementation that fuses the projection layers and causal 1d-convolutions. Here we find that DeltaProduct$_3$ with 805M parameters achieves 86% (compared to 75% without fusing projections and causal 1d-convolutions) of the training throughput of DeltaNet with 805M parameters at 4096 sequence length and a per-device batch size of 12 (on a NVIDIA H100-94GB). Moreover, we expect this to further improve with an optimized custom triton kernel (we are using DeltaNet’s).

Did DeltaProduct display any training instabilities?

No, we did not encounter any training instabilities with (Gated) DeltaProduct. We attribute this to the inherent stability of the recurrence, as the operator norm of the transition matrix is bounded by 1, which keeps the gradients well-behaved. The training loss curves for $n_h=1,2,3$ in Figure 14 of the supplementary material demonstrate this stability empirically.

In Figure 5 bottom 2nd and 4th plots, an increased number of layers seems to hurt performance on the state tracking tasks. You addressed it, but isn't this unexpected? Is there an explanation? Also, the number of layers is not consistent in that row, which reduces the readability of the figure as a whole.

Indeed, increasing layers from 8 to 10 reduces performance for the $S_4$ and $S_5$ groups. However, we see that these experiments are quite challenging and require multiple seeds to find a run that could fit the training context length of 128 permutations and show extrapolation to 512 permutations: we use 5 seeds per configuration yet sometimes still couldn't fit the context length. This observation reinforces our claim: increasing expressivity horizontally via $n_h$ offers a more parameter efficient and stable optimization path for these state-tracking tasks than increasing it vertically via depth. We thank you for pointing out the inconsistent legend labels in the colormap; we will fix the colormap in the revised figure to improve readability.

In Figure 9, the Gated DeltaNet effective ranks are extremely jagged. Why is this?

Both Gated DeltaNet and Gated DeltaProduct produce very jagged effective ranks. We found that both models learned to set the forget gate close to zero at the beginning of sequence tokens, effectively erasing the memory for the new sequence. Since the forget gate is a scalar for each hidden state per head the reset is very drastic.

line 42 should be an em-dash (https://en.wikipedia.org/wiki/Dash)

Thank you for pointing this out, we have fixed it in our manuscript.

Thank you again for helping us improve the paper.

References

S. Yang, B. Wang, Y. Zhang, Y. Shen, and Y. Kim. Parallelizing Linear Transformers with the Delta Rule over Sequence Length. In A. Globerson, L. Mackey, D. Belgrave, A. Fan, U. Pa- quet, J. Tomczak, and C. Zhang, editors, Proceedings of the 37th International Conference on Advances in Neural Information Processing Systems (NeurIPS’24), 2024

We thank the reviewer for the feedback and address each point below.

The model's core design using Householder products enforces a strict norm ||A(xi)|| <= 1 on its state-transition matrix. While this guarantees stability, it is also a fundamental expressive bottleneck that prevents the model from performing norm-increasing operations, which are crucial for amplifying signals and are a key reason why competing models like RWKV-7 achieve stronger theoretical guarantees.

We agree this is a crucial trade-off. We intentionally designed DeltaProduct with Householder products to prioritize guaranteed stability ($\Vert A(x_i)\Vert \leq 1$), a criterion for stable training on long sequences. This design choice directly contrasts with the potential instability of models like RWKV-7. As we demonstrate in Appendix B.6, the RWKV-7 recurrence becomes unstable in the very cases that grant it higher expressivity (i.e., when vectors $a_i$ vary per token); we show that the product of just two such matrices can have a spectral radius > 1.2.

However, the core product-of-matrices mechanism we propose is flexible. It can incorporate other matrix families, including the norm-increasing ones from RWKV-7 Indeed, a minor extension of our theory, building on Merrill et al. (2024) (Theorem 5.2) and Peng et al. (2025) (Lemma 3), shows that this generalized framework has stronger theoretical guarantees. For instance, a single-layer model using $n_h = n$ RWKV-7 matrix products can recognize any regular language accepted by a finite state automaton with $n$ states, while if $n_h \geq 2$, a three-layer model can recognize any regular language. We will add these results to the revised paper.

To validate this, we ran new experiments with products of RWKV-7 type state-transition matrices. We used $(\mathbf{I} - 2 * \mathbf{k} (\mathbf{k} \odot \mathbf{a})^T)$ as state-transition matrix where $\mathbf{k}$ is unit norm and $\mathbf{a}\in [0, 1]$ by transforming through a sigmoid. We find that a single layer using products of 2 RKWV-7 matrices obtains comparable results to DeltaProduct$_2$ on $S_3$. We will add new experiments for the other permutation groups using products of RWKV-7 state-transition matrices to a revised version of our manuscript.

Increasing the model's expressivity via the hyperparameter $n_h$ (the number of Householder products) is inefficient. It leads to a linear growth in both model parameters and, more critically, the sequential computation time of the recurrence, making it nh times slower per token. This architectural choice limits the practical scalability of using a large nh to achieve higher performance.

We see the linear increase in sequential computation with $n_h​$ as a tunable trade-off rather than a fixed limitation. Our experiments in Figure 8 show that a small increase in $n_h$​ (e.g., to 2 or 3) yields substantial gains in length extrapolation, making it a highly effective trade-off for many applications. Despite the theoretical linear increase in computational complexity, we find that in practice the training time increases sub-linearly with $n_h$ as shown in Figure 13 in Appendix C.4. The scaling-rate is also dependent on training sequence length and batch size and could be further improved by a custom triton kernel. After submission, we benchmarked the training throughput of the parameter-matched 805M parameter DeltaNet and DeltaProduct$_3$ on an improved implementation of our model with fused projections and causal 1d-convolutions (see also the response to reviewer RGKp). We find that DeltaProduct$_3$ achieves 86% (compared to our previous 75% without fusing projections and causal 1d-convolutions) of the training throughput of DeltaNet at 4096 token sequence length and a per-device batch size of 12.

The number of operations, nh, is a fixed, global hyperparameter. The model applies the same amount of costly computation to every token, regardless of its importance. This is wasteful for simple tokens (e.g., punctuation) and may be insufficient for complex ones (e.g., resolving a long-range dependency), representing a missed opportunity for adaptive computation that could allocate resources more efficiently.

Thank you for the suggestion. Adaptive computation is a promising direction that could greatly enhance the model's efficiency by allocating computation dynamically. While this falls outside the scope of our current work, which focuses on establishing the fundamental properties of the Householder product structure, we agree it is a valuable avenue for future research. We will add a discussion of adaptive computation, such as using a Mixture-of-Experts (MoE) approach to select $n_h$ per token, to our future work section.

References:

Merrill, W., Petty, J., & Sabharwal, A. (2024, July). The Illusion of State in State-Space Models. In International Conference on Machine Learning (pp. 35492-35506). PMLR.

Peng, B., Zhang, R., Goldstein, D., Alcaide, E., Du, X., Hou, H., ... & Zhou-Zheng, C. (2025). Rwkv-7" goose" with expressive dynamic state evolution. arXiv preprint arXiv:2503.14456.