Levinson-Durbin algorithm: O(N²) solver for symmetric Toeplitz systems. More...
#include "config.h"#include "utilities/expected.h"#include "matrix_traits.h"#include "fused/fused_vector.h"#include "math/math_utils.h"
Go to the source code of this file.
Namespaces | |
| namespace | matrix_algorithms |
Functions | |
| template<typename T , my_size_t N> | |
| Expected< FusedVector< T, N >, MatrixStatus > | matrix_algorithms::levinson_durbin (const FusedVector< T, N > &r, const FusedVector< T, N > &b) |
| Solve a symmetric Toeplitz system Tx = b using Levinson-Durbin. | |
Detailed Description
Levinson-Durbin algorithm: O(N²) solver for symmetric Toeplitz systems.
A symmetric Toeplitz matrix is fully defined by its first row: T(i,j) = r(|i−j|)
[ r(0) r(1) r(2) ... r(N-1) ] [ r(1) r(0) r(1) ... r(N-2) ] [ r(2) r(1) r(0) ... r(N-3) ] [ ... ... ... ... ... ] [r(N-1)r(N-2) ... r(0) ]
The Levinson-Durbin recursion solves Tx = b in O(N²) by exploiting the Toeplitz structure. It also produces the reflection coefficients (partial correlation coefficients) as a byproduct.
ALGORITHM
The algorithm proceeds in two stages:
Stage 1 — Levinson recursion (solve T·a = −r for prediction coefficients): Initialize: a(0) = −r(1)/r(0), err = r(0)·(1 − a(0)²) For m = 1 … N−2: λ = −(r(m+1) + Σ_{k=0}^{m-1} a(k)·r(m−k)) / err Update a(k) using λ and reversed previous a err *= (1 − λ²)
Stage 2 — Solve Tx = b using the Levinson solution: Uses the forward and backward vectors from Stage 1 to iteratively build the solution.
Complexity: O(N²) multiply-adds, O(N) storage.
USE CASES
- Autocorrelation-based filtering (Wiener filter, linear prediction)
- Autoregressive (AR) model estimation
- Spectral estimation (Burg's method)
- Any system where the matrix has constant diagonals
FAILURE MODES
- MatrixStatus::Singular — r(0) ≈ 0 (zero autocorrelation)
- MatrixStatus::NotConverged — prediction error drops to zero or negative (matrix is not positive definite)
