Working with eigendecomposition tools in TMB

[James-Thorson]

Hi all,

I've been looking through the functions distributed with TMB, but can't seem to work out how to use the options for computing eigendecomposition. Is there any example code I could look at?

I'm specifically looking to generate an unstructured p by p matrix B that is asymmetric but negative semi-definite with rank r, where r<=p. The kludge I envision is generating an asymmetric p by p matrix with rank r, and then fixing positive eigenvalues to zero via the posfun function. I welcome ideas on this applied problem in addition to help with the eigendecomposition tools in TMB.

cheers, jim

[skaug]

I do not have a complete answer to this, but you should seek for a subroutine inside the Eigen library, which is loaded by TMB. Note that library Eigen is doing much more than eigen-decompositon (in fact it is the basis for all the numerical linear algebra in TMB).

Maybe some of these routines will do the job? http://eigen.tuxfamily.org/dox/group__Eigenvalues__Module.html

If somebody finds a solution it would be nice if it could be posted to the code snippets: https://github.com/kaskr/adcomp/wiki/Code--snippets

[James-Thorson]

After taking another kick at this, I see that the following code works for extracting eigenvalues for dense, nonsymetric matrix Q_yy (it took a while because I hadn't seen much use of complex numbers in TMB, and e.g. the REPORT macro doesn't appear to work with complex numbers, etc).

  vector< std::complex< Type > > eigenvalues_Q_yy = Q_yy.eigenvalues();
  vector< Type > real_eigenvalues_Q_yy = eigenvalues_Q_yy.real();
  vector< Type > imag_eigenvalues_Q_yy = eigenvalues_Q_yy.imag();
  REPORT( real_eigenvalues_Q_yy );
  REPORT( imag_eigenvalues_Q_yy );

I still haven't figured out how the EigenSolver constructor works, or how to efficiently get eigenvectors. But I think the eigenvalues should be enough for my current purposes. If anyone figures out how to get eigenvectors, please chime in.

cheers, jim

[James-Thorson]

Hi all,

I'm still having trouble figuring out how to extract eigenvectors within TMB, and still have a problem where having access to eigenvectors would provide a more general solution than my current toolbox. The following example (from Eigen Doxygen documentation) works:

    MatrixXd A = MatrixXd::Random(6,6);
    EigenSolver<MatrixXd> es(A);
    es.eigenvalues();
    es.eigenvectors();

but the modified code for differentiable Type doesn't work:

    matrix<Type> A;
    A.setZero();
    EigenSolver< matrix<Type> > es(A);
    es.eigenvalues();
    es.eigenvectors();

Does anyone have advice for a general solution to calculating eigenvectors in TMB?

[calbertsen]

The documentation for EigenSolver states that the template parameter "[...] is expected to be an instantiation of the Matrix class template. Currently, only real matrices are supported."

It is not a general solution but the SelfAdjointEigenSolver does not require a real matrix and seems to work.

#include <TMB.hpp>
#include <Eigen/Eigenvalues>

template<class Type>
Type objective_function<Type>::operator() ()
{
  DATA_INTEGER(i);
  PARAMETER(mu);
  using namespace Eigen;
  Matrix<Type,Dynamic,Dynamic> A(i,i);
  A.setZero();
  SelfAdjointEigenSolver<Matrix<Type,Dynamic,Dynamic> > es(A);
  Rcout << "Eigen values: " << es.eigenvalues() << "\n\n";
  Rcout << "Eigen vectors: " << es.eigenvectors() << "\n\n";
  return mu;
}

[bhattsamir]

Hey all, did anyone ever get a definitive answer to this?

Using

SelfAdjointEigenSolver<Matrix<Type,Dynamic,Dynamic> > es(Kdi); 
matrix<Type> V = es.eigenvectors();
vector<Type> Ev = es.eigenvalues();

Is the only syntax that works for me. The problem is, when comparing the results to R's eigen function, the values match perfectly, but the vectors have some strange sign switching...

It is entirely possible (and likely) that I have messed something up.

Does anyone have a syntax that they have tested?

something along the lines of what james was doing

    matrix<Type> A;
    A.setZero();
    EigenSolver< matrix<Type> > es(A);
    es.eigenvalues();
    es.eigenvectors();

Thanks for the help, and thanks for TMB! fastest software out there!

[kaskr]

Eigenvectors are not unique: If x is an eigenvector then so is -x. You should be able to match each of es.eigenvectors() with one of R's eigenvectors multiplied by either 1 or -1.