Hankel.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] G. H. Golub and C. F. Van Loan, "Matrix Computations", 4th ed., Johns
//      Hopkins University Press, 2013, chapter 4.7 (Toeplitz, Hankel and
//      related systems) and chapter 4.8 (fast products via circulant embedding
//      and the FFT).
//  [2] R. H. Chan and X.-Q. Jin, "An Introduction to Iterative Toeplitz
//      Solvers", SIAM, 2007 (a Hankel matrix is a column-reversed Toeplitz
//      matrix, so its products and solves reduce to the Toeplitz machinery).
 
#ifndef EIGEN_STRUCTURED_HANKEL_H
#define EIGEN_STRUCTURED_HANKEL_H
 
// IWYU pragma: private
#include "./InternalHeaderCheck.h"
 
namespace Eigen {
 
template <typename Scalar_, int Rows_ = Dynamic, int Cols_ = Dynamic>
class Hankel;
 
namespace internal {
 
template <typename Scalar_, int Rows_, int Cols_>
struct traits<Hankel<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(n), negligible against the O(n log n) product evaluation.
  static constexpr unsigned int Flags = 0;
};
 
template <typename Scalar_, int Rows_, int Cols_>
struct evaluator_traits<Hankel<Scalar_, Rows_, Cols_>> {
  using Kind = IndexBased;
  using Shape = StructuredShape;
};
 
}  // namespace internal

template <typename Scalar_, int Rows_, int Cols_>
class Hankel : public EigenBase<Hankel<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 Hankel>'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.

  static constexpr int GenSizeAtCompileTime = (Rows_ == Dynamic || Cols_ == Dynamic) ? Dynamic : (Rows_ + Cols_ - 1);
  using GeneratorType = Matrix<Scalar, GenSizeAtCompileTime, 1>;

  template <typename ColDerived, typename RowDerived>
  Hankel(const MatrixBase<ColDerived>& col, const MatrixBase<RowDerived>& row)
      : m_rows(col.size()), m_cols(row.size()) {
    EIGEN_STATIC_ASSERT_VECTOR_ONLY(ColDerived)
    EIGEN_STATIC_ASSERT_VECTOR_ONLY(RowDerived)
    const Index m = col.size(), n = row.size();
    eigen_assert(m > 0 && n > 0 && "Hankel generators must be non-empty");
    m_h.resize(m + n - 1);
    m_h.head(m) = col;
    if (n > 1) m_h.tail(n - 1) = row.tail(n - 1);
    if (m > 1 && n > 1 && (m > internal::structured_direct_threshold() || n > internal::structured_direct_threshold()))
      m_symbol = computeSymbol();
  }
 
  EIGEN_DEVICE_FUNC Index rows() const { return m_rows.value(); }
  EIGEN_DEVICE_FUNC Index cols() const { return m_cols.value(); }

  const GeneratorType& generator() const { return m_h; }

  ColGeneratorType column() const { return m_h.head(rows()); }

  RowGeneratorType lastRow() const { return m_h.segment(rows() - 1, cols()); }

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

  Scalar coeff(Index row, Index col) const { return m_h.coeff(row + col); }

  Hankel<Scalar, Cols_, Rows_> transpose() const {
    return Hankel<Scalar, Cols_, Rows_>(m_h, cols(), rows(), transposedSymbol());
  }

  Hankel conjugate() const {
    return Hankel(m_h.conjugate(), rows(), cols(), internal::structured_reverse_symbol(m_symbol).conjugate());
  }

  Hankel<Scalar, Cols_, Rows_> adjoint() const { return conjugate().transpose(); }

  Toeplitz<Scalar, Rows_, Cols_> toToeplitz() const {
    ColGeneratorType c = m_h.segment(cols() - 1, rows());
    RowGeneratorType r = m_h.head(cols()).reverse();
    return Toeplitz<Scalar, Rows_, Cols_>(c, r);
  }

  static constexpr int SolveRowsAtCompileTime = Cols_ != Dynamic ? Cols_ : Rows_;

  template <typename Rhs>
  Matrix<Scalar, SolveRowsAtCompileTime, Rhs::ColsAtCompileTime> solve(const MatrixBase<Rhs>& b) const {
    EIGEN_STATIC_ASSERT(RowsAtCompileTime == Dynamic || Rhs::RowsAtCompileTime == Dynamic ||
                            int(RowsAtCompileTime) == int(Rhs::RowsAtCompileTime),
                        YOU_MIXED_MATRICES_OF_DIFFERENT_SIZES)
    eigen_assert(rows() == cols() && "Hankel::solve requires a square matrix");
    eigen_assert(b.rows() == rows() && "right-hand side has the wrong number of rows");
    LookAheadLevinson<Scalar> levinson(toToeplitz());
    eigen_assert(levinson.info() == Success &&
                 "Hankel::solve: the Levinson recursion broke down; "
                 "use LookAheadLevinson on toToeplitz() for diagnostics");
    return levinson.solve(b).colwise().reverse();
  }

  template <typename Dest>
  void evalTo(Dest& dst) const {
    EIGEN_IF_CONSTEXPR (Dest::IsRowMajor) {
      for (Index i = 0; i < rows(); ++i) dst.row(i) = m_h.segment(i, cols()).transpose();
      return;
    }
    for (Index j = 0; j < cols(); ++j) dst.col(j) = m_h.segment(j, rows());
  }

  template <typename Dest>
  void addTo(Dest& dst) const {
    EIGEN_IF_CONSTEXPR (Dest::IsRowMajor) {
      for (Index i = 0; i < rows(); ++i) dst.row(i) += m_h.segment(i, cols()).transpose();
      return;
    }
    for (Index j = 0; j < cols(); ++j) dst.col(j) += m_h.segment(j, rows());
  }

  template <typename Dest>
  void subTo(Dest& dst) const {
    EIGEN_IF_CONSTEXPR (Dest::IsRowMajor) {
      for (Index i = 0; i < rows(); ++i) dst.row(i) -= m_h.segment(i, cols()).transpose();
      return;
    }
    for (Index j = 0; j < cols(); ++j) dst.col(j) -= m_h.segment(j, rows());
  }

  template <typename Rhs>
  Product<Hankel, 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<Hankel, 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 == 1 || n == 1 ||
                       (m <= internal::structured_direct_threshold() && n <= internal::structured_direct_threshold());
    if (small) {
      directProduct(dst, rhs, alpha);
      return;
    }
    internal::structured_fft_apply(dst, m_symbol, m, rhs.colwise().reverse(), alpha);
  }
 
 private:
  template <typename OtherScalar, int OtherRows, int OtherCols>
  friend class Hankel;

  Hankel(const GeneratorType& h, Index rows, Index cols, const ComplexVector& symbol)
      : m_rows(rows), m_cols(cols), m_h(h), 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 (n == 1) {
      // The runtime-sized head keeps other fixed-size shapes well-formed.
      const ProductScalar xj = unitAlpha ? ProductScalar(rhs.coeff(0, k)) : ProductScalar(alpha * rhs.coeff(0, k));
      dst.col(k) += xj * m_h.head(m);
      return;
    }
 
    if (m == 1) {
      const ProductScalar acc = m_h.head(n).cwiseProduct(rhs.col(k)).sum();
      dst.coeffRef(0, k) += unitAlpha ? acc : ProductScalar(alpha * acc);
      return;
    }
 
    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;
    }
 
    // Column j is the contiguous, vectorizable slice h[j..j+m-1].
    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));
      dstCol += xj * m_h.segment(j, m);
    }
  }

  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_h.tail(m).template cast<Complex>();  // h[n-1 .. m+n-2]
    embedding.tail(n - 1) = m_h.head(n - 1).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;
  }

  ComplexVector transposedSymbol() const {
    ComplexVector sym = m_symbol;
    const Index p = sym.size();
    if (p > 0 && rows() != cols()) {
      Index s = (rows() - cols()) % p;
      if (s < 0) s += p;
      // Contiguous real buffers let sin/cos use packets; complex component views
      // alone have stride two and would force scalar trigonometric evaluation.
      Array<RealScalar, 64, 1> angles, values;
      Matrix<Complex, 64, 1> phases;
      Index fs = 0;  // f * s mod p
      for (Index f = 0; f < p;) {
        const Index count = numext::mini<Index>(64, p - f);
        for (Index i = 0; i < count; ++i) {
          angles[i] = RealScalar(2 * EIGEN_PI) * RealScalar(fs) / RealScalar(p);
          fs += s;
          if (fs >= p) fs -= p;
        }
        values.head(count) = angles.head(count).cos();
        phases.head(count).real() = values.head(count);
        values.head(count) = angles.head(count).sin();
        phases.head(count).imag() = values.head(count);
        sym.segment(f, count).array() *= phases.head(count).array();
        f += count;
      }
    }
    return sym;
  }
 
  // Empty (compile-time constant) when the corresponding dimension is fixed; the
  // dimensions cannot be recovered from m_h alone, whose length is rows()+cols()-1.
  internal::variable_if_dynamic<Index, RowsAtCompileTime> m_rows;
  internal::variable_if_dynamic<Index, ColsAtCompileTime> m_cols;
  GeneratorType m_h;
  ComplexVector m_symbol;
};

template <typename ColDerived, typename RowDerived>
Hankel<typename ColDerived::Scalar, ColDerived::SizeAtCompileTime, RowDerived::SizeAtCompileTime> makeHankel(
    const MatrixBase<ColDerived>& col, const MatrixBase<RowDerived>& row) {
  return Hankel<typename ColDerived::Scalar, ColDerived::SizeAtCompileTime, RowDerived::SizeAtCompileTime>(col, row);
}
 
namespace internal {
 
// Single product specialization covering every product tag; see the note in
// Circulant.h.
template <typename Scalar_, int Rows_, int Cols_, typename Rhs, int ProductTag>
struct generic_product_impl<Hankel<Scalar_, Rows_, Cols_>, Rhs, StructuredShape, DenseShape, ProductTag>
    : structured_product_impl<Hankel<Scalar_, Rows_, Cols_>, Rhs> {};
 
}  // namespace internal
 
}  // namespace Eigen
 
#endif  // EIGEN_STRUCTURED_HANKEL_H