SymEigen.h

Go to the documentation of this file.

// -*- c++ -*-

// Copyright (C) 2005,2009 Tom Drummond (twd20@cam.ac.uk)
//
// This file is part of the TooN Library.       This library is free
// software; 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, or (at your option)
// any later version.

// This library 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 General Public License for more details.

// You should have received a copy of the GNU General Public License along
// with this library; see the file COPYING.     If not, write to the Free
// Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307,
// USA.

// As a special exception, you may use this file as part of a free software
// library without restriction. Specifically, if other files instantiate
// templates or use macros or inline functions from this file, or you compile
// this file and link it with other files to produce an executable, this
// file does not by itself cause the resulting executable to be covered by
// the GNU General Public License.      This exception does not however
// invalidate any other reasons why the executable file might be covered by
// the GNU General Public License.

#ifndef __SYMEIGEN_H
#define __SYMEIGEN_H

#include <iostream>
#include <cassert>
#include <cmath>
#include <complex>
#include <TooN/lapack.h>

#include <TooN/TooN.h>

namespace TooN {
static const double root3=1.73205080756887729352744634150587236694280525381038062805580;

namespace Internal{

        using std::swap;

        static const double symeigen_condition_no=1e9;

        template <int Size> struct ComputeSymEigen {

                template<int Rows, int Cols, typename P, typename B>
                static inline void compute(const Matrix<Rows,Cols,P, B>& m, Matrix<Size,Size,P> & evectors, Vector<Size, P>& evalues) {

                        SizeMismatch<Rows, Cols>::test(m.num_rows(), m.num_cols());      //m must be square
                        SizeMismatch<Size, Rows>::test(m.num_rows(), evalues.size()); //m must be the size of the system
                        

                        evectors = m;
                        int N = evalues.size();
                        int lda = evalues.size();
                        int info;
                        int lwork=-1;
                        P size;

                        // find out how much space fortran needs
                        syev_((char*)"V",(char*)"U",&N,&evectors[0][0],&lda,&evalues[0], &size,&lwork,&info);
                        lwork = int(size);
                        Vector<Dynamic, P> WORK(lwork);

                        // now compute the decomposition
                        syev_((char*)"V",(char*)"U",&N,&evectors[0][0],&lda,&evalues[0], &WORK[0],&lwork,&info);

                        if(info!=0){
                                std::cerr << "In SymEigen<"<<Size<<">: " << info
                                                << " off-diagonal elements of an intermediate tridiagonal form did not converge to zero." << std::endl
                                                << "M = " << m << std::endl;
                        }
                }
        };

        template <> struct ComputeSymEigen<2> {

                template<typename P, typename B>
                static inline void compute(const Matrix<2,2,P,B>& m, Matrix<2,2,P>& eig, Vector<2, P>& ev) {
                        double trace = m[0][0] + m[1][1];
                        double det = m[0][0]*m[1][1] - m[0][1]*m[1][0];
                        double disc = trace*trace - 4 * det;
                        assert(disc>=0);
                        using std::sqrt;
                        double root_disc = sqrt(disc);
                        ev[0] = 0.5 * (trace - root_disc);
                        ev[1] = 0.5 * (trace + root_disc);
                        double a = m[0][0] - ev[0];
                        double b = m[0][1];
                        double magsq = a*a + b*b;
                        if (magsq == 0) {
                                eig[0][0] = 1.0;
                                eig[0][1] = 0;
                        } else {
                                eig[0][0] = -b;
                                eig[0][1] = a;
                                eig[0] *= 1.0/sqrt(magsq);
                        }
                        eig[1][0] = -eig[0][1];
                        eig[1][1] = eig[0][0];
                }
        };

        template <> struct ComputeSymEigen<3> {

                template<typename P, typename B>
                static inline void compute(const Matrix<3,3,P,B>& m, Matrix<3,3,P>& eig, Vector<3, P>& ev) {
            //method uses closed form solution of cubic equation to obtain roots of characteristic equation.

            //Polynomial terms of |a - l * Identity|
            //l^3 + a*l^2 + b*l + c

            const double& a11 = m[0][0];
            const double& a12 = m[0][1];
            const double& a13 = m[0][2];

            const double& a22 = m[1][1];
            const double& a23 = m[1][2];

            const double& a33 = m[2][2];

            //From matlab:
            double a = -a11-a22-a33;
            double b = a11*a22+a11*a33+a22*a33-a12*a12-a13*a13-a23*a23;
            double c = a11*(a23*a23)+(a13*a13)*a22+(a12*a12)*a33-a12*a13*a23*2.0-a11*a22*a33;

            //Using Cardano's method:
            double p = b - a*a/3;
            double q = c + (2*a*a*a - 9*a*b)/27;

            double alpha = -q/2;
            double beta_descriminant = q*q/4 + p*p*p/27;

            //beta_descriminant <= 0 for real roots!
            double beta = sqrt(-beta_descriminant);
            double r2 = alpha*alpha  - beta_descriminant;


            double cuberoot_r = pow(r2, 1./6);
            double theta3 = atan2(beta, alpha)/3;

            double A_plus_B = 2*cuberoot_r*cos(theta3);
            double A_minus_B = -2*cuberoot_r*sin(theta3);

            //calculate eigenvalues
            ev =  makeVector(A_plus_B, -A_plus_B/2 + A_minus_B * sqrt(3)/2, -A_plus_B/2 - A_minus_B * sqrt(3)/2) - Ones * a/3;

            if(ev[0] > ev[1])
                swap(ev[0], ev[1]);
            if(ev[1] > ev[2])
                swap(ev[1], ev[2]);
            if(ev[0] > ev[1])
                swap(ev[0], ev[1]);

            //calculate the eigenvectors
            eig[0][0]=a12 * a23 - a13 * (a22 - ev[0]);
            eig[0][1]=a12 * a13 - a23 * (a11 - ev[0]);
            eig[0][2]=(a11-ev[0])*(a22-ev[0]) - a12*a12;
            normalize(eig[0]);
            eig[1][0]=a12 * a23 - a13 * (a22 - ev[1]);
            eig[1][1]=a12 * a13 - a23 * (a11 - ev[1]);
            eig[1][2]=(a11-ev[1])*(a22-ev[1]) - a12*a12;
            normalize(eig[1]);
            eig[2][0]=a12 * a23 - a13 * (a22 - ev[2]);
            eig[2][1]=a12 * a13 - a23 * (a11 - ev[2]);
            eig[2][2]=(a11-ev[2])*(a22-ev[2]) - a12*a12;
            normalize(eig[2]);
                }
        };

};

template <int Size=Dynamic, typename Precision = double>
class SymEigen {
public:
        inline SymEigen(){}

        inline SymEigen(int m) : my_evectors(m,m), my_evalues(m) {}

        template<int R, int C, typename B>
        inline SymEigen(const Matrix<R, C, Precision, B>& m) : my_evectors(m.num_rows(), m.num_cols()), my_evalues(m.num_rows()) {
                compute(m);
        }

        template<int R, int C, typename B>
        inline void compute(const Matrix<R,C,Precision,B>& m){
                SizeMismatch<R, C>::test(m.num_rows(), m.num_cols());
                SizeMismatch<R, Size>::test(m.num_rows(), my_evectors.num_rows());
                Internal::ComputeSymEigen<Size>::compute(m, my_evectors, my_evalues);
        }

        template <int S, typename P, typename B>
        Vector<Size, Precision> backsub(const Vector<S,P,B>& rhs) const {
                return (my_evectors.T() * diagmult(get_inv_diag(Internal::symeigen_condition_no),(my_evectors * rhs)));
        }

        template <int R, int C, typename P, typename B>
        Matrix<Size,C, Precision> backsub(const Matrix<R,C,P,B>& rhs) const {
                return (my_evectors.T() * diagmult(get_inv_diag(Internal::symeigen_condition_no),(my_evectors * rhs)));
        }

        Matrix<Size, Size, Precision> get_pinv(const double condition=Internal::symeigen_condition_no) const {
                return my_evectors.T() * diagmult(get_inv_diag(condition),my_evectors);
        }

        Vector<Size, Precision> get_inv_diag(const double condition) const {
                Precision max_diag = -my_evalues[0] > my_evalues[my_evalues.size()-1] ? -my_evalues[0]:my_evalues[my_evalues.size()-1];
                Vector<Size, Precision> invdiag(my_evalues.size());
                for(int i=0; i<my_evalues.size(); i++){
                        if(fabs(my_evalues[i]) * condition > max_diag) {
                                invdiag[i] = 1/my_evalues[i];
                        } else {
                                invdiag[i]=0;
                        }
                }
                return invdiag;
        }

        Matrix<Size,Size,Precision>& get_evectors() {return my_evectors;}

        const Matrix<Size,Size,Precision>& get_evectors() const {return my_evectors;}


        Vector<Size, Precision>& get_evalues() {return my_evalues;}
        const Vector<Size, Precision>& get_evalues() const {return my_evalues;}

        bool is_posdef() const {
                for (int i = 0; i < my_evalues.size(); ++i) {
                        if (my_evalues[i] <= 0.0)
                                return false;
                }
                return true;
        }

        bool is_negdef() const {
                for (int i = 0; i < my_evalues.size(); ++i) {
                        if (my_evalues[i] >= 0.0)
                                return false;
                }
                return true;
        }

        Precision get_determinant () const {
                Precision det = 1.0;
                for (int i = 0; i < my_evalues.size(); ++i) {
                        det *= my_evalues[i];
                }
                return det;
        }

        Matrix<Size, Size, Precision> get_sqrtm () const {
                Vector<Size, Precision> diag_sqrt(my_evalues.size());
                // In the future, maybe throw an exception if an
                // eigenvalue is negative?
                for (int i = 0; i < my_evalues.size(); ++i) {
                        diag_sqrt[i] = std::sqrt(my_evalues[i]);
                }
                return my_evectors.T() * diagmult(diag_sqrt, my_evectors);
        }

        Matrix<Size, Size, Precision> get_isqrtm (const double condition=Internal::symeigen_condition_no) const {
                Vector<Size, Precision> diag_isqrt(my_evalues.size());

                // Because sqrt is a monotonic-preserving transformation,
                Precision max_diag = -my_evalues[0] > my_evalues[my_evalues.size()-1] ? (-std::sqrt(my_evalues[0])):(std::sqrt(my_evalues[my_evalues.size()-1]));
                // In the future, maybe throw an exception if an
                // eigenvalue is negative?
                for (int i = 0; i < my_evalues.size(); ++i) {
                        diag_isqrt[i] = std::sqrt(my_evalues[i]);
                        if(fabs(diag_isqrt[i]) * condition > max_diag) {
                                diag_isqrt[i] = 1/diag_isqrt[i];
                        } else {
                                diag_isqrt[i] = 0;
                        }
                }
                return my_evectors.T() * diagmult(diag_isqrt, my_evectors);
        }

private:
        // eigen vectors laid out row-wise so evectors[i] is the ith evector
        Matrix<Size,Size,Precision> my_evectors;

        Vector<Size, Precision> my_evalues;
};

}

#endif

---

libtoon
Author(s): Florian Weisshardt
autogenerated on Fri Jan 11 10:09:37 2013