Toeplitz.h

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
// SPDX-FileCopyrightText: The Eigen Authors
// SPDX-License-Identifier: MPL-2.0
 
// References:
//  [1] R. M. Gray, "Toeplitz and Circulant Matrices: A Review", Foundations and
//      Trends in Communications and Information Theory, 2(3), 2006.
//  [2] G. H. Golub and C. F. Van Loan, "Matrix Computations", 4th ed., Johns
//      Hopkins University Press, 2013, chapter 4.8 (fast Toeplitz products via
//      circulant embedding and the FFT).
 
#ifndef EIGEN_STRUCTURED_TOEPLITZ_H
#define EIGEN_STRUCTURED_TOEPLITZ_H
 
// IWYU pragma: private
#include "./InternalHeaderCheck.h"
 
namespace Eigen {
 
template <typename Scalar_, int Rows_ = Dynamic, int Cols_ = Dynamic>
class Toeplitz;
 
namespace internal {
 
template <typename Scalar_, int Rows_, int Cols_>
struct traits<Toeplitz<Scalar_, Rows_, Cols_>> {
  using Scalar = Scalar_;
  using StorageKind = Dense;
  using XprKind = MatrixXpr;
  using StorageIndex = int;
  static constexpr int RowsAtCompileTime = Rows_;
  static constexpr int ColsAtCompileTime = Cols_;
  static constexpr int MaxRowsAtCompileTime = Rows_;
  static constexpr int MaxColsAtCompileTime = Cols_;
  // Deliberately no NestByRefBit: transpose(), conjugate() and adjoint() return
  // owning temporaries, so Product must nest the operator by value for a
  // delayed-evaluated product expression to keep its left factor alive. The copy
  // is O(m+n), negligible against the O(n log n) product evaluation.
  static constexpr int Flags = 0;
};
 
template <typename Scalar_, int Rows_, int Cols_>
struct evaluator_traits<Toeplitz<Scalar_, Rows_, Cols_>> {
  using Kind = IndexBased;
  using Shape = StructuredShape;
};
 
}  // namespace internal

template <typename Scalar_, int Rows_, int Cols_>
class Toeplitz : public EigenBase<Toeplitz<Scalar_, Rows_, Cols_>> {
 public:
  using Scalar = Scalar_;
  using RealScalar = typename NumTraits<Scalar>::Real;
  using StorageIndex = int;
  using Complex = std::complex<RealScalar>;
  using ColGeneratorType = Matrix<Scalar, Rows_, 1>;
  using RowGeneratorType = Matrix<Scalar, Cols_, 1>;
  using ComplexVector = Matrix<Complex, Dynamic, 1>;
 
  static constexpr int RowsAtCompileTime = Rows_;
  static constexpr int ColsAtCompileTime = Cols_;
  static constexpr int MaxRowsAtCompileTime = Rows_;
  static constexpr int MaxColsAtCompileTime = Cols_;
  static constexpr int SizeAtCompileTime = internal::size_at_compile_time(Rows_, Cols_);
  static constexpr int MaxSizeAtCompileTime = SizeAtCompileTime;
  static constexpr bool IsRowMajor = false;
  // Deliberately no IsVectorAtCompileTime: Ref<const Toeplitz>'s default StrideType
  // argument reads it, so its absence makes internal::is_ref_compatible SFINAE to
  // false and keeps the iterative solvers on their matrix-free path.

  template <typename ColDerived, typename RowDerived>
  Toeplitz(const MatrixBase<ColDerived>& col, const MatrixBase<RowDerived>& row) : m_col(col), m_row(row) {
    EIGEN_STATIC_ASSERT_VECTOR_ONLY(ColDerived)
    EIGEN_STATIC_ASSERT_VECTOR_ONLY(RowDerived)
    eigen_assert(m_col.size() > 0 && m_row.size() > 0 && "Toeplitz generators must be non-empty");
    // Below the threshold every product takes the direct path and the symbol has
    // no consumer, so skip the embedding altogether.
    if (rows() > internal::structured_direct_threshold() || cols() > internal::structured_direct_threshold())
      m_symbol = computeSymbol();
  }
 
  EIGEN_DEVICE_FUNC Index rows() const { return m_col.size(); }
  EIGEN_DEVICE_FUNC Index cols() const { return m_row.size(); }

  const ColGeneratorType& column() const { return m_col; }
  const RowGeneratorType& row() const { return m_row; }

  ComplexVector symbol() const { return m_symbol.size() > 0 ? m_symbol : computeSymbol(); }

  Scalar coeff(Index row, Index col) const {
    Index k = row - col;
    return k >= 0 ? m_col.coeff(k) : m_row.coeff(-k);
  }

  Toeplitz<Scalar, Cols_, Rows_> transpose() const {
    using TransposeType = Toeplitz<Scalar, Cols_, Rows_>;
    const Index n = cols();
    typename TransposeType::ColGeneratorType col(n);
    col[0] = m_col[0];
    if (n > 1) col.tail(n - 1) = m_row.tail(n - 1);
    return TransposeType(col, m_col, internal::structured_reverse_symbol(m_symbol));
  }

  Toeplitz conjugate() const {
    return Toeplitz(m_col.conjugate(), m_row.conjugate(), internal::structured_reverse_symbol(m_symbol).conjugate());
  }

  Toeplitz<Scalar, Cols_, Rows_> adjoint() const {
    using AdjointType = Toeplitz<Scalar, Cols_, Rows_>;
    const Index n = cols();
    typename AdjointType::ColGeneratorType col(n);
    col[0] = numext::conj(m_col[0]);
    if (n > 1) col.tail(n - 1) = m_row.tail(n - 1).conjugate();
    return AdjointType(col, m_col.conjugate(), m_symbol.conjugate());
  }

  template <typename Dest>
  void evalTo(Dest& dst) const {
    const Index m = rows(), n = cols();
    EIGEN_IF_CONSTEXPR (Dest::IsRowMajor) {
      for (Index i = 0; i < m; ++i) {
        const Index h = numext::mini(i + 1, n);
        dst.row(i).head(h) = m_col.segment(i - h + 1, h).reverse().transpose();
        if (i + 1 < n) dst.row(i).tail(n - i - 1) = m_row.segment(1, n - i - 1).transpose();
      }
      return;
    }
    for (Index j = 0; j < n; ++j) {
      const Index h = numext::mini(j, m);
      dst.col(j).head(h) = m_row.segment(j - h + 1, h).reverse();
      if (j < m) dst.col(j).tail(m - j) = m_col.head(m - j);
    }
  }

  template <typename Dest>
  void addTo(Dest& dst) const {
    const Index m = rows(), n = cols();
    EIGEN_IF_CONSTEXPR (Dest::IsRowMajor) {
      for (Index i = 0; i < m; ++i) {
        const Index h = numext::mini(i + 1, n);
        dst.row(i).head(h) += m_col.segment(i - h + 1, h).reverse().transpose();
        if (i + 1 < n) dst.row(i).tail(n - i - 1) += m_row.segment(1, n - i - 1).transpose();
      }
      return;
    }
    for (Index j = 0; j < n; ++j) {
      const Index h = numext::mini(j, m);
      dst.col(j).head(h) += m_row.segment(j - h + 1, h).reverse();
      if (j < m) dst.col(j).tail(m - j) += m_col.head(m - j);
    }
  }

  template <typename Dest>
  void subTo(Dest& dst) const {
    const Index m = rows(), n = cols();
    EIGEN_IF_CONSTEXPR (Dest::IsRowMajor) {
      for (Index i = 0; i < m; ++i) {
        const Index h = numext::mini(i + 1, n);
        dst.row(i).head(h) -= m_col.segment(i - h + 1, h).reverse().transpose();
        if (i + 1 < n) dst.row(i).tail(n - i - 1) -= m_row.segment(1, n - i - 1).transpose();
      }
      return;
    }
    for (Index j = 0; j < n; ++j) {
      const Index h = numext::mini(j, m);
      dst.col(j).head(h) -= m_row.segment(j - h + 1, h).reverse();
      if (j < m) dst.col(j).tail(m - j) -= m_col.head(m - j);
    }
  }

  template <typename Rhs>
  Product<Toeplitz, Rhs> operator*(const MatrixBase<Rhs>& x) const {
    EIGEN_STATIC_ASSERT(ColsAtCompileTime == Dynamic || Rhs::RowsAtCompileTime == Dynamic ||
                            int(ColsAtCompileTime) == int(Rhs::RowsAtCompileTime),
                        INVALID_MATRIX_PRODUCT)
    eigen_assert(x.rows() == cols() && "invalid product: dimensions do not match");
    return Product<Toeplitz, Rhs>(*this, x.derived());
  }

  template <typename Dest, typename Rhs, typename ProductScalar>
  void addProduct(Dest& dst, const Rhs& rhs, const ProductScalar& alpha) const {
    const Index m = rows(), n = cols();
    eigen_assert(rhs.rows() == n && "invalid product: dimensions do not match");
    const bool small = m <= internal::structured_direct_threshold() && n <= internal::structured_direct_threshold();
    if (small)
      directProduct(dst, rhs, alpha);
    else
      internal::structured_fft_apply(dst, m_symbol, m, rhs, alpha);
  }
 
 private:
  // Grants transpose() and adjoint() access to the private constructor of the
  // dimension-swapped instantiation.
  template <typename OtherScalar, int OtherRows, int OtherCols>
  friend class Toeplitz;

  Toeplitz(const ColGeneratorType& col, const RowGeneratorType& row, const ComplexVector& symbol)
      : m_col(col), m_row(row), m_symbol(symbol) {}

  template <typename Dest, typename Rhs, typename ProductScalar>
  void directProductColumn(Dest& dst, const Rhs& rhs, Index k, const ProductScalar& alpha) const {
    const Index m = rows(), n = cols();
    // A unit alpha must not multiply: even the identity complex scalar (1,0)
    // pollutes an (Inf,0) value with NaN through the 0*Inf cross term.
    const bool unitAlpha = alpha == ProductScalar(1);
    if (m <= internal::structured_scalar_threshold() && n <= internal::structured_scalar_threshold()) {
      // Tiny sizes: a plain scalar loop beats the segment-based path below, whose
      // per-segment setup dominates when segments hold only a few entries.
      for (Index i = 0; i < m; ++i) {
        ProductScalar acc(0);
        for (Index j = 0; j < n; ++j) acc += coeff(i, j) * rhs.coeff(j, k);
        dst.coeffRef(i, k) += unitAlpha ? acc : ProductScalar(alpha * acc);
      }
      return;
    }
 
    // Split each column at its diagonal into contiguous generator slices.
    auto dstCol = dst.col(k);
    for (Index j = 0; j < n; ++j) {
      const ProductScalar xj = unitAlpha ? ProductScalar(rhs.coeff(j, k)) : ProductScalar(alpha * rhs.coeff(j, k));
      const Index h = numext::mini(j, m);
      dstCol.head(h) += xj * m_row.segment(j - h + 1, h).reverse();
      if (j < m) dstCol.tail(m - j) += xj * m_col.head(m - j);
    }
  }

  template <typename Dest, typename Rhs, typename ProductScalar>
  void directProduct(Dest& dst, const Rhs& rhs, const ProductScalar& alpha) const {
    for (Index k = 0; k < rhs.cols(); ++k) directProductColumn(dst, rhs, k, alpha);
  }

  ComplexVector computeSymbol() const {
    const Index m = rows(), n = cols();
    const Index p = internal::fft_next_good_size(m + n - 1);
    ComplexVector embedding = ComplexVector::Zero(p);
    embedding.head(m) = m_col.template cast<Complex>();
    embedding.tail(n - 1) = m_row.tail(n - 1).reverse().template cast<Complex>();
    if (p == 1) return embedding;  // the DFT of a single sample is the identity
    ComplexVector symbol(p);
    auto&& fft = internal::structured_fft_engine<RealScalar>();
    fft.fwd(symbol, embedding, p);
    return symbol;
  }
 
  ColGeneratorType m_col;
  RowGeneratorType m_row;
  ComplexVector m_symbol;
};

template <typename ColDerived, typename RowDerived>
Toeplitz<typename ColDerived::Scalar, ColDerived::SizeAtCompileTime, RowDerived::SizeAtCompileTime> makeToeplitz(
    const MatrixBase<ColDerived>& col, const MatrixBase<RowDerived>& row) {
  return Toeplitz<typename ColDerived::Scalar, ColDerived::SizeAtCompileTime, RowDerived::SizeAtCompileTime>(col, row);
}
 
namespace internal {
 
template <typename Scalar_, int Rows_, int Cols_, typename Rhs, int ProductTag>
struct generic_product_impl<Toeplitz<Scalar_, Rows_, Cols_>, Rhs, StructuredShape, DenseShape, ProductTag>
    : structured_product_impl<Toeplitz<Scalar_, Rows_, Cols_>, Rhs> {};
 
}  // namespace internal
 
}  // namespace Eigen
 
#endif  // EIGEN_STRUCTURED_TOEPLITZ_H