Re: [eigen] Eigenvalues of (lower) Hessenberg matrix

[Yixuan Qiu <yixuanq@xxxxxxxxx>]

Hi Ian,

Given that your matrices are at very small scale, I don't expect my class will have visible improvement. Actually you found it to be slower, which I think may be due to several reasons:

1. My class does some basic scaling of matrix before factorization.

2. Did you select "computeEigenvectors = false" in EigenSolver? My class assumes computing eigenvectors, so some calculations are unneeded in your case. You can try the reduced version attached.

3. Did you randomize the order of execution? Did you have enough replications of experiment? Was the difference statistically significant relative to the standard error? Personally I don't fully trust benchmark results measured by elapsed time at the scale of microseconds. I would recommend Callgrind as the profiler for such tasks.

Hope this helps.

Best,

Yixuan

2017-04-04 15:40 GMT-04:00 Ian Bell <ian.h.bell@xxxxxxxxx>:

Sadly, it seems to be consistently slower than the naive eigenvalues method of the MatrixXd class (see attached).  Any idea why?  Have you done any performance benchmarking with the UpperHessenberg code?

On Mon, Apr 3, 2017 at 8:16 PM, Ian Bell <ian.h.bell@xxxxxxxxx> wrote:

Amazing! I'll try that tomorrow at work.  From your testing, how long approximately would it take to find the eigenvalues of a 10x10 Hessenberg matrix?

On Mon, Apr 3, 2017 at 8:00 PM, Yixuan Qiu <yixuanq@xxxxxxxxx> wrote:

Eigen implements the upper Hessenberg eigen solver internally in the EigenSolver class.

I have extracted the relevant code and created an independent class in my Spectra library, so probably this is what you want: https://github.com/yixuan/spectra/blob/master/include/LinAlg/UpperHessenbergEigen..h

You can drop the namespace declaration if you want.

The usage of this class is similar to other eigen solvers in Eigen: https://github.com/yixuan/spectra/blob/master/test/Eigen.cpp#L27-L29

Best,

Yixuan

2017-04-03 20:59 GMT-04:00 Ian Bell <ian.h.bell@xxxxxxxxx>:

I have a matrix, that by its construction is known to be Hessenberg (rigorously, lower Hessenberg, but it doesn't matter because the transpose of a matrix has the same eigenvalues as the original matrix and all I care about is eigenvalues).  Is there any magic trick in Eigen that allows for more efficient evaluation of eigenvalues?  The standard method eigenvalues() doesn't seem to do anything too smart about checking for the Hessenberg-ness of the matrix.  Lapack has the function dhseqr, is there anything similar in Eigen?

Ian

--

Yixuan Qiu <yixuanq@xxxxxxxxx>
Department of Statistics,
Purdue University

--

Yixuan Qiu <yixuanq@xxxxxxxxx>
Department of Statistics,
Purdue University

// The code was adapted from Eigen/src/Eigenvaleus/EigenSolver.h
//
// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@xxxxxxxx>
// Copyright (C) 2010,2012 Jitse Niesen <jitse@xxxxxxxxxxxxxxxxx>
// Copyright (C) 2016-2017 Yixuan Qiu <yixuan.qiu@xxxxxxxx>
//
// 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/.

#ifndef UPPER_HESSENBERG_EIGEN_VAL_ONLY_H
#define UPPER_HESSENBERG_EIGEN_VAL_ONLY_H

#include <Eigen/Core>
#include <Eigen/Eigenvalues>
#include <stdexcept>

namespace Spectra {


template <typename Scalar = double>
class UpperHessenbergEigenValOnly
{
private:
    typedef Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> Matrix;
    typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1> Vector;

    typedef typename Matrix::Index Index;

    typedef Eigen::Ref<Matrix> GenericMatrix;
    typedef const Eigen::Ref<const Matrix> ConstGenericMatrix;

    typedef std::complex<Scalar> Complex;
    typedef Eigen::Matrix<Complex, Eigen::Dynamic, Eigen::Dynamic> ComplexMatrix;
    typedef Eigen::Matrix<Complex, Eigen::Dynamic, 1> ComplexVector;

    Index m_n;                             // Size of the matrix
    Eigen::RealSchur<Matrix> m_realSchur;  // Schur decomposition solver
    Matrix m_matT;                         // Schur T matrix
    ComplexVector m_eivalues;              // Eigenvalues

    bool m_computed;

public:

    UpperHessenbergEigenValOnly() :
        m_n(0), m_computed(false)
    {}

    UpperHessenbergEigenValOnly(ConstGenericMatrix& mat) :
        m_n(mat.rows()), m_computed(false)
    {
        compute(mat);
    }

    void compute(ConstGenericMatrix& mat)
    {
        using std::abs;
        using std::sqrt;

        if(mat.rows() != mat.cols())
            throw std::invalid_argument("UpperHessenbergEigenValOnly: matrix must be square");

        m_n = mat.rows();
        // Scale matrix prior to the Schur decomposition
        const Scalar scale = mat.cwiseAbs().maxCoeff();

        // Reduce to real Schur form
        Matrix Q = Matrix::Identity(m_n, m_n);
        m_realSchur.computeFromHessenberg(mat / scale, Q, false);
        if(m_realSchur.info() != Eigen::Success)
            throw std::runtime_error("UpperHessenbergEigenValOnly: eigen decomposition failed");

        m_matT = m_realSchur.matrixT();

        // Compute eigenvalues from matT
        m_eivalues.resize(m_n);
        Index i = 0;
        while(i < m_n)
        {
            // Real eigenvalue
            if(i == m_n - 1 || m_matT.coeff(i+1, i) == Scalar(0))
            {
                m_eivalues.coeffRef(i) = m_matT.coeff(i, i);
                ++i;
            }
            else  // Complex eigenvalues
            {
                Scalar p = Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i+1, i+1));
                Scalar z;
                // Compute z = sqrt(abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1)));
                // without overflow
                {
                    Scalar t0 = m_matT.coeff(i+1, i);
                    Scalar t1 = m_matT.coeff(i, i+1);
                    Scalar maxval = std::max(abs(p), std::max(abs(t0), abs(t1)));
                    t0 /= maxval;
                    t1 /= maxval;
                    Scalar p0 = p / maxval;
                    z = maxval * sqrt(abs(p0 * p0 + t0 * t1));
                }
                m_eivalues.coeffRef(i)   = Complex(m_matT.coeff(i+1, i+1) + p, z);
                m_eivalues.coeffRef(i+1) = Complex(m_matT.coeff(i+1, i+1) + p, -z);
                i += 2;
            }
        }

        // Scale eigenvalues back
        m_eivalues *= scale;

        m_computed = true;
    }

    ComplexVector eigenvalues()
    {
        if(!m_computed)
            throw std::logic_error("UpperHessenbergEigenValOnly: need to call compute() first");

        return m_eivalues;
    }
};


} // namespace Spectra

#endif // UPPER_HESSENBERG_EIGEN_VAL_ONLY_H