zilla_muf/structured/
toeplitz.rs

1use num_traits::Float;
2
3/// Causal Toeplitz matrix-vector product — the convolutional view of an
4/// SSM. Implicitly applies the lower-triangular matrix
5///   M[i][j] = kernel[i - j]  for j <= i,  0 otherwise
6/// to `x`, i.e. computes the causal convolution:
7///   y[i] = sum_{j=0}^{i} kernel[i - j] * x[j]
8///
9/// `kernel` and `x` must be the same length n; cost is O(n^2). This is
10/// the direct reference implementation — `fft_conv` computes the exact
11/// same quantity in O(n log n) and is tested against this function.
12pub fn toeplitz_matvec<T: Float>(kernel: &[T], x: &[T]) -> Vec<T> {
13	let n = x.len();
14	assert_eq!(kernel.len(), n, "kernel and x must be the same length");
15
16	// For each output i, sum kernel[i - j] * x[j] over j <= i. The
17	// `i - j` lag index is what makes this Toeplitz (entries are constant
18	// along each diagonal); the `j <= i` bound is what makes it causal
19	// (no future input can influence an earlier output).
20	(0..n)
21		.map(|i| {
22			(0..=i).fold(T::zero(), |acc, j| acc + kernel[i - j] * x[j])
23		})
24		.collect()
25}
26
27#[cfg(test)]
28mod tests {
29	use super::*;
30
31	/// Builds the dense n x n lower-triangular Toeplitz matrix directly
32	/// and does a naive matvec — an oracle independent of the
33	/// convolution-loop formulation used above.
34	fn dense_reference(kernel: &[f64], x: &[f64]) -> Vec<f64> {
35		let n = x.len();
36		let mut y = vec![0.0; n];
37		for i in 0..n {
38			for j in 0..=i {
39				y[i] += kernel[i - j] * x[j];
40			}
41		}
42		y
43	}
44
45	#[test]
46	fn matches_dense_reference() {
47		let n = 8;
48		let kernel: Vec<f64> = (0..n).map(|i| (i as f64 * 0.41).cos() * 0.8_f64.powi(i as i32)).collect();
49		let x: Vec<f64> = (0..n).map(|i| (i as f64 * 0.23).sin()).collect();
50
51		let expected = dense_reference(&kernel, &x);
52		let actual = toeplitz_matvec(&kernel, &x);
53		for (e, got) in expected.iter().zip(actual.iter()) {
54			assert!((e - got).abs() < 1e-9, "expected {e}, got {got}");
55		}
56	}
57
58	#[test]
59	fn impulse_kernel_passes_signal_through() {
60		// kernel = [1, 0, 0, ...] is the identity matrix -> y == x
61		let n = 5;
62		let mut kernel = vec![0.0; n];
63		kernel[0] = 1.0;
64		let x: Vec<f64> = (0..n).map(|i| i as f64 + 1.0).collect();
65		assert_eq!(toeplitz_matvec(&kernel, &x), x);
66	}
67}
zilla_muf/structured/toeplitz.rs

zilla_muf/structured/
toeplitz.rs
