Circulant.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] P. J. Davis, "Circulant Matrices", Wiley, 1979. Diagonalization of a
//      circulant matrix by the DFT and the closed-form eigenstructure used by
//      eigenvalues()/eigenvectors(); the SVD and pseudo-inverse below follow
//      from it by taking moduli/phases of the eigenvalues.
//  [2] R. M. Gray, "Toeplitz and Circulant Matrices: A Review", Foundations and
//      Trends in Communications and Information Theory, 2(3), 2006.
//  [3] G. H. Golub and C. F. Van Loan, "Matrix Computations", 4th ed., Johns
//      Hopkins University Press, 2013, chapter 4.8 (circulant systems and
//      FFT-based products) and chapter 5.4 (numerical rank conventions).
//  [4] J. J. Dongarra, J. R. Bunch, C. B. Moler and G. W. Stewart, "LINPACK
//      Users' Guide", SIAM, 1979. determinant()'s balanced accumulation follows
//      the convention of its xGEDI routines, which return determinants as a
//      (fraction, exponent) pair to avoid spurious overflow/underflow.
//  [5] P. H. Sterbenz, "Floating-Point Computation", Prentice-Hall, 1974.
//      Scaling by a power of two is exact, the property the balanced
//      accumulation relies on.
 
#ifndef EIGEN_STRUCTURED_CIRCULANT_H
#define EIGEN_STRUCTURED_CIRCULANT_H
 
// IWYU pragma: private
#include "./InternalHeaderCheck.h"
 
namespace Eigen {
 
template <typename Scalar_, int Size_ = Dynamic>
class Circulant;
 
namespace internal {
 
template <typename Scalar_, int Size_>
struct traits<Circulant<Scalar_, Size_>> {
  using Scalar = Scalar_;
  using StorageKind = Dense;
  using XprKind = MatrixXpr;
  using StorageIndex = int;
  static constexpr int RowsAtCompileTime = Size_;
  static constexpr int ColsAtCompileTime = Size_;
  static constexpr int MaxRowsAtCompileTime = Size_;
  static constexpr int MaxColsAtCompileTime = Size_;
  // 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 int Flags = 0;
};
 
template <typename Scalar_, int Size_>
struct evaluator_traits<Circulant<Scalar_, Size_>> {
  using Kind = IndexBased;
  using Shape = StructuredShape;
};
 
}  // namespace internal

template <typename Scalar_, int Size_>
class Circulant : public EigenBase<Circulant<Scalar_, Size_>> {
 public:
  using Scalar = Scalar_;
  using RealScalar = typename NumTraits<Scalar>::Real;
  using StorageIndex = int;
  using Complex = std::complex<RealScalar>;
  using GeneratorType = Matrix<Scalar, Size_, 1>;
  using ComplexVector = Matrix<Complex, Dynamic, 1>;
  using RealVector = Matrix<RealScalar, Size_, 1>;
  using ComplexMatrix = Matrix<Complex, Size_, Size_>;
 
  static constexpr int RowsAtCompileTime = Size_;
  static constexpr int ColsAtCompileTime = Size_;
  static constexpr int MaxRowsAtCompileTime = Size_;
  static constexpr int MaxColsAtCompileTime = Size_;
  static constexpr int SizeAtCompileTime = internal::size_at_compile_time(Size_, Size_);
  static constexpr int MaxSizeAtCompileTime = SizeAtCompileTime;
  static constexpr bool IsRowMajor = false;
  // Deliberately no IsVectorAtCompileTime: Ref<const Circulant>'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 Derived>
  explicit Circulant(const MatrixBase<Derived>& col) : m_col(col) {
    EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
    eigen_assert(m_col.size() > 0 && "Circulant generator must be non-empty");
    if (m_col.size() > internal::structured_direct_threshold()) {
      m_symbol = computeSymbol();
      m_prodSymbol = computeProdSymbol(m_col);
    }
  }
 
  EIGEN_DEVICE_FUNC Index rows() const { return m_col.size(); }
  EIGEN_DEVICE_FUNC Index cols() const { return m_col.size(); }

  const GeneratorType& column() const { return m_col; }

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

  Scalar coeff(Index row, Index col) const {
    Index k = row - col;
    if (k < 0) k += rows();
    return m_col.coeff(k);
  }

  Circulant transpose() const {
    const Index n = rows();
    GeneratorType col(n);
    col[0] = m_col[0];
    if (n > 1) col.tail(n - 1) = m_col.tail(n - 1).reverse();
    return Circulant(col, internal::structured_reverse_symbol(m_symbol),
                     internal::structured_reverse_symbol(m_prodSymbol));
  }

  Circulant conjugate() const {
    return Circulant(m_col.conjugate(), internal::structured_reverse_symbol(m_symbol).conjugate(),
                     internal::structured_reverse_symbol(m_prodSymbol).conjugate());
  }

  Circulant adjoint() const {
    const Index n = rows();
    GeneratorType col(n);
    col[0] = numext::conj(m_col[0]);
    if (n > 1) col.tail(n - 1) = m_col.tail(n - 1).reverse().conjugate();
    return Circulant(col, m_symbol.conjugate(), m_prodSymbol.conjugate());
  }

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

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

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

  template <typename Rhs>
  Product<Circulant, 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<Circulant, Rhs>(*this, x.derived());
  }

  template <typename Rhs>
  Matrix<Scalar, Size_, 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)
    const Index n = rows();
    eigen_assert(b.rows() == n && "right-hand side has the wrong number of rows");
    const ComplexVector s = symbol();
    const RealVector mods = s.cwiseAbs();
    const ComplexVector sinv = internal::structured_pinv_symbol(s, mods, internal::structured_rank_threshold(s, mods));
    Matrix<Scalar, Size_, Rhs::ColsAtCompileTime> x(n, b.cols());
    x.setZero();
    internal::structured_fft_apply(x, sinv, n, b.derived(), Scalar(1));
    return x;
  }

  Index rank() const {
    const ComplexVector s = symbol();
    const RealVector mods = s.cwiseAbs();
    return (!(mods.array() < internal::structured_rank_threshold(s, mods))).count();  // NaN entries count as non-zero
  }

  Circulant inverse() const {
    const Index n = rows();
    const ComplexVector sinv = symbol().cwiseInverse();
    GeneratorType col(n);
    if (n == 1) {
      col = internal::structured_scalar_part_impl<Scalar>::run(sinv);
    } else {
      auto&& fft = internal::structured_fft_engine<RealScalar>();
      ComplexVector ct(n);
      fft.inv(ct, sinv, n);
      col = internal::structured_scalar_part_impl<Scalar>::run(ct);
    }
    return Circulant(col, n > internal::structured_direct_threshold() ? sinv : ComplexVector(),
                     n > internal::structured_direct_threshold() ? computeProdSymbol(col) : ComplexVector());
  }

  Scalar determinant() const {
    const ComplexVector s = symbol();
    Complex det(1);
    Index exponent = 0;
    for (Index k = 0; k < s.size(); ++k)
      det = internal::structured_balance(Complex(det * internal::structured_balance(s[k], exponent)), exponent);
    return internal::structured_scalar_part_impl<Scalar>::run_scalar(internal::structured_ldexp_clamped(det, exponent));
  }

  ComplexVector eigenvalues() const { return symbol(); }

  ComplexMatrix eigenvectors() const {
    const Index n = rows();
    const ComplexVector roots = fourierRoots();
    ComplexMatrix F(n, n);
    for (Index k = 0; k < n; ++k) fourierColumn(F, roots, k, k);
    return F;
  }

  RealVector singularValues() const {
    ComplexVector s;
    RealVector mods;
    const std::vector<Index> perm = svdOrdering(s, mods);
    RealVector sv(s.size());
    for (Index t = 0; t < s.size(); ++t) sv[t] = mods[perm[t]];
    return sv;
  }

  ComplexMatrix matrixU() const {
    const Index n = rows();
    ComplexVector s;
    RealVector mods;
    const std::vector<Index> perm = svdOrdering(s, mods);
    const ComplexVector roots = fourierRoots();
    ComplexMatrix U(n, n);
    for (Index t = 0; t < n; ++t) {
      fourierColumn(U, roots, perm[t], t);
      const RealScalar a = mods[perm[t]];
      if (a > RealScalar(0)) U.col(t) *= s[perm[t]] / a;
    }
    return U;
  }

  ComplexMatrix matrixV() const {
    const Index n = rows();
    ComplexVector s;
    RealVector mods;
    const std::vector<Index> perm = svdOrdering(s, mods);
    const ComplexVector roots = fourierRoots();
    ComplexMatrix V(n, n);
    for (Index t = 0; t < n; ++t) fourierColumn(V, roots, perm[t], t);
    return V;
  }

  template <typename Dest, typename Rhs, typename ProductScalar>
  void addProduct(Dest& dst, const Rhs& rhs, const ProductScalar& alpha) const {
    const Index n = rows();
    eigen_assert(rhs.rows() == n && "invalid product: dimensions do not match");
    if (n <= internal::structured_direct_threshold()) {
      directProduct(dst, rhs, alpha);
      return;
    }
    // Products use the padded embedding symbol when the operator size is not
    // 5-smooth (see computeProdSymbol()); the exact-size transform of such a
    // size runs through kissfft's quadratic generic butterfly.
    const ComplexVector& s = m_prodSymbol.size() > 0 ? m_prodSymbol : m_symbol;
    internal::structured_fft_apply(dst, s, n, rhs, alpha);
  }
 
 private:
  Circulant(const GeneratorType& col, const ComplexVector& symbol, const ComplexVector& prodSymbol)
      : m_col(col), m_symbol(symbol), m_prodSymbol(prodSymbol) {}

  static ComplexVector computeProdSymbol(const GeneratorType& col) {
    const Index n = col.size();
    if (internal::fft_next_good_size(n) == n) return ComplexVector();
    const Index p = internal::fft_next_good_size(2 * n - 1);
    ComplexVector embedding = ComplexVector::Zero(p);
    embedding.head(n) = col.template cast<Complex>();
    embedding.tail(n - 1) = col.tail(n - 1).template cast<Complex>();
    ComplexVector symbol(p);
    auto&& fft = internal::structured_fft_engine<RealScalar>();
    fft.fwd(symbol, embedding, p);
    return symbol;
  }

  template <typename Dest, typename Rhs, typename ProductScalar>
  void directProductColumn(Dest& dst, const Rhs& rhs, Index k, const ProductScalar& alpha) const {
    const Index n = rows();
    // 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 <= 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 < n; ++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 c rotated downward by j; split it into two contiguous segments.
    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.head(j) += xj * m_col.tail(j);
      dstCol.tail(n - j) += xj * m_col.head(n - 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);
  }

  std::vector<Index> svdOrdering(ComplexVector& s, RealVector& mods) const {
    const Index n = rows();
    s = symbol();
    EIGEN_IF_CONSTEXPR (!NumTraits<Scalar>::IsComplex) {
      const Index pairs = (n - 1) / 2;
      mods.resize(n);
      mods.head(n - pairs) = s.head(n - pairs).cwiseAbs();
      mods.tail(pairs) = mods.segment(1, pairs).reverse();
    } else {
      mods = s.cwiseAbs();
    }
    return internal::structured_svd_permutation(mods);
  }
 
  // Cache the n distinct roots once per Fourier matrix, preserving polar()'s
  // bounded arguments and rounding while avoiding n^2 trigonometric evaluations.
  ComplexVector fourierRoots() const {
    const Index n = rows();
    const RealScalar scale = RealScalar(1) / numext::sqrt(RealScalar(n));
    ComplexVector roots(n);
    for (Index j = 0; j < n; ++j)
      roots[j] = std::polar(scale, RealScalar(2 * EIGEN_PI) * RealScalar(j) / RealScalar(n));
    return roots;
  }
 
  // Column k uses roots[j*k mod n]; incremental modular indexing avoids j*k overflow.
  void fourierColumn(ComplexMatrix& F, const ComplexVector& roots, Index k, Index dstCol) const {
    const Index n = rows();
    Index jk = 0;
    for (Index j = 0; j < n; ++j) {
      F(j, dstCol) = roots[jk];
      jk += k;
      if (jk >= n) jk -= n;
    }
  }

  ComplexVector computeSymbol() const {
    const Index n = m_col.size();
    const ComplexVector cc = m_col.template cast<Complex>();
    if (n == 1) return cc;  // the DFT of a single sample is the identity
    ComplexVector symbol(n);
    auto&& fft = internal::structured_fft_engine<RealScalar>();
    fft.fwd(symbol, cc, n);
    return symbol;
  }
 
  GeneratorType m_col;
  ComplexVector m_symbol;
  // The padded embedding symbol products use when the size is not 5-smooth;
  // empty otherwise. See computeProdSymbol().
  ComplexVector m_prodSymbol;
};

template <typename Derived>
Circulant<typename Derived::Scalar, Derived::SizeAtCompileTime> makeCirculant(const MatrixBase<Derived>& col) {
  return Circulant<typename Derived::Scalar, Derived::SizeAtCompileTime>(col);
}
 
namespace internal {
 
// StructuredShape avoids ambiguity with the stock dense product specializations.
template <typename Scalar_, int Size_, typename Rhs, int ProductTag>
struct generic_product_impl<Circulant<Scalar_, Size_>, Rhs, StructuredShape, DenseShape, ProductTag>
    : structured_product_impl<Circulant<Scalar_, Size_>, Rhs> {};
 
}  // namespace internal
 
}  // namespace Eigen
 
#endif  // EIGEN_STRUCTURED_CIRCULANT_H