How to get Eigenvalues and Eigenvectors of Nonsymmetric Matrices?

[apourz...@evidation.com]

Hi,

I have a nonsymmetric matrix Q, which is actually what's called a generator
matrix in continuous-time markov chain theory (see:
https://en.wikipedia.org/wiki/Continuous-time_Markov_chain). I need to get
the eigenvectors and eigenvalues of Q so that I can find probabilities of
state transitions computationally efficiently.

It seems like the Stan users removed the eigenvector() and eigenvalue()
functions. Why is this? Is there still a way to get the eigenvectors and
eigenvalues of a nonsymmetric matrix?

Best Regards,
Arya

[Ben Goodrich]

On Tuesday, June 23, 2015 at 8:26:43 PM UTC-4, apourz...@evidation.com wrote:

It seems like the Stan users removed the eigenvector() and eigenvalue()
functions. Why is this? 

In general, we can't guarantee that the eigenvalues are real and can't autodiff complex numbers.
 

Is there still a way to get the eigenvectors and
eigenvalues of a nonsymmetric matrix?

I have a branch that partially implements what Eigen calls a pseudo-eigendecomposition

http://eigen.tuxfamily.org/dox/classEigen_1_1EigenSolver.html#a4140972e2b45343d1ef1793c2824159c

that could be used in a wider variety of circumstances, including those where a nonsymmetric matrix happened to have all real eigenvalues.

Ben

[Arya Pourzanjani]

I need a little help understanding. My log posterior probability includes the sum of all the elements in a matrix P where P=U exp(tD) U^-1 and Q=UDU^-1. Q is my actual matrix of parameters I'm trying to learn and UDU^-1 is the eigendecomposition of Q. The unnormalized log posterior density seems easy to compute, but the derivative of P with respect to Q is only possible to symbolically differentiate when the eigenvalues of Q are all real? 

Is there any way I can get around this to implement my model in Stan to get the benefits of the NUTS? Is it possible to use HMC at all or do I have to use a sampler that doesn't rely on derivative like Metropolis Hastings?

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/QJe1TNioiyg/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Ben Goodrich]

On Wednesday, June 24, 2015 at 1:59:18 PM UTC-4, Arya Pourzanjani wrote:

I need a little help understanding. My log posterior probability includes the sum of all the elements in a matrix P where P=U exp(tD) U^-1 and Q=UDU^-1. Q is my actual matrix of parameters I'm trying to learn and UDU^-1 is the eigendecomposition of Q. The unnormalized log posterior density seems easy to compute, but the derivative of P with respect to Q is only possible to symbolically differentiate when the eigenvalues of Q are all real? 

Is there any way I can get around this to implement my model in Stan to get the benefits of the NUTS? Is it possible to use HMC at all or do I have to use a sampler that doesn't rely on derivative like Metropolis Hastings?

A slight limitation of Stan is that it cannot do autodifferentiation of complex numbers. So, that means either

• The log-posterior has to be composed entirely of operations on real numbers (so that autodifferentiation can be used)
• Any operation involving complex numbers would require a C++ function that demotes Stan's custom scalars to plain double precision numbers, calls a subfunction, calculates the necessary derivatives analytically in double precision, and puts the realizations of these derivatives into the right places (which is tedious)

There are some constructions such that Q is guaranteed to have real eigenvalues, in which case autodifferentiation would work fine if Stan implemented the eigendecomposition of an asymmetric matrix. But Stan doesn't implement the eigendecomposition of an asymmetric matrix because it generally can't verify from its elements that an asymmetric matrix has real eigenvalues only (I suppose we could implement it for some verifiable special cases like a triangular matrix or a right stochastic matrix).

If you know that your Q only has real eigenvalues by construction, then the easiest thing to do would be to implement the pseudo-eigendecomposition locally. Put something like

#ifndef STAN__MATH__MATRIX__PSEUDOEIGENDECOMPOSITION_HPP
#define STAN__MATH__MATRIX__PSEUDOEIGENDECOMPOSITION_HPP

#include <stan/math/prim/mat/fun/Eigen.hpp>
#include <Eigen/Eigenvalues>
#include <stan/math/prim/scal/err/check_square.hpp>
#include <stan/math/prim/scal/err/check_nonzero_size.hpp>
#include <stan/math/prim/scal/err/check_finite.hpp>

namespace stan {
  namespace math {

    template <typename T>
    std::vector<Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> >
    pseudoeigendecomposition(const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m) {
      typedef Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> matrix_t;
      stan::math::check_nonzero_size("pseudoeigendecomposition", "m", m);
      stan::math::check_square("pseudoeigendecomposition", "m", m);
      stan::math::check_finite("pseudoeigendecomposition", "m", m);
      Eigen::EigenSolver<matrix_t> es(m);
      std::vector<matrix_t> out(2);
      out[0] = es.pseudoEigenvectors();
      out[1] = es.pseudoEigenvalueMatrix();
      return out;
    }
  }
}
#endif

into a file called lib/stan_math_2.7.0/stan/math/prim/mat/fun/pseudoeigendecomposition.hpp once Stan 2.7.0 comes out (doing this for Stan 2.6.x would involve a different path). And then change src/stan/lang/function_signatures.h to include the line

add("pseudoeigendecomposition",MATRIX_T,expr_type(MATRIX_T,1));

That should work (modulo typos) and the pseudoEigenvalueMatrix() will be diagonal since all the eigenvalues are real.

Ben

[Ben Goodrich]

On Wednesday, June 24, 2015 at 3:13:00 PM UTC-4, Ben Goodrich wrote:

into a file called lib/stan_math_2.7.0/stan/math/prim/mat/fun/pseudoeigendecomposition.hpp once Stan 2.7.0 comes out (doing this for Stan 2.6.x would involve a different path). And then change src/stan/lang/function_signatures.h to include the line

add("pseudoeigendecomposition",MATRIX_T,expr_type(MATRIX_T,1)); 

Also, add

#include <stan/math/prim/mat/fun/pseudoeigendecomposition.hpp>

to lib/stan_math_2.7.0/stan/math.hpp

Ben

[Ben Goodrich]

On Wed, Jun 24, 2015 at 4:55 PM, Arya Pourzanjani <apourz...@evidation.com> wrote:

Unfortunately, I can't guarantee my Q has real eigenvalues.

In that case, you have a deeper problem that I described at

https://groups.google.com/d/msg/stan-users/KUKqOnIzMGk/u59Y-OukYpEJ

where you are essentially putting positive prior density on proposals for the parameters that are not admissible if Q having real eigenvalues is a constraint. (If complex eigenvalues are admissible, then you have the more manageable problem of handling the real and complex parts separately).

 

Is the problem specifically that Stan's auto-differentiation module doesn't support derivatives of complex number

Yes

 

or that auto-differentiation of complex numbers if not possible?

I wouldn't be surprised if it is possible, but I have not seen any software that does so.

 

For example would PyMC3 which uses Theano for auto-differentiation be able to accomplish what I want?

I don't think Theano knows derivatives of eigenvalues and eigenvectors, especially for an asymmetric matrix.

 

I'd much rather prefer to use Stan, but I guess my only option would be to do the second bullet point you mentioned. How can I get an idea for how difficult that would that be? If it means I get to use Stan with all its optimizations and easy interface it could be worth it.

If having eigenvalues with non-zero imaginary parts is admissible, then what I wrote in my previous post is not quite sufficient because Eigen utilizes a few unnecessary operations on complex numbers when it does a pseudo-eigendecomposition. I do have some local C++ code that reimplemented that in Stan with only real operations, but it would probably be better if I proposed that in a patch to the Eigen developers. But then the question remains if you had V and D such that QV = VD where D has (some) 2x2 blocks along its diagonal, could you write your log-posterior handling the real and imaginary parts of the eigenvalues separately?

Ben

[Ben Goodrich]

On Wed, Jun 24, 2015 at 5:58 PM, Arya Pourzanjani <apourz...@evidation.com> wrote:

Q having complex eigenvalues is indeed admissible. In fact, in my toy example, the TRUE value of Q does indeed have complex eigenvalues.

> Q

     [,1] [,2] [,3] [,4] [,5]

[1,]   -6    2    2    1    1

[2,]    1   -4    0    1    2

[3,]    1    0   -4    2    1

[4,]    2    1    0   -3    0

[5,]    1    1    1    1   -4

> eigen(Q)

$values

[1] -6.453264e+00+1.63114e-01i -6.453264e+00-1.63114e-01i -4.046736e+00+7.02911e-01i -4.046736e+00-7.02911e-01i -5.907916e-16+0.00000e+00i

But then the question remains if you had V and D such that QV = VD where D has (some) 2x2 blocks along its diagonal, could you write your log-posterior handling the real and imaginary parts of the eigenvalues separately?

Not sure what you mean here. Could you elaborate? Maybe it's because I don't understand what the "pseudo-eigendecomposition" is. I don't think 2x2 along the diagonal would work, because I need to compute a matrix exponential and that's only computationally feasible when the matrix is diagonal.  

Eigen doesn't seem to provide a reference for it other than it being real matrices V and D where D may have 2x2 blocks along its diagonal such that QV = VD and the eigenvalues of Q are the eigenvalues of D. It is a function of the real Schur decomposition.

But a pseudo-eigendecomposition (if it were modified to only use real operations internally) may be adequate in your case because the matrix exponential of a block diagonal matrix is a block diagonal matrix of matrix exponentials and the matrix exponential of a 2x2 matrix is known.

Ben

[Ben Goodrich]

On Wed, Jun 24, 2015 at 6:54 PM, Arya Pourzanjani <apourz...@evidation.com> wrote:

That's good news! Does this just mean I can download your branch that implements pseudo-eigendecomposition and go to work? I would have to create a function in my Stan model file that computes the matrix exponential of the block diagonal matrix, but that doesn't sound too difficult.

Not quite. We really need an unwritten patch to Eigen first because my local implementation was for a generalized eigenvalue problem rather than a plain eigenvalue problem.

Ben

[Sebastian Weber]

Hi!

It would be really great if a function could be included in Stan which calculates the matrix exponential for cases when the eigenvalues are real. I understand the difficulty that you cannot validate the user-input, but how about giving huge Warning messages during model startup or mode compilation?

I am interested in compartmental PK systems and these are guaranteed to have real eigenvalues for all the practical cases I would be interested. If I could calculated the repeated exp(At) function - this would be great.

Best,

Sebastian

... using a patched Stan version for my projects is not an option...

[Ben Goodrich]

On Thursday, June 25, 2015 at 2:34:04 PM UTC-4, Sebastian Weber wrote:

It would be really great if a function could be included in Stan which calculates the matrix exponential for cases when the eigenvalues are real.

I think a pseudoeigendecomposition will be fine for many purposes (whether the eigenvalues are real or not), presuming we can get Eigen to accept the patch.
 

I understand the difficulty that you cannot validate the user-input, but how about giving huge Warning messages during model startup or mode compilation?

I am not a fan of how much mandatory input validation Stan does, but in this case, Stan would segfault even if the imaginary parts were zero. We have to change Eigen or reimplement it to not call Eigen.

Ben

 

[Arya Pourzanjani]

If having eigenvalues with non-zero imaginary parts is admissible, then what I wrote in my previous post is not quite sufficient because Eigen utilizes a few unnecessary operations on complex numbers when it does a pseudo-eigendecomposition.

Are you saying that Eigen is too inefficient at pseudo-eigendecomposition involving complex numbers or that Eigen simply can't do it for some reason?

 

I do have some local C++ code that reimplemented that in Stan with only real operations, but it would probably be better if I proposed that in a patch to the Eigen developers. But then the question remains if you had V and D such that QV = VD where D has (some) 2x2 blocks along its diagonal, could you write your log-posterior handling the real and imaginary parts of the eigenvalues separately?

What do you mean by handling the real and imaginary parts of the log-posterior separately?  

 

Not quite. We really need an unwritten patch to Eigen first because my local implementation was for a generalized eigenvalue problem rather than a plain eigenvalue problem.

 I'm still not sure what the Eigen patch does exactly.

[Ben Goodrich]

On Thursday, June 25, 2015 at 3:35:58 PM UTC-4, Arya Pourzanjani wrote:

If having eigenvalues with non-zero imaginary parts is admissible, then what I wrote in my previous post is not quite sufficient because Eigen utilizes a few unnecessary operations on complex numbers when it does a pseudo-eigendecomposition.

Are you saying that Eigen is too inefficient at pseudo-eigendecomposition involving complex numbers or that Eigen simply can't do it for some reason?

It uses complex numbers internally when it doesn't need to, even though the input matrix is real and the output matrices are real. So, we need to avoid that to make it not segfault with Stan.
 

I do have some local C++ code that reimplemented that in Stan with only real operations, but it would probably be better if I proposed that in a patch to the Eigen developers. But then the question remains if you had V and D such that QV = VD where D has (some) 2x2 blocks along its diagonal, could you write your log-posterior handling the real and imaginary parts of the eigenvalues separately?

What do you mean by handling the real and imaginary parts of the log-posterior separately?  

You don't really need the eigenvalues to calculate exp(D) but if you did, then for any complex number

c = a + i * b

you only get the real (a) and imaginary (b) parts. So, you can manually do some things like abs(c) = sqrt(a^2 + b^2) by handling the real and imaginary parts separately.

Ben

[Christopher Fisher]

Hi all-

I need to use pseudoeigendecomposition function but I am having trouble getting it work. I modified the paths to reflect the change in folder structure, but Stan still cannot find the function. Would someone be able to provide updated instructions?

Thank you,

Chris 

On Wednesday, June 24, 2015 at 3:13:00 PM UTC-4, Ben Goodrich wrote:

[Ben Goodrich]

On Wednesday, March 15, 2017 at 2:52:22 PM UTC-4, Christopher Fisher wrote:

I need to use pseudoeigendecomposition function but I am having trouble getting it work. I modified the paths to reflect the change in folder structure, but Stan still cannot find the function. Would someone be able to provide updated instructions?

There are instructions for how to add a function in the wiki, but in this particular case, thos instructions would not get you very far. See

http://eigen.tuxfamily.org/bz/show_bug.cgi?id=1029

for why. I have a reimplementation implementation on a branch.

Ben

[Bob Carpenter]

A branc on stan-dev/math? Is there anything holding up
merging it?

- Bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.

[Ben Goodrich]

On Sunday, March 19, 2017 at 7:28:17 PM UTC-4, Bob Carpenter wrote:

A branc on stan-dev/math?  Is there anything holding up
merging it?

Pushing, testing, dealing with fwd mode, free time, etc.