Added a Hessenberg decomposition class for both real and complex matrices.

This is the first step towards a non-selfadjoint eigen solver.
Notes:
 - We might consider merging Tridiagonalization and Hessenberg toghether ?
 - Or we could factorize some code into a Householder class (could also be shared with QR)

Eigen/QR

@@ -9,6 +9,7 @@ namespace Eigen {
 #include "src/QR/Tridiagonalization.h"
 #include "src/QR/EigenSolver.h"
 #include "src/QR/SelfAdjointEigenSolver.h"
+#include "src/QR/HessenbergDecomposition.h"
 
 } // namespace Eigen
 

Eigen/src/QR/HessenbergDecomposition.h

@@ -0,0 +1,243 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_HESSENBERGDECOMPOSITION_H
+#define EIGEN_HESSENBERGDECOMPOSITION_H
+
+/** \class HessenbergDecomposition
+  *
+  * \brief Reduces a squared matrix to an Hessemberg form
+  *
+  * \param MatrixType the type of the matrix of which we are computing the Hessenberg decomposition
+  *
+  * This class performs an Hessenberg decomposition of a matrix \f$ A \f$ such that:
+  * \f$ A = Q H Q^* \f$ where \f$ Q \f$ is unitary and \f$ H \f$ a Hessenberg matrix.
+  *
+  * \sa class Tridiagonalization, class Qr
+  */
+template<typename _MatrixType> class HessenbergDecomposition
+{
+  public:
+
+    typedef _MatrixType MatrixType;
+    typedef typename MatrixType::Scalar Scalar;
+    typedef typename NumTraits<Scalar>::Real RealScalar;
+
+    enum {
+      Size = MatrixType::RowsAtCompileTime,
+      SizeMinusOne = MatrixType::RowsAtCompileTime==Dynamic
+                        ? Dynamic
+                        : MatrixType::RowsAtCompileTime-1};
+
+    typedef Matrix<Scalar, SizeMinusOne, 1> CoeffVectorType;
+    typedef Matrix<RealScalar, Size, 1> DiagonalType;
+    typedef Matrix<RealScalar, SizeMinusOne, 1> SubDiagonalType;
+
+    typedef typename NestByValue<DiagonalCoeffs<MatrixType> >::RealReturnType DiagonalReturnType;
+
+    typedef typename NestByValue<DiagonalCoeffs<
+        NestByValue<Block<
+          MatrixType,SizeMinusOne,SizeMinusOne> > > >::RealReturnType SubDiagonalReturnType;
+
+    HessenbergDecomposition()
+    {}
+
+    HessenbergDecomposition(int rows, int cols)
+      : m_matrix(rows,cols), m_hCoeffs(rows-1)
+    {}
+
+    HessenbergDecomposition(const MatrixType& matrix)
+      : m_matrix(matrix),
+        m_hCoeffs(matrix.cols()-1)
+    {
+      _compute(m_matrix, m_hCoeffs);
+    }
+
+    /** Computes or re-compute the Hessenberg decomposition for the matrix \a matrix.
+      *
+      * This method allows to re-use the allocated data.
+      */
+    void compute(const MatrixType& matrix)
+    {
+      m_matrix = matrix;
+      m_hCoeffs.resize(matrix.rows()-1,1);
+      _compute(m_matrix, m_hCoeffs);
+    }
+
+    /** \returns the householder coefficients allowing to
+      * reconstruct the matrix Q from the packed data.
+      *
+      * \sa packedMatrix()
+      */
+    CoeffVectorType householderCoefficients(void) const { return m_hCoeffs; }
+
+    /** \returns the internal result of the decomposition.
+      *
+      * The returned matrix contains the following information:
+      *  - the upper part and lower sub-diagonal represent the Hessenberg matrix H
+      *  - the rest of the lower part contains the Householder vectors that, combined with
+      *    Householder coefficients returned by householderCoefficients(),
+      *    allows to reconstruct the matrix Q as follow:
+      *       Q = H_{N-1} ... H_1 H_0
+      *    where the matrices H are the Householder transformation:
+      *       H_i = (I - h_i * v_i * v_i')
+      *    where h_i == householderCoefficients()[i] and v_i is a Householder vector:
+      *       v_i = [ 0, ..., 0, 1, M(i+2,i), ..., M(N-1,i) ]
+      *
+      * See LAPACK for further details on this packed storage.
+      */
+    const MatrixType& packedMatrix(void) const { return m_matrix; }
+
+    MatrixType matrixQ(void) const;
+    MatrixType matrixH(void) const;
+
+  private:
+
+    static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs);
+
+  protected:
+    MatrixType m_matrix;
+    CoeffVectorType m_hCoeffs;
+};
+
+
+/** \internal
+  * Performs a tridiagonal decomposition of \a matA in place.
+  *
+  * \param matA the input selfadjoint matrix
+  * \param hCoeffs returned Householder coefficients
+  *
+  * The result is written in the lower triangular part of \a matA.
+  *
+  * Implemented from Golub's "Matrix Computations", algorithm 8.3.1.
+  *
+  * \sa packedMatrix()
+  */
+template<typename MatrixType>
+void HessenbergDecomposition<MatrixType>::_compute(MatrixType& matA, CoeffVectorType& hCoeffs)
+{
+  assert(matA.rows()==matA.cols());
+  int n = matA.rows();
+  for (int i = 0; i<n-2; ++i)
+  {
+    // let's consider the vector v = i-th column starting at position i+1
+
+    // start of the householder transformation
+    // squared norm of the vector v skipping the first element
+    RealScalar v1norm2 = matA.col(i).end(n-(i+2)).norm2();
+
+    if (ei_isMuchSmallerThan(v1norm2,static_cast<Scalar>(1)))
+    {
+      hCoeffs.coeffRef(i) = 0.;
+    }
+    else
+    {
+      Scalar v0 = matA.col(i).coeff(i+1);
+      RealScalar beta = ei_sqrt(ei_abs2(v0)+v1norm2);
+      if (ei_real(v0)>=0.)
+        beta = -beta;
+      matA.col(i).end(n-(i+2)) *= (Scalar(1)/(v0-beta));
+      matA.col(i).coeffRef(i+1) = beta;
+      Scalar h = (beta - v0) / beta;
+      // end of the householder transformation
+
+      // Apply similarity transformation to remaining columns,
+      // i.e., A = H' A H where H = I - h v v' and v = matA.col(i).end(n-i-1)
+      matA.col(i).coeffRef(i+1) = 1;
+
+      // first let's do A = H A
+      matA.corner(BottomRight,n-i-1,n-i-1) -= ((ei_conj(h) * matA.col(i).end(n-i-1)) *
+        (matA.col(i).end(n-i-1).adjoint() * matA.corner(BottomRight,n-i-1,n-i-1)).lazy()).lazy();
+
+      // now let's do A = A H
+      matA.corner(BottomRight,n,n-i-1) -= ((matA.corner(BottomRight,n,n-i-1) * matA.col(i).end(n-i-1)).lazy() *
+        (h * matA.col(i).end(n-i-1).adjoint())).lazy();
+
+      matA.col(i).coeffRef(i+1) = beta;
+      hCoeffs.coeffRef(i) = h;
+    }
+  }
+  if (NumTraits<Scalar>::IsComplex)
+  {
+    // Householder transformation on the remaining single scalar
+    int i = n-2;
+    Scalar v0 = matA.coeff(i+1,i);
+
+    RealScalar beta = ei_sqrt(ei_abs2(v0));
+    if (ei_real(v0)>=0.)
+      beta = -beta;
+    Scalar h = (beta - v0) / beta;
+    hCoeffs.coeffRef(i) = h;
+
+    // A = H* A
+    matA.corner(BottomRight,n-i-1,n-i) -= ei_conj(h) * matA.corner(BottomRight,n-i-1,n-i);
+
+    // A = A H
+    matA.col(n-1) -= h * matA.col(n-1);
+  }
+  else
+  {
+    hCoeffs.coeffRef(n-2) = 0;
+  }
+}
+
+/** reconstructs and returns the matrix Q */
+template<typename MatrixType>
+typename HessenbergDecomposition<MatrixType>::MatrixType
+HessenbergDecomposition<MatrixType>::matrixQ(void) const
+{
+  int n = m_matrix.rows();
+  MatrixType matQ = MatrixType::identity(n,n);
+  for (int i = n-2; i>=0; i--)
+  {
+    Scalar tmp = m_matrix.coeff(i+1,i);
+    m_matrix.const_cast_derived().coeffRef(i+1,i) = 1;
+
+    matQ.corner(BottomRight,n-i-1,n-i-1) -=
+      ((m_hCoeffs.coeff(i) * m_matrix.col(i).end(n-i-1)) *
+      (m_matrix.col(i).end(n-i-1).adjoint() * matQ.corner(BottomRight,n-i-1,n-i-1)).lazy()).lazy();
+
+    m_matrix.const_cast_derived().coeffRef(i+1,i) = tmp;
+  }
+  return matQ;
+}
+
+/** constructs and returns the matrix H.
+  * Note that the matrix H is equivalent to the upper part of the packed matrix
+  * (including the lower sub-diagonal). Therefore, it might be often sufficient
+  * to directly use the packed matrix instead of creating a new one.
+  */
+template<typename MatrixType>
+typename HessenbergDecomposition<MatrixType>::MatrixType
+HessenbergDecomposition<MatrixType>::matrixH(void) const
+{
+  // FIXME should this function (and other similar) rather take a matrix as argument
+  // and fill it (avoids temporaries)
+  int n = m_matrix.rows();
+  MatrixType matH = m_matrix;
+  matH.corner(BottomLeft,n-2, n-2).template part<Lower>().setZero();
+  return matH;
+}
+
+#endif // EIGEN_HESSENBERGDECOMPOSITION_H

Eigen/src/QR/Tridiagonalization.h

@@ -29,7 +29,7 @@
   *
   * \brief Trigiagonal decomposition of a selfadjoint matrix
   *
-  * \param MatrixType the type of the matrix of which we are computing the eigen decomposition
+  * \param MatrixType the type of the matrix of which we are performing the tridiagonalization
   *
   * This class performs a tridiagonal decomposition of a selfadjoint matrix \f$ A \f$ such that:
   * \f$ A = Q T Q^* \f$ where \f$ Q \f$ is unitatry and \f$ T \f$ a real symmetric tridiagonal matrix
@@ -81,7 +81,7 @@ template<typename _MatrixType> class Tridiagonalization
     void compute(const MatrixType& matrix)
     {
       m_matrix = matrix;
-      m_hCoeffs.resize(matrix.rows()-1);
+      m_hCoeffs.resize(matrix.rows()-1, 1);
       _compute(m_matrix, m_hCoeffs);
     }
 
@@ -111,6 +111,7 @@ template<typename _MatrixType> class Tridiagonalization
     const MatrixType& packedMatrix(void) const { return m_matrix; }
 
     MatrixType matrixQ(void) const;
+    MatrixType matrixT(void) const;
     const DiagonalReturnType diagonal(void) const;
     const SubDiagonalReturnType subDiagonal(void) const;
 
@@ -252,6 +253,25 @@ Tridiagonalization<MatrixType>::subDiagonal(void) const
     .nestByValue().diagonal().nestByValue().real();
 }
 
+/** constructs and returns the tridiagonal matrix T.
+  * Note that the matrix T is equivalent to the diagonal and sub-diagonal of the packed matrix.
+  * Therefore, it might be often sufficient to directly use the packed matrix, or the vector
+  * expressions returned by diagonal() and subDiagonal() instead of creating a new matrix.
+  */
+template<typename MatrixType>
+typename Tridiagonalization<MatrixType>::MatrixType
+Tridiagonalization<MatrixType>::matrixT(void) const
+{
+  // FIXME should this function (and other similar) rather take a matrix as argument
+  // and fill it (avoids temporaries)
+  int n = m_matrix.rows();
+  MatrixType matT = m_matrix;
+  matT.corner(TopRight,n-1, n-1).diagonal() = subDiagonal().conjugate();
+  matT.corner(TopRight,n-2, n-2).template part<Upper>().setZero();
+  matT.corner(BottomLeft,n-2, n-2).template part<Lower>().setZero();
+  return matT;
+}
+
 /** Performs a full decomposition in place */
 template<typename MatrixType>
 void Tridiagonalization<MatrixType>::decomposeInPlace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)

test/eigensolver.cpp

@@ -54,9 +54,11 @@ template<typename MatrixType> void eigensolver(const MatrixType& m)
 void test_eigensolver()
 {
   for(int i = 0; i < 1; i++) {
+    // very important to test a 3x3 matrix since we provide a special path for it
     CALL_SUBTEST( eigensolver(Matrix3f()) );
     CALL_SUBTEST( eigensolver(Matrix4d()) );
     CALL_SUBTEST( eigensolver(MatrixXd(7,7)) );
     CALL_SUBTEST( eigensolver(MatrixXcd(6,6)) );
+    CALL_SUBTEST( eigensolver(MatrixXcd(3,3)) );
   }
 }

test/nomalloc.cpp

@@ -22,7 +22,9 @@
 // License and a copy of the GNU General Public License along with
 // Eigen. If not, see <http://www.gnu.org/licenses/>.
 
+// this hack is needed to make this file compiles with -pedantic (gcc)
 #define throw(X)
+// discard vectorization since operator new is not called in that case
 #define EIGEN_DONT_VECTORIZE 1
 #include "main.h"
 

test/qr.cpp

@@ -39,17 +39,31 @@ template<typename MatrixType> void qr(const MatrixType& m)
 
   MatrixType a = MatrixType::random(rows,cols);
   QR<MatrixType> qrOfA(a);
-
   VERIFY_IS_APPROX(a, qrOfA.matrixQ() * qrOfA.matrixR());
   VERIFY_IS_NOT_APPROX(a+MatrixType::identity(rows, cols), qrOfA.matrixQ() * qrOfA.matrixR());
+
+  SquareMatrixType b = a.adjoint() * a;
+
+  // check tridiagonalization
+  Tridiagonalization<SquareMatrixType> tridiag(b);
+  VERIFY_IS_APPROX(b, tridiag.matrixQ() * tridiag.matrixT() * tridiag.matrixQ().adjoint());
+
+  // check hessenberg decomposition
+  HessenbergDecomposition<SquareMatrixType> hess(b);
+  VERIFY_IS_APPROX(b, hess.matrixQ() * hess.matrixH() * hess.matrixQ().adjoint());
+  VERIFY_IS_APPROX(tridiag.matrixT(), hess.matrixH());
+  b = SquareMatrixType::random(cols,cols);
+  hess.compute(b);
+  VERIFY_IS_APPROX(b, hess.matrixQ() * hess.matrixH() * hess.matrixQ().adjoint());
 }
 
 void test_qr()
 {
   for(int i = 0; i < 1; i++) {
     CALL_SUBTEST( qr(Matrix2f()) );
-    CALL_SUBTEST( qr(Matrix3d()) );
+    CALL_SUBTEST( qr(Matrix4d()) );
     CALL_SUBTEST( qr(MatrixXf(12,8)) );
-//     CALL_SUBTEST( qr(MatrixXcd(17,7)) ); // complex numbers are not supported yet
+    CALL_SUBTEST( qr(MatrixXcd(5,5)) );
+    CALL_SUBTEST( qr(MatrixXcd(7,3)) );
   }
 }