Levinson recursion (Levinson-Durbin algorithm)
(linear_algebra_matrix/levinson.hpp)

• View this file on GitHub
• Last update: 2021-06-13 19:08:25+09:00
• Include: #include "linear_algebra_matrix/levinson.hpp"
• Link: https://en.wikipedia.org/wiki/Levinson_recursion

Solve Toeplitz system $\mathbf{T} \mathbf{x} = \mathbf{y}$, where $\mathbf{T}$ is Topelitz and every its principal submatrix is nonsingular.

Complexity

$O(N^2)$, where $N$ is the dimension of $\mathbf{T}$

References

• [1] N. Levinson, “The Wiener RMS error criterion in filter design and prediction,” J. Math. and Phys., 25, 261-298, 1947.
• [2] J. Durbin, “The fitting of time series models,”, Rev. Inst. Int. Stat. 28, 233-243, 1960.

Code

#pragma once
#include <cassert>
#include <vector>

// CUT begin
// solve Mx = y based on Levinson recursion
//
// M = [
// [ t_(N-1) t_(N-2) ... t_1 t_0     ]
// [ t_N     t_(N-1) ... t_2 t_1     ]
// [ ...                             ]
// [ t_(N*2-2)       ...     t_(N-1) ]]
//
// Complexity: O(N^2)
// Reference: https://en.wikipedia.org/wiki/Levinson_recursion
// return: - x = N-dim vector (if all principal minors of M are nonsingular)
//         - empty vector     (otherwise)
template <typename T>
std::vector<T> levinson_recursion(const std::vector<T> &toeplitz, const std::vector<T> &y) {
    const int N = y.size();
    assert(int(toeplitz.size()) == N * 2 - 1);
    std::vector<T> x(N);
    std::vector<T> fn(1), en(1);
    if (toeplitz[N - 1] == 0) { return {}; }
    fn[0] = en[0] = toeplitz[N - 1].inv();
    x[0] = fn[0] * y[0];
    for (int n = 1; n < N; n++) {
        T ef = 0, eb = 0;
        for (int i = 0; i < n; i++) {
            ef += fn[i] * toeplitz[N - 1 + n - i];
            eb += en[i] * toeplitz[N - 2 - i];
        }
        if (ef * eb == 1) { return {}; }
        auto nxtfn = fn, nxten = en;
        nxtfn.emplace_back(0), nxten.insert(nxten.begin(), 0);
        for (int i = 0; i < n; i++) {
            nxtfn[i + 1] -= ef * en[i];
            nxten[i] -= eb * fn[i];
        }
        T c = (1 - ef * eb).inv();
        for (int i = 0; i < n + 1; i++) { nxtfn[i] *= c, nxten[i] *= c; }
        fn = nxtfn, en = nxten;

        T z = y[n];
        for (int i = 0; i < n; i++) z -= x[i] * toeplitz[N - 1 + n - i];
        for (int i = 0; i < n + 1; i++) x[i] += en[i] * z;
    }
    return x;
}

#line 2 "linear_algebra_matrix/levinson.hpp"
#include <cassert>
#include <vector>

// CUT begin
// solve Mx = y based on Levinson recursion
//
// M = [
// [ t_(N-1) t_(N-2) ... t_1 t_0     ]
// [ t_N     t_(N-1) ... t_2 t_1     ]
// [ ...                             ]
// [ t_(N*2-2)       ...     t_(N-1) ]]
//
// Complexity: O(N^2)
// Reference: https://en.wikipedia.org/wiki/Levinson_recursion
// return: - x = N-dim vector (if all principal minors of M are nonsingular)
//         - empty vector     (otherwise)
template <typename T>
std::vector<T> levinson_recursion(const std::vector<T> &toeplitz, const std::vector<T> &y) {
    const int N = y.size();
    assert(int(toeplitz.size()) == N * 2 - 1);
    std::vector<T> x(N);
    std::vector<T> fn(1), en(1);
    if (toeplitz[N - 1] == 0) { return {}; }
    fn[0] = en[0] = toeplitz[N - 1].inv();
    x[0] = fn[0] * y[0];
    for (int n = 1; n < N; n++) {
        T ef = 0, eb = 0;
        for (int i = 0; i < n; i++) {
            ef += fn[i] * toeplitz[N - 1 + n - i];
            eb += en[i] * toeplitz[N - 2 - i];
        }
        if (ef * eb == 1) { return {}; }
        auto nxtfn = fn, nxten = en;
        nxtfn.emplace_back(0), nxten.insert(nxten.begin(), 0);
        for (int i = 0; i < n; i++) {
            nxtfn[i + 1] -= ef * en[i];
            nxten[i] -= eb * fn[i];
        }
        T c = (1 - ef * eb).inv();
        for (int i = 0; i < n + 1; i++) { nxtfn[i] *= c, nxten[i] *= c; }
        fn = nxtfn, en = nxten;

        T z = y[n];
        for (int i = 0; i < n; i++) z -= x[i] * toeplitz[N - 1 + n - i];
        for (int i = 0; i < n + 1; i++) x[i] += en[i] * z;
    }
    return x;
}

Back to top page