Fast convolution algorithm for state space models

Gregory Beylkin

REVIEW 2 major objections 4 minor 12 references

Fast convolution algorithm for state space models

T0 review · 2 major / 4 minor · reviewed 2026-08-12 · deepseek-v4-flash

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗ Review JSON

Pith's one-line read An unconditionally stable cascade algorithm applies state-space transfer functions in time at about twice the cost of the standard recursion.

desk verdict A correct variant of the author's old FIR-cascade trick for stable matrices, with an advertised extension to |λ|>1 that is false and a complexity claim off by a log factor. read the letter →

arxiv 2411.17729 v3 pith:7IAITW6F submitted 2024-11-22 math.NA cs.AIcs.NA

Gregory Beylkin This is my paper · ORCID

classification math.NAcs.AIcs.NA MSC 65F6093C05

keywords statespacemodelslineartime-invariantsystemstransferfunctioncascadealgorithmmatrixpolynomialstructuredmatricesHiPPOIIR/FIRapproximation

verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

This paper is about applying the transfer function of a linear time-invariant state space model as a convolution in time, without running the recursion that becomes unstable when the state matrix has an eigenvalue close to or above 1 in modulus. The proposed algorithm truncates the infinite product (I - $z^{{-1}}$A)^{-1} = ∏_{n≥0}(I + ($z^{{-1}}$A)^{2^n}) to N+1 factors, forming a matrix polynomial of degree $2^{{N+1}}$-1. In the time domain this polynomial acts as a cascade that reuses intermediate matrix-vector products, so L output samples cost at most about 2L matrix-vector multiplications, compared with L for the standard recursion. The paper's advertised payoff is unconditional stability: structured approximations of the state matrix that would break the eigenvalue condition no longer have to be rejected, so cheaper O(m) or O(m log m) matrix-vector multiplication schemes become available. A numerical example with a HiPPO matrix of size 100 shows that N=14 achieves degree 32767 with 15 matrix-vector multiplications.

What carries the argument

The load-bearing object is the infinite product identity (I - $z^{{-1}}$A)^{-1} = ∏_{n=0}^{∞} (I + ($z^{{-1}}$A)^{2^n}), valid on the unit circle when the spectral radius of A is below 1. The algorithm truncates this product after N+1 factors, giving a matrix polynomial of degree $2^{{N+1}}$-1, and applies it via a cascade: at level n, the matrix $A^{{2^{n-1}}$} acts only on entries spaced $2^{{n-1}}$ apart, reusing earlier intermediate vectors. That re-use is what caps the total work at about 2L matrix-vector products for L time steps, independent of the polynomial degree.

What would settle it

For scalar systems, set a = 1.01 (eigenvalue above 1), B = C = 1, D = 0, and compare the cascade output H_N with the true stable target H(z) = 1/(1 - $z^{{-1}}$ a_0) for a_0 = 0.99. The tail (a/z)^{$2^{{N+1}}$} grows with N, so the error ‖H - H_N‖ on |z| = 1 will not decrease as N increases; observing that would refute the unconditional-stability claim.

Watch

Extended reading notes

Core claim

The paper's central claim is that the rational transfer function H(z) = C(I - $z^{{-1}}$A)^{-1}B + D of a discrete-time LTI system can be replaced, to any user-selected accuracy, by the matrix polynomial H_N(z) = C ∏_{n=0}^{N}(I + ($z^{{-1}}$A)^{2^n})B + D, and that this polynomial can be applied as a cascade in the time domain using only about 2L matrix-vector multiplications for L output samples. Because the cascade only ever multiplies by powers of A and never iterates a recursion with a potentially unstable state matrix, the paper asserts the algorithm is unconditionally stable: even if a structured approximation of A used inside the product has eigenvalues with magnitude above 1, the convolution remains bounded. This removes the constraint that approximate state matrices must be confined to the unit disk, and thereby widens the range of structured matrix representations that can accelerate the matrix-vector product.

Load-bearing premise

The argument rests on the assumption that the state matrix has all singular values below 1, which makes the cascade product converge; the paper's promise that eigenvalues above 1 are acceptable is not backed by that convergence proof.

Editorial extensions

If this is right

The eigenvalue restriction that confines standard SSM recursions to the open unit disk can be lifted in the cascade formulation, so structured approximations such as diagonal-plus-low-rank, partitioned low-rank, or wavelet non-standard forms can be used with less fear of instability.
The cost of the cascade is at most a factor of two above the direct recursion, and the algorithm is parallel-friendly since each level's updates act on disjoint index sets.
The approximation accuracy is controlled by N: the truncation error decays like γ^{2^{N+1}}/(1-γ), so any user-selected relative error can be met by choosing N of order log_2 log(1/ε).
For long sequences, the cascade's L + L/2 + L/4 + ... pattern of matrix-vector products amounts to about 2L operations, independent of the degree of the approximating polynomial.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

The unconditional-stability claim as stated in the abstract goes beyond what Lemma 1 proves, since the error bound requires the largest singular value γ to stay below 1; for an approximate eigenvalue above 1 the tail (z^{-1}A)^{2^{N+1}} grows, so a separate error analysis would be needed to support the claim.
For normal state matrices the cascade identity reduces to the scalar geometric identity on each eigenvalue, so stability depends only on the spectral radius; for highly non-normal matrices the singular-value condition is the binding constraint, and the cascade may be preferable to diagonalization when eigenvectors are ill-conditioned.
The cascade's reuse pattern suggests natural extensions to quantized or randomized matrix-vector products, since stability no longer couples to the spectrum of the approximate state matrix.
A direct numerical test on scalar systems with a = 1.01 would separate the paper's stability claim from its accuracy claim; the paper's example stays within γ<1 and so does not exercise the advertised beyond-unit-circle regime.

Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Desk Editor's Note

read the letter

Greg's paper is a short note that rewrites his earlier FIR-cascade approximation of IIR filters in state-space language. The core construction is correct for the stable case: when the largest singular value γ of A is below 1, Lemma 1's error bound holds and the cascade update rules in Section 3 are a clean way to apply the truncated polynomial. That part is not new—it is the same product factorization from his own [6,7]—but it is clearly presented.

The abstract, however, makes a stronger claim: the algorithm is unconditionally stable and 'it is not necessary to assure that an approximate state matrix has all eigenvalues with absolute values strictly less than 1.' That claim is false. Lemma 1 requires γ<1. For a scalar eigenvalue λ>1, the truncated product grows like λ^{2^{N+1}}, so no N gives a fixed accuracy; the tail term cannot be small. The stress-test counterexample is right. The paper does not offer a separate argument that covers that case.

There is also a complexity error. The update rules apply A^{2^{n−1}} to L−2^{n−1} columns at stage n. Summing over n gives O(L log L) matrix-vector products, not the promised 2L. Section 6's own description—A to every entry, A^2 to every second, A^4 to every fourth—confirms the log factor.

The paper has no experiments or code; the only example is a brief Mathematica check. The references are appropriate and include the author's prior work. So the stable-case material is fine, but the advertised advantage for unstable approximations is the main point, and it fails. A revised version limited to the stable case might be publishable as a short note, but the current claims cannot stand.

My recommendation: this does not deserve referee time as is. Desk reject or return for major revision with the unsupported claims removed.

Referee Report

2 major / 4 minor

Summary. The paper proposes a cascade algorithm for applying the transfer function H(z) = C(I_m - z^{-1}A)^{-1}B + D of a linear time-invariant system as a convolution in the time domain. The method replaces the inverse (I_m - z^{-1}A)^{-1} by the finite product \prod_{n=0}^N (I + z^{-2^n} A^{2^n}), whose time-domain realization updates a sequence by repeatedly applying the powers A, A^2, A^4, ... . The abstract and Section 6 claim that the resulting algorithm is unconditionally stable, that it requires no more than 2L matrix-vector multiplications to compute L states, and that it remains accurate even when an approximate state matrix has some eigenvalues with absolute value greater than 1. Section 4 reports a numerical example with a discretized HiPPO matrix whose eigenvalues are all below 1.

Significance. The underlying algebraic identity is elegant and correct in the stable case: for a matrix A with spectral radius below 1, the telescoping product gives a uniform polynomial approximation of (I - z^{-1}A)^{-1} on the unit circle, and Lemma 1 quantifies this approximation in terms of the largest singular value \gamma < 1. If the advertised complexity and stability properties were true, the method would be a useful alternative to direct recurrence for long-range state-space models. However, the central advertised property -- accuracy with eigenvalues exceeding 1 -- is false, and the complexity count in the abstract and Section 6 is inconsistent with the algorithm's own update rules. The paper contains no code, no reproducible experiments, and no numerical test of the claimed eigenvalue-above-1 regime, so the positive aspects are limited to the derivation of Lemma 1 and the cascade identity itself.

major comments (2)

[Section 2, Lemma 1; abstract; Section 6] The central claim of unconditional accuracy with eigenvalues above 1 is unsupported and, in fact, false. Lemma 1 assumes that the largest singular value \gamma of A is strictly less than 1, and its tail bound \gamma^{2^{N+1}}/(1-\gamma) exists only in that regime. The abstract and Section 6, however, state that within the desired accuracy some eigenvalues of an approximate state matrix may have absolute value exceeding 1 and that the cascade remains an unconditionally stable approximation of H(z). For a scalar eigenvalue \lambda > 1 the telescoping identity gives H_N(1) = \sum_{k=0}^{2^{N+1}-1} \lambda^k = (\lambda^{2^{N+1}}-1)/(\lambda-1), whereas H(1) = 1/(1-\lambda). With \lambda = 1.01, H(1) = -100 while H_N(1) grows like 100 \cdot 1.01^{2^{N+1}}, so no choice of N yields any fixed accuracy. Thus the claimed benefit of allowing structured approximations with eigenvalues above 1 fails exactly where Lemma 1's hypothesis is violated, and this is a load-bearing error for the paper's main message.
[Section 3, Eqs. (3.2)-(3.4); abstract; Section 6] The claimed arithmetic complexity is inconsistent with the stated algorithm. At stage n the update rule applies A^{2^{n-1}} to L - 2^{n-1} columns, so the total number of matrix-vector multiplications is \sum_{n} (L - 2^{n-1}) over the relevant range, which is O(L log L), not O(L) and not the abstract's 'no more than 2L'. For example, for L = 8 the numbers are 7 + 6 + 4 = 17 matrix-vector multiplications, already exceeding 16. Section 6's description that 'A is applied to every entry, A^2 to every second entry, A^4 to every fourth entry, etc.' does not match Eqs. (3.3)-(3.4), where A^2 is applied to L-2 columns and A^4 to L-4 columns. The complexity claims in the abstract and Section 3 consequently need to be corrected, and the claimed 'factor not exceeding two' advantage over the standard recursion is not supported.

minor comments (4)

[Eq. (3.1)] Equation (3.1) writes the direct convolution as y_\ell = \sum_{k=0}^\ell C A^{\ell-k} B u_k + D u_k; the last term should be D u_\ell, not D u_k.
[Section 2, after Eq. (2.2)] The sentence that applying the matrix polynomial 'as a cascade in time domain requires only N+1 matrix-vector multiplications' is at best ambiguous and, taken literally, is inconsistent with the algorithm in Section 3, which requires O(L log L) such multiplications for L outputs.
[Section 4] The numerical example only treats a stable discretized HiPPO matrix whose eigenvalues are all below 1; it provides no evidence for the advertised scenario of an approximate state matrix with eigenvalues greater than 1 in modulus.
[Section 5] The structured approximation proposals (partitioned low-rank, wavelet non-standard form, multiwavelets) are described only as possibilities; no concrete algorithm, error analysis, or numerical demonstration is provided for any of them.

Circularity Check

0 steps flagged · score 0.0 of 10

No circularity: the cascade convolution is a direct application of a standard matrix product identity; the advertised robustness to eigenvalues above 1 is a mathematical gap, not a circular reduction.

full rationale

The paper's central derivation is the infinite product identity (I - z^{-1}A)^{-1} = ∏_{n=0}^∞ (I + (z^{-1}A)^{2^n}), which is used to define H_N(z) as a truncated matrix polynomial and then implemented as a cascade in the time domain. This identity and the tail bound in Lemma 1 are standard Neumann-series algebra; Lemma 1 is stated under the explicit hypothesis that the largest singular value γ of A is less than 1, and its proof is said to be by induction, with references [6,7] cited only for similarity of technique. No fitted parameters, no data-calibrated constants, and no statistically forced predictions enter the derivation, so the central claim does not reduce to its own inputs by construction. The self-citations to the author's earlier FIR-filter papers are background for the cascade approximation method, not the load-bearing justification of the product identity. The advertised claim that eigenvalues may exceed 1 is not a circular step; it is an unsupported extension beyond Lemma 1's hypothesis, and the stated complexity estimate is an arithmetic inconsistency, both of which are correctness risks rather than circularity. Therefore no circular step is present.

Assumptions & free parameters 0 free parameters · 4 assumptions · 0 invented entities

No fitted constants appear. The derivation relies on a classical matrix series and on a domain assumption on the spectrum of A. The paper's extension to eigenvalues above 1 is an ad hoc assumption without proof, and the complexity bound is not derived from the stated algorithm.

assumptions (4)

domain assumption The state matrix A has all eigenvalues strictly inside the unit circle.
Stated in Section 1 and required for the infinite product and geometric series to converge. The paper later claims this condition is not necessary.
standard math The factorization (I - z^{-1}A)^{-1} = product_{n>=0}(I + (z^{-1}A)^{2^n}) holds uniformly on |z|=1.
Follows from the matrix geometric series when the spectral radius is below 1. Used in Section 2 to define H_N.
standard math The truncation error is bounded by gamma^{2^{N+1}}/(1-gamma) ||C|| ||B|| for the largest singular value gamma.
Lemma 1 in Section 2. The bound is valid only when gamma < 1.
ad hoc to paper A structured approximation of A with some eigenvalues above 1 in modulus still yields powers that track those of A closely enough.
Asserted in the abstract and Section 5 without proof. For |lambda| > 1, the powers lambda^{2^n} grow, so the assumption cannot hold.

how reviews work

0 comments

Cite this review

Pith. "Pith review of Fast convolution algorithm for state space models." pith.science (2026). https://pith.science/paper/7IAITW6F

@misc{pith2026241117729,
  author       = {Pith},
  title        = {Pith review of: Fast convolution algorithm for state space models},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/7IAITW6F}},
  note         = {Machine review of arXiv:2411.17729}
}

abstract

We present an unconditionally stable algorithm for applying matrix transfer function of a linear time invariant system (LTI) in time domain. The state matrix of an LTI system used for modeling long range dependencies in state space models (SSMs) has eigenvalues close to $1$. The standard recursion defining LTI system becomes unstable if the $m\times m$ state matrix has just one eigenvalue with absolute value even slightly greater than 1. This may occur when approximating a state matrix by a structured matrix to reduce the cost of matrix-vector multiplication from $\mathcal{O}\left(m^{2}\right)$ to $\mathcal{O}\left(m\right)$ or $\mathcal{O}\left(m\log m\right).$ We introduce an unconditionally stable algorithm that uses an approximation of the rational transfer function in the z-domain by a matrix polynomial of degree $2^{N+1}-1$, where $N$ is chosen to achieve any user-selected accuracy. Using a cascade implementation in time domain, applying such transfer function to compute $L$ states requires no more than $2L$ matrix-vector multiplications (whereas the standard recursion requires $L$ matrix-vector multiplications). However, using unconditionally stable algorithm, it is not necessary to assure that an approximate state matrix has all eigenvalues with absolute values strictly less than 1 i.e., within the desired accuracy, the absolute value of some eigenvalues may possibly exceed $1$. Consequently, this algorithm allows one to use a wider variety of structured approximations to reduce the cost of matrix-vector multiplication and we briefly describe several of them to be used for this purpose.

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

12 extracted references · 8 canonical work pages

[1]

Efﬁciently modeling long seque nces with structured state spaces,

A. Gu, K. Goel, and C. Ré, “Efﬁciently modeling long seque nces with structured state spaces,” arXiv preprint arXiv:2111.00396, 2021

arXiv 2021
[2]

Gu, Modeling Sequences with Structured State Spaces

A. Gu, Modeling Sequences with Structured State Spaces . Ph.D. Thesis, Stanford University, 2023

work page 2023
[3]

Mamba: Linear-time sequence modeling w ith selective state spaces,

A. Gu and T. Dao, “Mamba: Linear-time sequence modeling w ith selective state spaces,” arXiv preprint arXiv:2312.00752 , 2024

arXiv 2024
[4]

Robustifying state-space models for long sequence s via approximate diagonalization,

A. Y u, A. Nigmetov, D. Morozov, M. W. Mahoney, and N. B. Eri chson, “Robustifying state-space models for long sequence s via approximate diagonalization,” arXiv preprint arXiv:2310.01698, 2023

arXiv 2023
[5]

Simple hardware-efﬁcient long convolutions for sequence modeling,

D. Fu, E. L. Epstein, E. Nguyen, A. W. Thomas, M. Zhang, T. D ao, A. Rudra, and C. Ré, “Simple hardware-efﬁcient long convolutions for sequence modeling,” in International Conference on Machine Learning , pp. 10373–10391, PMLR, 2023

work page 2023
[6]

On factored FIR approximation of IIR ﬁlters ,

G. Beylkin, “On factored FIR approximation of IIR ﬁlters ,” Appl. Comput. Harmon. Anal. , vol. 2, no. 3, pp. 293–298, 1995

work page 1995
[7]

On the design of hi ghly accurate and efﬁcient IIR and FIR ﬁlters,

G. Beylkin, R. D. Lewis, and L. Monzón, “On the design of hi ghly accurate and efﬁcient IIR and FIR ﬁlters,” IEEE Trans. Signal Process., vol. 60, no. 8, pp. 4045–4054, 2012

work page 2012
[8]

How t o train your hippo: State space models with generalized orthogonal basis projections,

A. Gu, I. Johnson, A. Timalsina, A. Rudra, and C. Ré, “How t o train your hippo: State space models with generalized orthogonal basis projections,” arXiv preprint arXiv:2206.12037, 2022

arXiv 2022

Show all 12 references

[9]

Hippo: Recurr ent memory with optimal polynomial projections,

A. Gu, T. Dao, S. Ermon, A. Rudra, and C. Ré, “Hippo: Recurr ent memory with optimal polynomial projections,” Advances in neural information processing systems , vol. 33, pp. 1474–1487, 2020

2020
[10]

Wave propagation using bas es for bandlimited functions,

G. Beylkin and K. Sandberg, “Wave propagation using bas es for bandlimited functions,” W ave Motion, vol. 41, no. 3, pp. 263– 291, 2005

2005
[11]

Fast wavelet tr ansforms and numerical algorithms, I,

G. Beylkin, R. Coifman, and V . Rokhlin, “Fast wavelet tr ansforms and numerical algorithms, I,” Comm. Pure Appl. Math. , vol. 44, no. 2, pp. 141–183, 1991. Yale Univ. Technical Repor t Y ALEU/DCS/RR-696, August 1989. 5 FAST CONVOLUTION ALGORITHM FOR STA TE SPACE MODELS. GREGO...

1991
[12]

A class of bases in L2 for the sparse representation of integral operators,

B. Alpert, “A class of bases in L2 for the sparse representation of integral operators,” SIAM J. Math. Anal , vol. 24, no. 1, pp. 246–262, 1993. DEPARTMENT OF APPLIED MATHEMATICS , U NIVERSITY OF COLORADO AT BOULDER , UCB 526, B OULDER , CO 80309- 0526 6

1993

Pith tools

Reviewed August 12, 2026 · model on record in the stance chip above.

[1] [1]

Efﬁciently modeling long seque nces with structured state spaces,

A. Gu, K. Goel, and C. Ré, “Efﬁciently modeling long seque nces with structured state spaces,” arXiv preprint arXiv:2111.00396, 2021

arXiv 2021

[2] [2]

Gu, Modeling Sequences with Structured State Spaces

A. Gu, Modeling Sequences with Structured State Spaces . Ph.D. Thesis, Stanford University, 2023

work page 2023

[3] [3]

Mamba: Linear-time sequence modeling w ith selective state spaces,

A. Gu and T. Dao, “Mamba: Linear-time sequence modeling w ith selective state spaces,” arXiv preprint arXiv:2312.00752 , 2024

arXiv 2024

[4] [4]

Robustifying state-space models for long sequence s via approximate diagonalization,

A. Y u, A. Nigmetov, D. Morozov, M. W. Mahoney, and N. B. Eri chson, “Robustifying state-space models for long sequence s via approximate diagonalization,” arXiv preprint arXiv:2310.01698, 2023

arXiv 2023

[5] [5]

Simple hardware-efﬁcient long convolutions for sequence modeling,

D. Fu, E. L. Epstein, E. Nguyen, A. W. Thomas, M. Zhang, T. D ao, A. Rudra, and C. Ré, “Simple hardware-efﬁcient long convolutions for sequence modeling,” in International Conference on Machine Learning , pp. 10373–10391, PMLR, 2023

work page 2023

[6] [6]

On factored FIR approximation of IIR ﬁlters ,

G. Beylkin, “On factored FIR approximation of IIR ﬁlters ,” Appl. Comput. Harmon. Anal. , vol. 2, no. 3, pp. 293–298, 1995

work page 1995

[7] [7]

On the design of hi ghly accurate and efﬁcient IIR and FIR ﬁlters,

G. Beylkin, R. D. Lewis, and L. Monzón, “On the design of hi ghly accurate and efﬁcient IIR and FIR ﬁlters,” IEEE Trans. Signal Process., vol. 60, no. 8, pp. 4045–4054, 2012

work page 2012

[8] [8]

How t o train your hippo: State space models with generalized orthogonal basis projections,

A. Gu, I. Johnson, A. Timalsina, A. Rudra, and C. Ré, “How t o train your hippo: State space models with generalized orthogonal basis projections,” arXiv preprint arXiv:2206.12037, 2022

arXiv 2022