HessenbergDecomposition.h

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr>
// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
//
// 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-License-Identifier: MPL-2.0
 
#ifndef EIGEN_HESSENBERGDECOMPOSITION_H
#define EIGEN_HESSENBERGDECOMPOSITION_H
 
// IWYU pragma: private
#include "./InternalHeaderCheck.h"
 
namespace Eigen {
 
namespace internal {
 
template <typename MatrixType>
struct HessenbergDecompositionMatrixHReturnType;
template <typename MatrixType>
struct traits<HessenbergDecompositionMatrixHReturnType<MatrixType>> {
  using ReturnType = typename MatrixType::PlainObject;
};
 
template <typename MatrixType, typename CoeffVectorType, typename WorkspaceType>
void hessenberg_decomposition_inplace(MatrixType& matA, CoeffVectorType& hCoeffs, WorkspaceType& temp);
 
}  // namespace internal

template <typename MatrixType_>
class HessenbergDecomposition {
 public:
  using MatrixType = MatrixType_;
 
  enum {
    Size = MatrixType::RowsAtCompileTime,
    SizeMinusOne = Size == Dynamic ? Dynamic : Size - 1,
    Options = internal::plain_object_options<MatrixType>::value,
    MaxSize = MatrixType::MaxRowsAtCompileTime,
    MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : MaxSize - 1
  };

  using Scalar = typename MatrixType::Scalar;
  using Index = Eigen::Index;  

  using CoeffVectorType = Matrix<Scalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1>;

  using HouseholderSequenceType =
      HouseholderSequence<MatrixType, internal::remove_all_t<typename CoeffVectorType::ConjugateReturnType>>;
 
  using MatrixHReturnType = internal::HessenbergDecompositionMatrixHReturnType<MatrixType>;

  explicit HessenbergDecomposition(Index size = Size == Dynamic ? 2 : Size)
      : m_matrix(size, size), m_temp(size), m_isInitialized(false) {
    if (size > 1) m_hCoeffs.resize(size - 1);
  }

  template <typename InputType>
  explicit HessenbergDecomposition(const EigenBase<InputType>& matrix)
      : m_matrix(matrix.derived()), m_temp(matrix.rows()), m_isInitialized(false) {
    computeInPlace();
  }

  template <typename InputType>
  explicit HessenbergDecomposition(EigenBase<InputType>& matrix)
      : m_matrix(matrix.derived()), m_temp(matrix.rows()), m_isInitialized(false) {
    computeInPlace();
  }

  template <typename InputType>
  HessenbergDecomposition& compute(const EigenBase<InputType>& matrix) {
    m_matrix = matrix.derived();
    computeInPlace();
    return *this;
  }

  const CoeffVectorType& householderCoefficients() const {
    eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
    return m_hCoeffs;
  }

  const MatrixType& packedMatrix() const {
    eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
    return m_matrix;
  }

  HouseholderSequenceType matrixQ() const {
    eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
    return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate()).setLength(m_matrix.rows() - 1).setShift(1);
  }

  MatrixHReturnType matrixH() const {
    eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
    return MatrixHReturnType(*this);
  }
 
 private:
  using VectorType = Matrix<Scalar, 1, Size, int(Options) | int(RowMajor), 1, MaxSize>;
  using RealScalar = typename NumTraits<Scalar>::Real;
 
  void computeInPlace() {
    if (m_matrix.rows() < 2) {
      m_isInitialized = true;
      return;
    }
    m_hCoeffs.resize(m_matrix.rows() - 1, 1);
    internal::hessenberg_decomposition_inplace(m_matrix, m_hCoeffs, m_temp);
    m_isInitialized = true;
  }
 
 protected:
  MatrixType m_matrix;
  CoeffVectorType m_hCoeffs;
  VectorType m_temp;
  bool m_isInitialized;
};
 
namespace internal {

template <typename MatrixType, typename CoeffVectorType, typename WorkspaceType>
void hessenberg_decomposition_inplace(MatrixType& matA, CoeffVectorType& hCoeffs, WorkspaceType& temp) {
  using Scalar = typename MatrixType::Scalar;
  using RealScalar = typename NumTraits<Scalar>::Real;
  eigen_assert(matA.rows() == matA.cols());
  Index n = matA.rows();
  temp.resize(n);
  for (Index i = 0; i < n - 1; ++i) {
    // let's consider the vector v = i-th column starting at position i+1
    Index remainingSize = n - i - 1;
    RealScalar beta;
    Scalar h;
    auto householder = matA.col(i).tail(remainingSize);
    const RealScalar tailSqNorm =
        remainingSize == 1 ? RealScalar(0) : householder.tail(remainingSize - 1).unwind().squaredNorm();
    const RealScalar tol = (std::numeric_limits<RealScalar>::min)();
    // Preserve negligible subdiagonal entries instead of rotating them into much larger matrix coefficients.  The
    // latter can erase small eigenvalues through cancellation even though the reflector itself is accurate.
    if (tailSqNorm <= tol && numext::abs2(numext::imag(householder.coeff(0))) <= tol) {
      h = Scalar(0);
      beta = numext::real(householder.coeff(0));
      householder.tail(remainingSize - 1).setZero();
    } else {
      householder.makeHouseholderInPlace(h, beta);
    }
    matA.col(i).coeffRef(i + 1) = beta;
    hCoeffs.coeffRef(i) = h;
 
    // Apply similarity transformation to remaining columns,
    // i.e., compute A = H A H'
 
    // A = H A
    matA.bottomRightCorner(remainingSize, remainingSize)
        .applyHouseholderOnTheLeft(matA.col(i).tail(remainingSize - 1), h, &temp.coeffRef(0));
 
    // A = A H'
    matA.rightCols(remainingSize)
        .applyHouseholderOnTheRight(matA.col(i).tail(remainingSize - 1), numext::conj(h), &temp.coeffRef(0));
  }
}

template <typename MatrixType, typename CoeffVectorType, typename WorkspaceType, typename MatrixQType>
void hessenberg_decomposition_inplace(MatrixType& matA, CoeffVectorType& hCoeffs, WorkspaceType& temp,
                                      MatrixQType& matQ, bool computeQ) {
  using HouseholderSequenceType =
      HouseholderSequence<MatrixType, remove_all_t<typename CoeffVectorType::ConjugateReturnType>>;
  const Index n = matA.rows();
  hCoeffs.resize(n - 1);
  hessenberg_decomposition_inplace(matA, hCoeffs, temp);
  if (computeQ) {
    Index firstReflector = 0;
    while (firstReflector < n - 1 && numext::is_exactly_zero(hCoeffs.coeff(firstReflector))) ++firstReflector;
    if (firstReflector == n - 1) {
      // All tau_i = 0: Q = I, without allocating block Householder factors.
      matQ.setIdentity(n, n);
    } else {
      HouseholderSequenceType(matA, hCoeffs.conjugate()).setLength(n - 1).setShift(1).evalTo(matQ, temp);
    }
  }
  if (n > 2) matA.bottomLeftCorner(n - 2, n - 2).template triangularView<Lower>().setZero();
}

template <typename MatrixType>
struct HessenbergDecompositionMatrixHReturnType
    : public ReturnByValue<HessenbergDecompositionMatrixHReturnType<MatrixType>> {
 public:
  HessenbergDecompositionMatrixHReturnType(const HessenbergDecomposition<MatrixType>& hess) : m_hess(hess) {}

  template <typename ResultType>
  inline void evalTo(ResultType& result) const {
    result = m_hess.packedMatrix();
    Index n = result.rows();
    if (n > 2) result.bottomLeftCorner(n - 2, n - 2).template triangularView<Lower>().setZero();
  }
 
  Index rows() const { return m_hess.packedMatrix().rows(); }
  Index cols() const { return m_hess.packedMatrix().cols(); }
 
 protected:
  const HessenbergDecomposition<MatrixType>& m_hess;
};
 
}  // end namespace internal
 
}  // end namespace Eigen
 
#endif  // EIGEN_HESSENBERGDECOMPOSITION_H