yule_walker_solver.cc

#include <vtestbed/base/for_loop.hh>
#include <vtestbed/config/real_type.hh>
#include <vtestbed/core/debug.hh>
#include <vtestbed/core/yule_walker_solver.hh>
#include <vtestbed/linalg/linear_algebra.hh>
#include <vtestbed/profile/profile.hh>

#include <algorithm>
#include <cassert>
#include <fstream>
#include <stdexcept>
#include <unordered_map>

//#define  VTB_DEBUG_YULE_WALKER_SOLVER
#if defined(VTB_DEBUG_YULE_WALKER_SOLVER)
#include <iostream>
#endif

namespace {

    inline blitz::TinyVector<int,1>
    vector_shape(int a) {
        return blitz::shape(3*a*a + 3*a + 1);
    }

    inline blitz::TinyVector<int,2>
    matrix_shape(int a, int b) {
        return blitz::shape(3*a*a + 3*a + 1, 3*b*b + 3*b + 1);
    }

    inline blitz::TinyVector<int,2>
    submatrix_shape(int a, int b) {
        return blitz::shape(6*a, 3*b*b + 3*b + 1);
    }

    template <class T>
    inline vtb::linalg::Vector<T>
    rho_vector(int i, int j, int k, int d, const blitz::Array<T,3>& acf) {
        #define ACF(x,y,z) acf(std::abs(x), std::abs(y), std::abs(z))
        typedef vtb::linalg::Vector<T> vector_type;
        using blitz::Range;
        vector_type rho_ijkd(vector_shape(d));
        T* data = rho_ijkd.data();
        *data++ = ACF(i-d, j-d, k-d);
        for (int e=1; e<=d; ++e) {
            for (int idx=0; idx<e; ++idx) {
                *data++ = ACF(i-d+e, j-d+idx, k-d);
            }
            for (int idx=0; idx<e; ++idx) {
                *data++ = ACF(i-d+e-idx, j-d+e, k-d);
            }
            for (int idx=0; idx<e; ++idx) {
                *data++ = ACF(i-d, j-d+e, k-d+idx);
            }
            for (int idx=0; idx<e; ++idx) {
                *data++ = ACF(i-d, j-d+e-idx, k-d+e);
            }
            for (int idx=0; idx<e; ++idx) {
                *data++ = ACF(i-d+idx, j-d, k-d+e);
            }
            for (int idx=0; idx<e; ++idx) {
                *data++ = ACF(i-d+e, j-d, k-d+e-idx);
            }
        }
        #undef ACF
        return rho_ijkd;
    }

    template <class T>
    inline T
    rho_scalar(int i, int j, int k, const blitz::Array<T,3>& acf) {
        #define ACF(x,y,z) acf(std::abs(x), std::abs(y), std::abs(z))
        return ACF(i, j, k);
        #undef ACF
    }

    template <class T>
    inline vtb::linalg::Matrix<T>
    R_matrix(int a, int b, const blitz::Array<T,3>& acf) {
        using blitz::Range;
        typedef vtb::linalg::Matrix<T> matrix_type;
        matrix_type R_ab(matrix_shape(a, b));
        R_ab(0, Range::all()) = rho_vector(a, a, a, b, acf);
        int offset = 1;
        for (int e=1; e<=a; ++e) {
            for (int idx=0; idx<e; ++idx) {
                R_ab(offset++, Range::all()) = rho_vector(a-e, a-idx, a, b, acf);
            }
            for (int idx=0; idx<e; ++idx) {
                R_ab(offset++, Range::all()) = rho_vector(a-e+idx, a-e, a, b, acf);
            }
            for (int idx=0; idx<e; ++idx) {
                R_ab(offset++, Range::all()) = rho_vector(a, a-e, a-idx, b, acf);
            }
            for (int idx=0; idx<e; ++idx) {
                R_ab(offset++, Range::all()) = rho_vector(a, a-e+idx, a-e, b, acf);
            }
            for (int idx=0; idx<e; ++idx) {
                R_ab(offset++, Range::all()) = rho_vector(a-idx, a, a-e, b, acf);
            }
            for (int idx=0; idx<e; ++idx) {
                R_ab(offset++, Range::all()) = rho_vector(a-e, a, a-e+idx, b, acf);
            }
        }
        return R_ab;
    }

    template <class T>
    inline vtb::linalg::Vector<T>
    R_vector_b0(int a, const blitz::Array<T,3>& acf) {
        using blitz::Range;
        using blitz::shape;
        typedef vtb::linalg::Vector<T> vector_type;
        vector_type R_ab(vector_shape(a));
        T* data = R_ab.data();
        *data++ = rho_scalar(a, a, a, acf);
        for (int e=1; e<=a; ++e) {
            for (int idx=0; idx<e; ++idx) {
                *data++ = rho_scalar(a-e, a-idx, a, acf);
            }
            for (int idx=0; idx<e; ++idx) {
                *data++ = rho_scalar(a-e+idx, a-e, a, acf);
            }
            for (int idx=0; idx<e; ++idx) {
                *data++ = rho_scalar(a, a-e, a-idx, acf);
            }
            for (int idx=0; idx<e; ++idx) {
                *data++ = rho_scalar(a, a-e+idx, a-e, acf);
            }
            for (int idx=0; idx<e; ++idx) {
                *data++ = rho_scalar(a-idx, a, a-e, acf);
            }
            for (int idx=0; idx<e; ++idx) {
                *data++ = rho_scalar(a-e, a, a-e+idx, acf);
            }
        }
        return R_ab;
    }

    template <class T>
    inline vtb::linalg::Vector<T>
    R_vector_a0(int b, blitz::Array<T,3> acf) {
        return rho_vector(0, 0, 0, b, acf);
    }

    template <class T>
    inline blitz::Array<T,3>
    result_array(blitz::Array<vtb::linalg::Vector<T>,1> beta, int order) {
        blitz::Array<T,3> result(blitz::shape(order, order, order));
        result = 0;
        result(0,0,0) = 0;
        for (int d=1; d<order; ++d) {
            const T* data = beta(d).data();
            result(d,d,d) = *data++;
            for (int e=1; e<=d; ++e) {
                for (int idx=0; idx<e; ++idx) {
                    result(d-e,d-idx,d) = *data++;
                }
                for (int idx=0; idx<e; ++idx) {
                    result(d-e+idx,d-e,d) = *data++;
                }
                for (int idx=0; idx<e; ++idx) {
                    result(d,d-e,d-idx) = *data++;
                }
                for (int idx=0; idx<e; ++idx) {
                    result(d,d-e+idx,d-e) = *data++;
                }
                for (int idx=0; idx<e; ++idx) {
                    result(d-idx,d,d-e) = *data++;
                }
                for (int idx=0; idx<e; ++idx) {
                    result(d-e,d,d-e+idx) = *data++;
                }
            }
        }
        return result;
    }

    /*
    template <class T>
    inline blitz::Array<T,2>
    operator*(blitz::Array<T,2> lhs, blitz::Array<T,2> rhs) {
        return linalg::multiply(lhs, rhs);
    }
    */

    template <class T, int N>
    inline bool
    is_square(const blitz::Array<T,N>& rhs) {
        const int first = rhs.extent(0);
        bool result = true;
        for (int i=1; i<rhs.dimensions(); ++i) {
            if (rhs.extent(i) != first) {
                result = false;
                break;
            }
        }
        return result;
    }

    template <class T>
    void
    append_column_block(
        vtb::linalg::Matrix<T>& lhs,
        const vtb::linalg::Matrix<T>& rhs
    ) {
        if (lhs.numElements() == 0) {
            lhs.resize(rhs.shape());
            lhs = rhs;
        } else {
            using blitz::Range;
            #if defined(VTB_DEBUG_YULE_WALKER_SOLVER)
            if (lhs.rows() != rhs.rows()) {
                throw std::runtime_error{"bad shape"};
            }
            #endif
            const int old_cols = lhs.columns();
            lhs.resizeAndPreserve(lhs.rows(), old_cols + rhs.columns());
            lhs(Range::all(), Range(old_cols, blitz::toEnd)) = rhs;
        }
    }

    template <class T>
    void
    append_row_block(
        vtb::linalg::Matrix<T>& lhs,
        const vtb::linalg::Matrix<T>& rhs
    ) {
        if (lhs.numElements() == 0) {
            lhs.resize(rhs.shape());
            lhs = rhs;
        } else {
            using blitz::Range;
            #if defined(VTB_DEBUG_YULE_WALKER_SOLVER)
            if (lhs.cols() != rhs.cols()) {
                throw std::runtime_error{"bad shape"};
            }
            #endif
            const int old_rows = lhs.rows();
            lhs.resizeAndPreserve(old_rows + rhs.rows(), lhs.columns());
            lhs(Range(old_rows, blitz::toEnd), Range::all()) = rhs;
        }
    }

}

namespace std {

    template <>
    struct hash<std::tuple<int,int>>: public std::hash<int> {
        inline size_t
        operator()(const std::tuple<int,int>& rhs) const noexcept {
            typedef std::hash<int> base;
            return base::operator()(std::get<0>(rhs)) ^
                   base::operator()(std::get<1>(rhs));
        }
    };

    template <>
    struct hash<std::tuple<int,int,int>>: public std::hash<int> {
        inline size_t
        operator()(const std::tuple<int,int,int>& rhs) const noexcept {
            typedef std::hash<int> base;
            return base::operator()(std::get<0>(rhs)) ^
                   base::operator()(std::get<1>(rhs)) ^
                   base::operator()(std::get<2>(rhs));
        }
    };

}

template <class T>
blitz::TinyVector<int,3>
vtb::core::chop_right(const blitz::Array<T,3>& rhs, T eps) {
    using blitz::abs;
    using blitz::all;
    using blitz::Range;
    const auto& _ = Range::all();
    // shrink third dimension
    int k = rhs.extent(2) - 1;
    while (k >= 1 && all(abs(rhs(_,_,k)) < eps)) {
        --k;
    }
    // shrink second dimension
    int j = rhs.extent(1) - 1;
    while (j >= 1 && all(abs(rhs(_, j, Range(0,k))) < eps)) {
        --j;
    }
    // shrink first dimension
    int i = rhs.extent(0) - 1;
    while (i >= 1 && all(abs(rhs(i, Range(0,j), Range(0,k))) < eps)) {
        --i;
    }
    return {i+1,j+1,k+1};
}

template <class T>
blitz::TinyVector<int,4>
vtb::core::chop_right(const blitz::Array<T,4>& rhs, T eps) {
    using blitz::abs;
    using blitz::all;
    using blitz::Range;
    const auto& _ = Range::all();
    // shrink fourth dimension
    int l = rhs.extent(3) - 1;
    while (l >= 1 && all(abs(rhs(_,_,_,l)) < eps)) {
        --l;
    }
    // shrink third dimension
    int k = rhs.extent(2) - 1;
    while (k >= 1 && all(abs(rhs(_, _, k, Range(0,l)) < eps))) {
        --k;
    }
    // shrink second dimension
    int j = rhs.extent(1) - 1;
    while (j >= 1 && all(abs(rhs(_, j, Range(0,k), Range(0,l))) < eps)) {
        --j;
    }
    // shrink first dimension
    int i = rhs.extent(0) - 1;
    while (i >= 1 && all(abs(rhs(i, Range(0,j), Range(0,k), Range(0,l))) < eps)) {
        --i;
    }
    return {i+1,j+1,k+1,l+1};
}

template <class T, int N>
auto
vtb::core::Choi_yule_walker_solver<T,N>
::solve(array_type acf_in) -> array_type {
    VTB_PROFILE_BLOCK("Choi_yule_walker_solver");
    #define PHI(i,j,k) Phi[std::make_tuple(i,j,k)]
    typedef vtb::linalg::Matrix<T> matrix_type;
    typedef vtb::linalg::Vector<T> vector_type;
    typedef std::tuple<int,int,int> key_type;
    typedef std::unordered_map<key_type,matrix_type> sparse_array_type;
    typedef blitz::Array<vector_type,1> vector_array_type;
    typedef blitz::Array<matrix_type,1> matrix_array_type;
    using blitz::Range;
    using blitz::shape;
    using blitz::any;
    if (!is_square(acf_in)) {
        throw std::invalid_argument("ACF is not square");
    }
    const T variance = acf_variance(acf_in);
    if (variance <= 0) {
        throw std::invalid_argument("ACF variance <= 0");
    }
    array_type acf(acf_in/variance);
    int max_order = this->max_order();
    if (max_order == 0) {
        max_order = blitz::max(acf.shape())-1;
    }
    sparse_array_type Phi;
    vector_array_type Pi_l(shape(max_order+2));
    vector_array_type Pi_lp1(shape(max_order+2));
    matrix_array_type Theta(shape(max_order+2));
    matrix_type R_1_1 = R_matrix(1, 1, acf);
    vector_type R_1_0 = R_vector_b0(1, acf);
    vector_type R_0_1 = R_vector_a0(1, acf);
//  std::clog << "R_1_1=" << R_1_1 << std::endl;
//  std::clog << "R_1_0=" << R_1_0 << std::endl;
//  std::clog << "R_0_1=" << R_0_1 << std::endl;
    matrix_type R_sup_1_1(R_1_1.shape());
    R_sup_1_1 = R_1_1;
    #if defined(VTB_DEBUG_YULE_WALKER_SOLVER)
    if (!R_sup_1_1.is_symmetric()) {
        throw std::runtime_error{"R matrix is not symmetric"};
    }
    #endif
    R_sup_1_1.inverse_self();
//  std::clog << "R_1_1^{-1}=" << R_sup_1_1 << std::endl;
//  std::clog << "R_1_1*R_1_1^{-1}=" << (R_1_1*R_sup_1_1) << std::endl;
    vector_type Pi_2_1 = R_sup_1_1 * R_1_0;
//  std::clog << "Pi_2_1=" << Pi_2_1 << std::endl;
    T lambda = T(1) - blitz::dot(R_0_1, Pi_2_1);
    T var0 = variance;
    T var = variance*lambda;
    if (var <= 0) {
        throw std::runtime_error("white noise variance <= 0");
    }
    #if defined(VTB_DEBUG_YULE_WALKER_SOLVER)
    std::clog << __func__ << ':'
        << "order=" << 1
        << ",var=" << var
        << ",lambda=" << lambda
        << '\n';
    #endif
    Pi_l(1).reference(Pi_2_1);
    PHI(2,1,0).reference(matrix_type(R_sup_1_1*R_matrix(1, 2, acf)));
    bool changed = false;
    int l = 1;
    do {
        ++l;
        var0 = var;
        matrix_type R_l_l = R_matrix(l, l, acf);
        vector_type R_l_0 = R_vector_b0(l, acf);
        matrix_type sum1(R_l_l.shape());
        sum1 = 0;
        vector_type sum2(R_l_0.shape());
        sum2 = 0;
        for (int m=1; m<=l-1; ++m) {
            matrix_type R_l_m(R_matrix(l, m, acf));
            sum1 += R_l_m*PHI(l,m,0);
//          sum2 += linalg::multiply_by_column_vector(R_l_m, Pi_l(m));
//          R_l_m.transpose_self();
            sum2 += R_l_m.transpose() * Pi_l(m);
        }
        matrix_type Theta_l(R_l_l - sum1);
        vector_type h_l(R_l_0 - sum2);
//      std::clog << "h_l=" << h_l << std::endl;
//      std::clog << "Theta_l=" << Theta_l << std::endl;
        Theta_l.inverse_self();
        Theta(l).reference(Theta_l);
//      { std::ifstream("/tmp/z") >> Theta_l; }
//      std::clog << "Theta_l^{-1}=" << Theta_l << std::endl;
        Pi_lp1(l).reference(Theta_l * h_l);
//      std::clog << "Pi(l+1,l)=" << Pi_lp1(l) << std::endl;
        for (int a=1; a<=l-1; ++a) {
            Pi_lp1(a).reference(vector_type(
//              Pi_l(a) - linalg::multiply_by_column_vector(PHI(l,a,0), Pi_lp1(l))
                Pi_l(a) - PHI(l,a,0).transpose() * Pi_lp1(l)
            ));
        }
        lambda -= blitz::dot(h_l, Pi_lp1(l));
        var = variance*lambda;
        if (var <= 0) {
            throw std::runtime_error("white noise variance <= 0");
        }
        if (l < max_order) {
            PHI(2,1,l-1).reference(R_sup_1_1*R_matrix(1, l+1, acf));
            for (int m=2; m<=l; ++m) {
                matrix_type R_m_lp1(R_matrix(m,l+1,acf));
                matrix_type sum3(R_m_lp1.shape());
                sum3 = 0;
                for (int n=1; n<=m-1; ++n) {
                    sum3 += R_matrix(m,n,acf)*PHI(m,n,l-m+1);
                }
                PHI(m+1,m,l-m).reference(Theta(m)*matrix_type(R_m_lp1 - sum3));
                for (int a=1; a<=m-1; ++a) {
                    PHI(m+1,a,l-m).reference(matrix_type(
                        PHI(m,a,l-m+1) - PHI(m,a,0)*PHI(m+1,m,l-m)
                    ));
                }
            }
        }
        blitz::swap(Pi_l, Pi_lp1);
        #if defined(VTB_DEBUG_YULE_WALKER_SOLVER)
        std::clog
            << __func__ << ':'
            << "order=" << l
            << ",var=" << var
            << ",lambda=" << lambda
            << '\n';
        #endif
        changed = !this->variance_has_not_changed_much(var, var0);
    } while (l < max_order && changed);
    #undef PHI
    array_type result;
    this->white_noise_variance(var);
    result.reference(result_array(Pi_lp1, l));
    if (this->_chop) {
        result.resizeAndPreserve(chop_right(result, this->_chopepsilon));
    }
    return result;
}

template blitz::TinyVector<int,3>
vtb::core::chop_right(
    const blitz::Array<VTB_REAL_TYPE,3>& rhs,
    VTB_REAL_TYPE eps
);


template <class T>
auto
vtb::core::Autocovariance_matrix_generator<T,3>::
AC_matrix_block(int i0, int j0) -> matrix_type {
    const int n = _arorder(2);
    matrix_type block(blitz::shape(n, n));
    for (int k=0; k<n; ++k) {
        for (int j=0; j<n; ++j) {
            block(k,j) = this->_acf(i0, j0, std::abs(k - j));
        }
    }
    return block;
}

template <class T>
auto
vtb::core::Autocovariance_matrix_generator<T,3>::
AC_matrix_block(int i0) -> matrix_type {
    const int n = _arorder(1);
    matrix_type result;
    for (int i = 0; i < n; ++i) {
        matrix_type row;
        for (int j = 0; j < n; ++j) {
            matrix_type tmp = AC_matrix_block(i0, std::abs(i - j));
            append_column_block(row, tmp);
        }
        append_row_block(result, row);
    }
    return result;
}

template <class T>
auto
vtb::core::Autocovariance_matrix_generator<T,3>::
operator()() -> matrix_type {
    const int n = _arorder(0);
    matrix_type result;
    for (int i = 0; i < n; ++i) {
        matrix_type row;
        for (int j = 0; j < n; ++j) {
            matrix_type tmp = AC_matrix_block(std::abs(i - j));
            append_column_block(row, tmp);
        }
        append_row_block(result, row);
    }
    return result;
}

template <class T, int N>
auto
vtb::core::Gauss_yule_walker_solver<T,N>
::solve(array_type acf) -> array_type {
    VTB_PROFILE_BLOCK("Gauss_yule_walker_solver");
    using blitz::Range;
    using blitz::toEnd;
    using vtb::linalg::Positive_definite_matrix;
    using vtb::linalg::Square_matrix;
    using vtb::linalg::Symmetric_matrix;
    using vtb::linalg::Vector;
    using vtb::linalg::solve;
    blitz::TinyVector<int,N> order;
    order = where(this->order() == 0, acf.shape()-1, this->order());
    #if defined(VTB_DEBUG_YULE_WALKER_SOLVER)
    std::clog << "order=" << order << std::endl;
    #endif
    Autocovariance_matrix_generator<T,N> gen(acf, order);
    Square_matrix<T> acm(gen());
    const int m = acm.rows() - 1;
    // eliminate the first equation and move the first column of the
    // remaining matrix to the right-hand side of the system
    Vector<T> rhs(m);
    rhs = acm(Range(1, toEnd), 0);
    // lhs is the autocovariance matrix without first
    // column and row
    //Positive_definite_matrix<T> lhs(blitz::shape(m, m));
    Symmetric_matrix<T> lhs(blitz::shape(m, m));
    lhs = acm(Range(1, toEnd), Range(1, toEnd));
    #if defined(VTB_DEBUG_YULE_WALKER_SOLVER)
    if (lhs.extent(0) != m ||
        lhs.extent(1) != m ||
        rhs.extent(0) != m ||
        !lhs.is_symmetric()) {
        throw std::runtime_error{"bad autocovariance matrix"};
    }
    dbg::write_matrix_mathematica("lhs", lhs);
    dbg::write_vector_mathematica("rhs", rhs);
    #endif
    Vector<T> x = solve(lhs, rhs);
    #if defined(VTB_DEBUG_YULE_WALKER_SOLVER)
    const T residual = max(abs(rhs - lhs*x));
    std::clog << "residual=" << residual << std::endl;
    #endif
    array_type result(order);
    result.data()[0] = 0;
    std::copy_n(x.data(), x.size(), result.data()+1);
    this->white_noise_variance(
        autoregression_white_noise_variance(acf, result)
    );
    if (this->_chop) {
        result.resizeAndPreserve(chop_right(result, this->_chopepsilon));
    }
    return result;
}

template <class T>
auto
vtb::core::Autocovariance_matrix_generator<T,4>::
AC_matrix_block(int i0, int j0, int k0) -> matrix_type {
    const int n = _arorder(3);
    matrix_type block(blitz::shape(n, n));
    for (int l=0; l<n; ++l) {
        for (int k=0; k<n; ++k) {
            block(l,k) = this->_acf(i0, j0, k0, std::abs(l - k));
        }
    }
    return block;
}

template <class T>
auto
vtb::core::Autocovariance_matrix_generator<T,4>::
AC_matrix_block(int i0, int j0) -> matrix_type {
    const int n = _arorder(2);
    matrix_type result;
    for (int k = 0; k < n; ++k) {
        matrix_type row;
        for (int j = 0; j < n; ++j) {
            matrix_type tmp = AC_matrix_block(i0, j0, std::abs(k - j));
            append_column_block(row, tmp);
        }
        append_row_block(result, row);
    }
    return result;
}

template <class T>
auto
vtb::core::Autocovariance_matrix_generator<T,4>::
AC_matrix_block(int i0) -> matrix_type {
    const int n = _arorder(1);
    matrix_type result;
    for (int i = 0; i < n; ++i) {
        matrix_type row;
        for (int j = 0; j < n; ++j) {
            matrix_type tmp = AC_matrix_block(i0, std::abs(i - j));
            append_column_block(row, tmp);
        }
        append_row_block(result, row);
    }
    return result;
}

template <class T>
auto
vtb::core::Autocovariance_matrix_generator<T,4>::
operator()() -> matrix_type {
    const int n = _arorder(0);
    matrix_type result;
    for (int i = 0; i < n; ++i) {
        matrix_type row;
        for (int j = 0; j < n; ++j) {
            matrix_type tmp = AC_matrix_block(std::abs(i - j));
            append_column_block(row, tmp);
        }
        append_row_block(result, row);
    }
    return result;
}

template <class T, int N>
auto
vtb::core::Fixed_point_iteration_yule_walker_solver<T,N>
::solve(array_type acf) -> array_type {
    VTB_PROFILE_BLOCK("Fixed_point_iteration_yule_walker_solver");
    typedef blitz::TinyVector<int,N> index;
    using blitz::RectDomain;
    using blitz::abs;
    using blitz::all;
    using blitz::isfinite;
    using blitz::max;
    using blitz::sum;
    using std::abs;
    if (!all(this->_order <= acf.shape())) {
        throw std::invalid_argument("bad order");
    }
    const int max_iterations = this->_maxiterations;
    const T min_var_wn = this->_minvarwn;
    const T eps = this->_epsvarwn;
    const T max_residual = this->_maxresidual;
    const shape_type& order = this->_order;
    array_type theta(order);
    theta = 0;
    const int ni = order(0);
    const int nj = order(1);
    const int nk = order(2);
    T var_wn = acf(0,0,0);
    T old_var_wn = 0;
    T residual = 0;
    int it = 0;
    do {
        theta(0,0,0) = 0;
        for (int i = ni - 1; i >= 0; --i) {
            for (int j = nj - 1; j >= 0; --j) {
                for (int k = nk - 1; k >= 0; --k) {
                    const shape_type ijk{i,j,k};
                    RectDomain<3> sub1(ijk, order - 1);
                    RectDomain<3> sub2({0}, order - ijk - 1);
                    theta(i,j,k) = -acf(i,j,k) / var_wn +
                        sum(theta(sub1)*theta(sub2));
                }
            }
        }
        theta(0,0,0) = -1;
        if (!all(isfinite(theta))) {
            throw std::runtime_error("bad MA model coefficients");
        }
        theta(0,0,0) = -1;
        residual = std::numeric_limits<T>::min();
        for_loop<3>(
            order,
            [&] (const index& idx) {
                RectDomain<3> sub1(idx, order-1);
                RectDomain<3> sub2({0}, order-idx-1);
                T new_residual =
                    abs(acf(idx) - sum(theta(sub1)*theta(sub2))*var_wn);
                if (residual < new_residual) {
                    residual = new_residual;
                }
            }
        );
        old_var_wn = var_wn;
        var_wn = moving_average_white_noise_variance(acf, theta);
        if (var_wn <= min_var_wn) {
            std::stringstream msg;
            msg << "bad white noise variance: " << var_wn;
            throw std::runtime_error(msg.str());
        }
        #if defined(VTB_DEBUG_YULE_WALKER_SOLVER)
        std::clog << __func__ << ':' << "Iteration=" << it
                  << ",var_wn=" << var_wn
                  << ",resudual=" << residual
                  << ",theta(0,0,0)=" << theta(0,0,0)
                  << std::endl;
        #endif
        ++it;
    } while ((it < max_iterations) &&
             abs(var_wn - old_var_wn) > eps &&
             (residual > max_residual));
    #if defined(VTB_DEBUG_YULE_WALKER_SOLVER)
    std::clog << "Calculated MA model coefficients:"
              << "\n\tno. of iterations = " << it
              << "\n\twhite noise variance = " << var_wn
              << "\n\tmax(theta) = " << max(abs(theta))
              << "\n\tmax(residual) = " << residual
              << std::endl;
    #endif
    this->white_noise_variance(var_wn);
    return theta;
}

template class vtb::core::Choi_yule_walker_solver<VTB_REAL_TYPE,3>;
template class vtb::core::Gauss_yule_walker_solver<VTB_REAL_TYPE,3>;
template class vtb::core::Gauss_yule_walker_solver<VTB_REAL_TYPE,4>;
template class vtb::core::Autocovariance_matrix_generator<VTB_REAL_TYPE,3>;
template class vtb::core::Autocovariance_matrix_generator<VTB_REAL_TYPE,4>;
template class vtb::core::Fixed_point_iteration_yule_walker_solver<VTB_REAL_TYPE,3>;
#if defined(VTB_REAL_TYPE_FLOAT)
template class vtb::core::Choi_yule_walker_solver<double,3>;
template class vtb::core::Gauss_yule_walker_solver<double,3>;
template class vtb::core::Gauss_yule_walker_solver<double,4>;
template class vtb::core::Autocovariance_matrix_generator<double,3>;
template class vtb::core::Autocovariance_matrix_generator<double,4>;
#endif