From ea7d6f5aa12f420ec91792cbcdadd6b31c358b14 Mon Sep 17 00:00:00 2001 From: Sebastian Jiro Schlecht <45119263+SebastianJiroSchlecht@users.noreply.github.com> Date: Thu, 1 Oct 2026 13:28:58 +0200 Subject: [PATCH 1/2] Add FDN recurrence example: can a lossless FDN unmix itself? Runs small lossless 4-line FDNs for up to 1M samples in closed form (Schur diagonalisation of the dss_to_ss transition matrix) and tracks the participation ratio of the state. Compares Hadamard, Householder and Haar-random orthogonal feedback matrices across four delay sets. Co-Authored-By: Claude Opus 5.5 --- examples/example_fdn_recurrence.py | 240 +++++++++++++++++++++++++++++ 1 file changed, 240 insertions(+) create mode 100644 examples/example_fdn_recurrence.py diff --git a/examples/example_fdn_recurrence.py b/examples/example_fdn_recurrence.py new file mode 100644 index 00000000..b84994c5 --- /dev/null +++ b/examples/example_fdn_recurrence.py @@ -0,0 +1,240 @@ +# gallery_category: Feedback Matrices +# gallery_description: Run a tiny lossless FDN for millions of samples to see whether it unmixes back into a sparse response, and how the feedback matrix decides that. + +import marimo + +__generated_with = "0.23.9" +app = marimo.App() + + +@app.cell +def _(): + import marimo as mo + + return (mo,) + + +@app.cell(hide_code=True) +def _(mo): + mo.md(r""" + # Can an FDN unmix itself? + + A lossless FDN is an orthogonal map on its $M = \sum_i m_i$ delay-line samples: energy is never lost, only shuffled. Poincaré recurrence then says the state must come back arbitrarily close to where it started, infinitely often. If it does, the impulse response should briefly collapse from dense noise back into a few sparse echoes. + + The catch is *how long* that takes. The state moves on a torus spanned by the eigenphases of the $M \times M$ state-transition matrix. Poles that are exact roots of unity are periodic and come back on their own; only the rationally independent rest have to line up by chance, and the waiting time grows roughly like $(1/\varepsilon)^{d}$ in the number $d$ of those free modes. So the question becomes: how many free modes does a given feedback matrix leave? + + Four short delays keep $M = 26$, small enough to diagonalise exactly and evaluate millions of samples in closed form. + """) + return + + +@app.cell +def _(): + import matplotlib.pyplot as plt + import numpy as np + from scipy.linalg import hadamard, schur + + import pyFDN + + delays = np.array([3, 5, 7, 11]) + n = len(delays) + horizon = 1_000_000 + return delays, hadamard, horizon, n, np, plt, pyFDN, schur + + +@app.cell(hide_code=True) +def _(mo): + mo.md(r""" + ## Running it for a million samples + + `pyFDN.dss_to_ss` expands the FDN into its $M$-state transition matrix $\mathbf{T}$. Because $\mathbf{T}$ is orthogonal, its complex Schur form is diagonal: $\mathbf{T} = \mathbf{Z}\boldsymbol{\Lambda}\mathbf{Z}^H$, so the state after $k$ samples is $\mathbf{x}(k) = \mathbf{Z}\boldsymbol{\Lambda}^k\mathbf{Z}^H\mathbf{x}(0)$ — no sample-by-sample loop, and no accumulated rounding error. + + Sparsity of the state is measured by its participation ratio $\|\mathbf{x}\|_2^4 / \|\mathbf{x}\|_4^4$, the effective number of occupied delay-line taps. The impulse enters all four lines equally, so it starts at 4; a fully mixed, Gaussian-like state sits near $M/3 \approx 9$; a value near 1 means all the energy is back in a single tap. A pole counts as *periodic* when $\lambda^k = 1$ for some $k \le 5000$. + """) + return + + +@app.cell +def _(n, np, pyFDN, schur): + def recurrence(delays, A, horizon, chunk=100_000): + """Participation ratio, output and periodic poles of a lossless FDN.""" + T, b, c, _ = pyFDN.dss_to_ss( + delays, A, np.ones((n, 1)) / np.sqrt(n), np.ones((1, n)), np.zeros((1, 1)) + ) + S, Z = schur(T, output="complex") + lam = np.diag(S) + modes = Z.conj().T @ b[:, 0] + pr = np.empty(horizon) + y = np.empty(horizon) + for start in range(0, horizon, chunk): + k = np.arange(start, min(start + chunk, horizon)) + x = (Z @ (modes[:, None] * lam[:, None] ** k)).real + pr[k] = (x**2).sum(0) ** 2 / (x**4).sum(0) + y[k] = c[0] @ x + orders = np.arange(1, 5001) + hits = np.abs(lam[:, None] ** orders - 1) < 1e-8 + periodic = hits.any(1) + period = int(np.lcm.reduce(orders[hits[periodic].argmax(1)], initial=1)) + return pr, y, int(periodic.sum()), period + + return (recurrence,) + + +@app.cell(hide_code=True) +def _(mo): + mo.md(r""" + ## Hadamard against random orthogonal + + The $4 \times 4$ Hadamard matrix against a batch of Haar-random orthogonal matrices, all on the same delays. The Householder reflection $\mathbf{I} - \tfrac{1}{2}\mathbf{1}\mathbf{1}^\top$ is a useful control: like the Hadamard matrix (which it equals up to signs and a row permutation) its eigenvalues are only $\pm 1$, so if eigenvalue structure of the matrix alone were what mattered, the two should behave alike. + """) + return + + +@app.cell +def _(delays, hadamard, horizon, n, np, pyFDN, recurrence): + rng = np.random.default_rng(0) + matrices = { + "Hadamard": hadamard(n) / np.sqrt(n), + "Householder": np.eye(n) - 2 / n, + } + matrices |= {f"random {i}": pyFDN.random_orthogonal(n, rng) for i in range(12)} + results = {name: recurrence(delays, A, horizon) for name, A in matrices.items()} + return matrices, results + + +@app.cell +def _(delays, mo, np, results): + _settle = 1000 # skip the initial transient, where the state is still sparse + mo.ui.table( + [ + { + "matrix": name, + "periodic poles": periodic, + "free poles": int(delays.sum()) - periodic, + "common period": period if periodic else None, + "min participation": round(float(pr[_settle:].min()), 2), + "at sample": int(np.argmin(pr[_settle:]) + _settle), + } + for name, (pr, _y, periodic, period) in results.items() + ], + selection=None, + ) + return + + +@app.cell(hide_code=True) +def _(mo): + mo.md(r""" + ## Over time + + The running minimum of the participation ratio in blocks of 2000 samples. Every curve falls below the starting value of 4 at some point — partial unmixing happens for any matrix — but only the Hadamard FDN gets close to a single tap, and the Householder FDN barely unmixes at all. + """) + return + + +@app.cell +def _(horizon, np, plt, results): + _block = 2000 + _t = np.arange(0, horizon, _block) + _block / 2 + _fig, _ax = plt.subplots(figsize=(9, 4)) + _highlight = {"Hadamard": "C3", "Householder": "C0"} + for _name, (_pr, *_) in results.items(): + _env = _pr.reshape(-1, _block).min(1) + if _name in _highlight: + _ax.plot(_t, _env, color=_highlight[_name], lw=1.2, label=_name, zorder=3) + else: + _ax.plot(_t, _env, color="0.6", lw=0.5) + _ax.plot([], [], color="0.6", lw=0.5, label="random orthogonal") + _ax.axhline(4, color="k", ls="--", lw=0.8, label="initial state") + _ax.set_xlabel("Time (samples)") + _ax.set_ylabel("Min participation ratio per block") + _ax.legend(loc="upper right") + _fig.tight_layout() + _fig.gca() + return + + +@app.cell(hide_code=True) +def _(mo): + mo.md(r""" + ## The response at its sparsest + + The output around the Hadamard FDN's sparsest moment, next to a random matrix at *its* sparsest. Sample 0 marks the minimum. + """) + return + + +@app.cell +def _(np, plt, results): + _fig, _axes = plt.subplots(2, 1, figsize=(9, 4.5), sharex=True) + for _ax, _name in zip(_axes, ["Hadamard", "random 0"], strict=True): + _pr, _y, *_ = results[_name] + _k = int(np.argmin(_pr[1000:]) + 1000) + _lags = np.arange(-60, 61) + _ax.stem(_lags, _y[_k + _lags], basefmt=" ") + _ax.set_title(f"{_name}: around sample {_k} (participation {_pr[_k]:.2f})") + _ax.set_ylabel("Output") + _axes[-1].set_xlabel("Samples relative to the sparsest moment") + _fig.tight_layout() + _fig.gca() + return + + +@app.cell(hide_code=True) +def _(mo): + mo.md(r""" + ## Other delays + +One delay set proves little, so the same comparison on three more, over a shorter horizon of 300 000 samples. Each row is one delay set: the minimum participation ratio reached by the Hadamard and Householder FDNs, and the best and median of eight random orthogonal matrices. + """) + return + + +@app.cell +def _(hadamard, matrices, mo, n, np, pyFDN, recurrence): + _rng = np.random.default_rng(1) + _rows = [] + for _delays in ([3, 5, 7, 11], [5, 7, 11, 13], [4, 6, 9, 13], [7, 11, 13, 17]): + _min = { + _name: recurrence(_delays, matrices[_name], 300_000)[0][1000:].min() + for _name in ("Hadamard", "Householder") + } + _rand = [ + recurrence(_delays, pyFDN.random_orthogonal(n, _rng), 300_000)[0][ + 1000: + ].min() + for _ in range(8) + ] + _rows.append( + { + "delays": str(_delays), + "M": sum(_delays), + "Hadamard": round(float(_min["Hadamard"]), 2), + "Householder": round(float(_min["Householder"]), 2), + "random best": round(float(min(_rand)), 2), + "random median": round(float(np.median(_rand)), 2), + } + ) + mo.ui.table(_rows, selection=None) + return + + +@app.cell(hide_code=True) +def _(mo): + mo.md(r""" + ## Why Hadamard + + The Hadamard FDN gets sparsest on every delay set, ahead of the best random matrix each time. The structure of the matrix clearly matters. Its eigenvalues are not the reason, though. The Householder matrix has the same $\pm 1$ eigenvalues and the same entry magnitudes, and it unmixes *less* than a typical random matrix. Only the sign pattern differs. + + Two observations point at the mechanism, without settling it: + + - **Focusing.** One sample before the Hadamard FDN's sparsest moment, almost all the energy sits in four taps of nearly equal magnitude reaching the matrix together. The Sylvester Hadamard matrix maps an equal-magnitude vector with the right signs (any of its columns) onto a single line, so a near-return of that pattern collapses into one echo. The Householder matrix only flips the sign of $\mathbf{1}$, its own eigenvector, so the same moment stays four echoes wide. + - **Periodic poles.** On the first delay set, 12 of the Hadamard FDN's 26 poles are exact roots of unity with a common period of 42 samples. They come back on their own, which leaves fewer free modes that have to line up by chance. That cannot be the whole story: on other delays the Hadamard FDN has few periodic poles and still wins. + + Two caveats keep this a curiosity rather than an audible effect. First, it scales badly: the minima creep up as $M$ grows, even over these few tens of states. A practical FDN has thousands of free modes, and recurrence time grows exponentially with their number. Second, a real FDN is lossy: by the time even this tiny network unmixes, an RT of a second has pushed the response hundreds of decibels below the start. + """) + return + + +if __name__ == "__main__": + app.run() From 11283d77d26fb2f0f87253756bfb44fd3e6d80c8 Mon Sep 17 00:00:00 2001 From: Sebastian Jiro Schlecht <45119263+SebastianJiroSchlecht@users.noreply.github.com> Date: Thu, 1 Oct 2026 13:39:34 +0200 Subject: [PATCH 2/2] FDN recurrence example: scale to M = 1000 with a block time-domain recursion Closed-form Schur evaluation is O(M^2) per sample, too slow at M = 1000. The block recursion over delay-line inputs is exact (checked against dss_to_impz) and computes the state participation ratio from sliding window sums. Adds a size sweep from M = 26 to 1000. Co-Authored-By: Claude Opus 5.5 --- examples/example_fdn_recurrence.py | 208 +++++++++++++++++------------ 1 file changed, 119 insertions(+), 89 deletions(-) diff --git a/examples/example_fdn_recurrence.py b/examples/example_fdn_recurrence.py index b84994c5..58edff10 100644 --- a/examples/example_fdn_recurrence.py +++ b/examples/example_fdn_recurrence.py @@ -1,5 +1,5 @@ # gallery_category: Feedback Matrices -# gallery_description: Run a tiny lossless FDN for millions of samples to see whether it unmixes back into a sparse response, and how the feedback matrix decides that. +# gallery_description: Run a lossless FDN for millions of samples to see whether it unmixes back into a sparse response, and how the feedback matrix and size decide that. import marimo @@ -21,9 +21,7 @@ def _(mo): A lossless FDN is an orthogonal map on its $M = \sum_i m_i$ delay-line samples: energy is never lost, only shuffled. Poincaré recurrence then says the state must come back arbitrarily close to where it started, infinitely often. If it does, the impulse response should briefly collapse from dense noise back into a few sparse echoes. - The catch is *how long* that takes. The state moves on a torus spanned by the eigenphases of the $M \times M$ state-transition matrix. Poles that are exact roots of unity are periodic and come back on their own; only the rationally independent rest have to line up by chance, and the waiting time grows roughly like $(1/\varepsilon)^{d}$ in the number $d$ of those free modes. So the question becomes: how many free modes does a given feedback matrix leave? - - Four short delays keep $M = 26$, small enough to diagonalise exactly and evaluate millions of samples in closed form. + The catch is *how long* that takes. The state moves quasi-periodically on a torus spanned by the eigenphases of the $M \times M$ state-transition matrix, and the waiting time for all of them to line up again grows exponentially with the number of independent modes. Four delay lines with $M = 1000$ states are long enough to watch the mixing happen, and short enough to run for a minute of audio. """) return @@ -32,60 +30,83 @@ def _(mo): def _(): import matplotlib.pyplot as plt import numpy as np - from scipy.linalg import hadamard, schur + from scipy.linalg import hadamard import pyFDN - delays = np.array([3, 5, 7, 11]) + delays = np.array([227, 241, 251, 281]) # M = 1000 n = len(delays) - horizon = 1_000_000 - return delays, hadamard, horizon, n, np, plt, pyFDN, schur + horizon = 2_000_000 + return delays, hadamard, horizon, n, np, plt, pyFDN @app.cell(hide_code=True) def _(mo): mo.md(r""" - ## Running it for a million samples + ## Running it for two million samples - `pyFDN.dss_to_ss` expands the FDN into its $M$-state transition matrix $\mathbf{T}$. Because $\mathbf{T}$ is orthogonal, its complex Schur form is diagonal: $\mathbf{T} = \mathbf{Z}\boldsymbol{\Lambda}\mathbf{Z}^H$, so the state after $k$ samples is $\mathbf{x}(k) = \mathbf{Z}\boldsymbol{\Lambda}^k\mathbf{Z}^H\mathbf{x}(0)$ — no sample-by-sample loop, and no accumulated rounding error. + The state of the FDN at time $k$ is the last $m_i$ samples written into each delay line $i$. Writing $s_i(k)$ for the signal entering line $i$, the lossless recursion is $\mathbf{s}(k) = \mathbf{A}\,[s_1(k - m_1), \dots, s_N(k - m_N)]^\top$, with the impulse $\mathbf{s}(0) = \mathbf{1}/\sqrt{N}$. No line reads anything younger than the shortest delay, so whole blocks of that length are computed with one matrix product. - Sparsity of the state is measured by its participation ratio $\|\mathbf{x}\|_2^4 / \|\mathbf{x}\|_4^4$, the effective number of occupied delay-line taps. The impulse enters all four lines equally, so it starts at 4; a fully mixed, Gaussian-like state sits near $M/3 \approx 9$; a value near 1 means all the energy is back in a single tap. A pole counts as *periodic* when $\lambda^k = 1$ for some $k \le 5000$. + Sparsity of the state is measured by its participation ratio $\|\mathbf{x}\|_2^4 / \|\mathbf{x}\|_4^4$, the effective number of occupied delay-line taps. Both norms are sums of $s_i^2$ and $s_i^4$ over a sliding window of length $m_i$ per line, which cumulative sums give at every sample. The impulse enters all four lines equally, so the ratio starts at 4; a fully mixed, Gaussian-like state sits near $M/3 \approx 333$; a value of 1 means all the energy is in a single tap. """) return @app.cell -def _(n, np, pyFDN, schur): - def recurrence(delays, A, horizon, chunk=100_000): - """Participation ratio, output and periodic poles of a lossless FDN.""" - T, b, c, _ = pyFDN.dss_to_ss( - delays, A, np.ones((n, 1)) / np.sqrt(n), np.ones((1, n)), np.zeros((1, 1)) +def _(np): + def recurrence(delays, A, horizon): + """Delay-line input signals and state participation ratio of a lossless FDN.""" + delays = np.asarray(delays) + n = len(delays) + pad = delays.max() + s = np.zeros((pad + horizon, n)) + s[pad] = 1 / np.sqrt(n) + for start in range(pad + 1, pad + horizon, delays.min()): + stop = min(start + delays.min(), pad + horizon) + heads = [s[start - m : stop - m, i] for i, m in enumerate(delays)] + s[start:stop] = np.stack(heads, 1) @ A.T + sums = [np.cumsum(np.vstack([np.zeros(n), s**p]), 0) for p in (2, 4)] + t = np.arange(pad, pad + horizon) + 1 # state after writing sample t - 1 + e2, e4 = ( + sum(c[t, i] - c[t - m, i] for i, m in enumerate(delays)) for c in sums ) - S, Z = schur(T, output="complex") - lam = np.diag(S) - modes = Z.conj().T @ b[:, 0] - pr = np.empty(horizon) - y = np.empty(horizon) - for start in range(0, horizon, chunk): - k = np.arange(start, min(start + chunk, horizon)) - x = (Z @ (modes[:, None] * lam[:, None] ** k)).real - pr[k] = (x**2).sum(0) ** 2 / (x**4).sum(0) - y[k] = c[0] @ x - orders = np.arange(1, 5001) - hits = np.abs(lam[:, None] ** orders - 1) < 1e-8 - periodic = hits.any(1) - period = int(np.lcm.reduce(orders[hits[periodic].argmax(1)], initial=1)) - return pr, y, int(periodic.sum()), period + return s[pad:], e2**2 / e4 return (recurrence,) +@app.cell(hide_code=True) +def _(mo): + mo.md(r""" + The output is the sum of the four delay-line outputs. It matches `pyFDN.dss_to_impz` over the first 5000 samples, which checks the block recursion. + """) + return + + +@app.cell +def _(delays, hadamard, n, np, pyFDN, recurrence): + def output(s): + return sum(np.r_[np.zeros(m), s[:-m, i]] for i, m in enumerate(delays)) + + _A = hadamard(n) / np.sqrt(n) + _ref = pyFDN.dss_to_impz( + delays, + _A, + np.ones((n, 1)) / np.sqrt(n), + np.ones((1, n)), + np.zeros((1, 1)), + 5000, + )[:, 0, 0] + assert np.allclose(output(recurrence(delays, _A, 5000)[0]), _ref) + return (output,) + + @app.cell(hide_code=True) def _(mo): mo.md(r""" ## Hadamard against random orthogonal - The $4 \times 4$ Hadamard matrix against a batch of Haar-random orthogonal matrices, all on the same delays. The Householder reflection $\mathbf{I} - \tfrac{1}{2}\mathbf{1}\mathbf{1}^\top$ is a useful control: like the Hadamard matrix (which it equals up to signs and a row permutation) its eigenvalues are only $\pm 1$, so if eigenvalue structure of the matrix alone were what mattered, the two should behave alike. + The $4 \times 4$ Hadamard matrix against a batch of Haar-random orthogonal matrices, all on the same delays. The Householder reflection $\mathbf{I} - \tfrac{1}{2}\mathbf{1}\mathbf{1}^\top$ is a useful control: it equals the Hadamard matrix up to signs and a row permutation, so it has the same $\pm 1$ eigenvalues and the same entry magnitudes. """) return @@ -103,19 +124,19 @@ def _(delays, hadamard, horizon, n, np, pyFDN, recurrence): @app.cell -def _(delays, mo, np, results): - _settle = 1000 # skip the initial transient, where the state is still sparse +def _(mo, np, results): + _settle = 100_000 # past the initial mixing mo.ui.table( [ { "matrix": name, - "periodic poles": periodic, - "free poles": int(delays.sum()) - periodic, - "common period": period if periodic else None, - "min participation": round(float(pr[_settle:].min()), 2), + "after 1000": round(float(pr[1000]), 1), + "after 5000": round(float(pr[5000]), 1), + "late median": round(float(np.median(pr[_settle:])), 1), + "late min": round(float(pr[_settle:].min()), 1), "at sample": int(np.argmin(pr[_settle:]) + _settle), } - for name, (pr, _y, periodic, period) in results.items() + for name, (_s, pr) in results.items() ], selection=None, ) @@ -125,30 +146,38 @@ def _(delays, mo, np, results): @app.cell(hide_code=True) def _(mo): mo.md(r""" - ## Over time + ## Mixing, then nothing - The running minimum of the participation ratio in blocks of 2000 samples. Every curve falls below the starting value of 4 at some point — partial unmixing happens for any matrix — but only the Hadamard FDN gets close to a single tap, and the Householder FDN barely unmixes at all. + Left: the first 50 000 samples at full resolution, as the participation ratio climbs from 4 towards the fully mixed plateau. Right: the minimum per block of 4000 samples over the whole run. It fluctuates below the plateau but never comes back down towards 4: over 42 seconds at 48 kHz, no matrix gets below 190 taps. + + The Hadamard FDN behaves differently from the start. It reaches its plateau fastest, within about 5000 samples, but that plateau sits lower: its median of about 300 taps is below every random matrix, so it never mixes fully. The Householder FDN mixes almost as fast and then sits right at $M/3$. """) return @app.cell -def _(horizon, np, plt, results): - _block = 2000 - _t = np.arange(0, horizon, _block) + _block / 2 - _fig, _ax = plt.subplots(figsize=(9, 4)) +def _(horizon, n, np, plt, results): _highlight = {"Hadamard": "C3", "Householder": "C0"} - for _name, (_pr, *_) in results.items(): - _env = _pr.reshape(-1, _block).min(1) + _block = 4000 + _fig, (_early, _late) = plt.subplots(1, 2, figsize=(11, 4), sharey=True) + for _name, (_s, _pr) in results.items(): + _style = {"color": _highlight.get(_name, "0.6"), "lw": 0.5} if _name in _highlight: - _ax.plot(_t, _env, color=_highlight[_name], lw=1.2, label=_name, zorder=3) - else: - _ax.plot(_t, _env, color="0.6", lw=0.5) - _ax.plot([], [], color="0.6", lw=0.5, label="random orthogonal") - _ax.axhline(4, color="k", ls="--", lw=0.8, label="initial state") - _ax.set_xlabel("Time (samples)") - _ax.set_ylabel("Min participation ratio per block") - _ax.legend(loc="upper right") + _style |= {"label": _name, "lw": 1.0, "zorder": 3} + _early.plot(_pr[:50_000], **_style) + _late.plot( + np.arange(0, horizon, _block) + _block / 2, + _pr.reshape(-1, _block).min(1), + **_style, + ) + _early.plot([], [], color="0.6", lw=0.5, label="random orthogonal") + for _ax in (_early, _late): + _ax.axhline(n, color="k", ls="--", lw=0.8) + _ax.set_xlabel("Time (samples)") + _ax.set_yscale("log") + _early.set_ylabel("Participation ratio (taps)") + _late.set_ylabel("Min per block") + _early.legend(loc="lower right") _fig.tight_layout() _fig.gca() return @@ -159,22 +188,24 @@ def _(mo): mo.md(r""" ## The response at its sparsest - The output around the Hadamard FDN's sparsest moment, next to a random matrix at *its* sparsest. Sample 0 marks the minimum. + The output of the Hadamard FDN at the start, and around its sparsest moment after mixing. Even that moment looks like dense noise. """) return @app.cell -def _(np, plt, results): +def _(np, output, plt, results): + _s, _pr = results["Hadamard"] + _y = output(_s) + _k = int(np.argmin(_pr[100_000:]) + 100_000) _fig, _axes = plt.subplots(2, 1, figsize=(9, 4.5), sharex=True) - for _ax, _name in zip(_axes, ["Hadamard", "random 0"], strict=True): - _pr, _y, *_ = results[_name] - _k = int(np.argmin(_pr[1000:]) + 1000) - _lags = np.arange(-60, 61) - _ax.stem(_lags, _y[_k + _lags], basefmt=" ") - _ax.set_title(f"{_name}: around sample {_k} (participation {_pr[_k]:.2f})") + _axes[0].plot(_y[:2000], lw=0.7) + _axes[0].set_title("Start of the response") + _axes[1].plot(_y[_k - 1000 : _k + 1000], lw=0.7) + _axes[1].set_title(f"Around sample {_k} (participation {_pr[_k]:.0f})") + for _ax in _axes: _ax.set_ylabel("Output") - _axes[-1].set_xlabel("Samples relative to the sparsest moment") + _axes[-1].set_xlabel("Samples in window") _fig.tight_layout() _fig.gca() return @@ -183,36 +214,38 @@ def _(np, plt, results): @app.cell(hide_code=True) def _(mo): mo.md(r""" - ## Other delays + ## From 26 to 1000 states -One delay set proves little, so the same comparison on three more, over a shorter horizon of 300 000 samples. Each row is one delay set: the minimum participation ratio reached by the Hadamard and Householder FDNs, and the best and median of eight random orthogonal matrices. + Small networks do unmix. The same comparison over a range of sizes, one million samples each, shows the deepest dip after mixing as a fraction of the fully mixed $M/3$: Hadamard, Householder, and the best and median of eight random orthogonal matrices. """) return @app.cell -def _(hadamard, matrices, mo, n, np, pyFDN, recurrence): +def _(matrices, mo, n, np, pyFDN, recurrence): _rng = np.random.default_rng(1) _rows = [] - for _delays in ([3, 5, 7, 11], [5, 7, 11, 13], [4, 6, 9, 13], [7, 11, 13, 17]): - _min = { - _name: recurrence(_delays, matrices[_name], 300_000)[0][1000:].min() - for _name in ("Hadamard", "Householder") - } - _rand = [ - recurrence(_delays, pyFDN.random_orthogonal(n, _rng), 300_000)[0][ - 1000: - ].min() - for _ in range(8) - ] + for _delays in ( + [3, 5, 7, 11], + [7, 11, 13, 17], + [17, 23, 29, 31], + [53, 59, 67, 71], + [227, 241, 251, 281], + ): + + def _dip(A, d=_delays): + pr = recurrence(d, A, 1_000_000)[1] + return pr[100 * sum(d) :].min() / (sum(d) / 3) + + _rand = [_dip(pyFDN.random_orthogonal(n, _rng)) for _ in range(8)] _rows.append( { "delays": str(_delays), "M": sum(_delays), - "Hadamard": round(float(_min["Hadamard"]), 2), - "Householder": round(float(_min["Householder"]), 2), - "random best": round(float(min(_rand)), 2), - "random median": round(float(np.median(_rand)), 2), + "Hadamard": round(float(_dip(matrices["Hadamard"])), 3), + "Householder": round(float(_dip(matrices["Householder"])), 3), + "random best": round(float(min(_rand)), 3), + "random median": round(float(np.median(_rand)), 3), } ) mo.ui.table(_rows, selection=None) @@ -222,16 +255,13 @@ def _(hadamard, matrices, mo, n, np, pyFDN, recurrence): @app.cell(hide_code=True) def _(mo): mo.md(r""" - ## Why Hadamard - - The Hadamard FDN gets sparsest on every delay set, ahead of the best random matrix each time. The structure of the matrix clearly matters. Its eigenvalues are not the reason, though. The Householder matrix has the same $\pm 1$ eigenvalues and the same entry magnitudes, and it unmixes *less* than a typical random matrix. Only the sign pattern differs. + ## What decides it - Two observations point at the mechanism, without settling it: + **Size.** With 26 states, the Hadamard FDN gets back to almost a single tap within a million samples (12% of $M/3$, about one tap). Every added state is another phase that has to line up, so the deepest dip climbs steadily towards the fully mixed value: by $M = 1000$, no matrix gets below about 58% of it. A practical FDN, with several thousand states, will not unmix in any realistic listening time. Losses make the point moot anyway: an RT of a second leaves the response hundreds of decibels down long before then. - - **Focusing.** One sample before the Hadamard FDN's sparsest moment, almost all the energy sits in four taps of nearly equal magnitude reaching the matrix together. The Sylvester Hadamard matrix maps an equal-magnitude vector with the right signs (any of its columns) onto a single line, so a near-return of that pattern collapses into one echo. The Householder matrix only flips the sign of $\mathbf{1}$, its own eigenvector, so the same moment stays four echoes wide. - - **Periodic poles.** On the first delay set, 12 of the Hadamard FDN's 26 poles are exact roots of unity with a common period of 42 samples. They come back on their own, which leaves fewer free modes that have to line up by chance. That cannot be the whole story: on other delays the Hadamard FDN has few periodic poles and still wins. + **Matrix.** At every size the Hadamard FDN dips deepest. The Householder FDN, with identical eigenvalues and entry magnitudes, does no better than a typical random matrix and usually worse. So the eigenvalues of $\mathbf{A}$ are not the reason; its sign pattern is. Poles that are exact roots of unity are not the reason either: at $M = 1000$ the Hadamard FDN has only the trivial ones at $\pm 1$, like a random matrix. - Two caveats keep this a curiosity rather than an audible effect. First, it scales badly: the minima creep up as $M$ grows, even over these few tens of states. A practical FDN has thousands of free modes, and recurrence time grows exponentially with their number. Second, a real FDN is lossy: by the time even this tiny network unmixes, an RT of a second has pushed the response hundreds of decibels below the start. + A plausible mechanism, seen at $M = 26$: just before the sparsest moment, the energy sits in four nearly equal taps reaching the matrix together, and the Sylvester Hadamard matrix maps any of its $\pm 1$ columns onto a single line. The Householder matrix only flips the sign of $\mathbf{1}$, its own eigenvector, so that moment stays four echoes wide. The same focusing would explain why the large Hadamard FDN never quite mixes fully. """) return