computing eigenvalues with Theano

I want to use theano.tensor.nlinalg.eig to compute the eigenvalues of a non-symmetric square matrix.

The question is: How do I get the complex values of the eigenvalues/eigenvectors?

It seems that theano.tensor.nlinalg.eig only returns the real part of them (see code below).

Thanks all.

code:

import numpy as np
import theano
rng = np.random
T = theano.tensor

# Create an asymmetric (random) square matrix
N = 3
asymm_Matrix = rng.randn(N,N)

# Compute eigenv* with numpy
np_eigenvalues, np_eigenvectors = np.linalg.eig(asymm_Matrix)

# Define a theano node that computes eigenv*
symMatrix = T.dmatrix("symMatrix")
symEigenvalues, eigenvectors = T.nlinalg.eig(symMatrix)
get_Eigen = theano.function([symMatrix], [symEigenvalues, eigenvectors] )

# Compute eigenv* with theano
theano_eigenvalues, theano_eigenvectors = get_Eigen(asymm_Matrix)

print("---- asymm_Matrix :")
print(asymm_Matrix)
print()
print("---- np.linalg.eig")
print("eigenvalues :")
print(np_eigenvalues)
print("eigenvectors :")
print(np_eigenvectors)
print()
print("---- T.nlinalg.eig")
print("eigenvalues :")
print(theano_eigenvalues)
print("eigenvectors :")
print(theano_eigenvectors)

output:

 ---- asymm_Matrix :
 [[-0.163  -0.2099  1.1227]
  [-1.132  -0.9667 -1.0436]
  [-0.0974 -0.4395 -0.3839]]

 ---- np.linalg.eig
 eigenvalues :
 [-1.0136+0.1346j -1.0136-0.1346j  0.5136+0.j    ]
 eigenvectors :
 [[-0.3948-0.1604j -0.3948+0.1604j  0.6638+0.j    ]
  [ 0.7736+0.j      0.7736-0.j     -0.6977+0.j    ]
  [ 0.4630+0.0742j  0.4630-0.0742j  0.2696+0.j    ]]

 ---- T.nlinalg.eig
 eigenvalues :
 [-1.0136 -1.0136  0.5136]
 eigenvectors :
 [[-0.3948 -0.3948  0.6638]
  [ 0.7736  0.7736 -0.6977]
  [ 0.463   0.463   0.2696]]

Solution

I did not find a clear documentation of how T.nlinalg.eig works but it seems that its output has always the same type as the inputs (see line 321 in nlinalg.eig code https://github.com/Theano/Theano/blob/master/theano/tensor/nlinalg.py).

So, to get a complex output, I simply defined your input symMatrix as complex: symMatrix = T.cmatrix("symMatrix").

As symMatrix is now complex, when you call get_Eigen, the argument asymm_Matrix also has to be complex, so I converted it accordingly into a new variable: c_asymm_Matrix = asymm_Matrix.astype('complex64').

The following code returns the imaginary parts as well:

import numpy as np
import theano
rng = np.random
T = theano.tensor

# Create an asymmetric (random) square matrix
N = 3
asymm_Matrix = rng.randn(N,N)

# Compute eigenv* with numpy
np_eigenvalues, np_eigenvectors = np.linalg.eig(asymm_Matrix)

# Define a theano node that computes eigenv*
# symMatrix defined as complex so that nlinalg.eig outputs are complex as well
symMatrix = T.cmatrix("symMatrix") 
symEigenvalues, eigenvectors = T.nlinalg.eig(symMatrix)
get_Eigen = theano.function([symMatrix], [symEigenvalues, eigenvectors] ) 

# Compute eigenv* with theano
# Our input for get_Eigen has to be complex too now
c_asymm_Matrix = asymm_Matrix.astype('complex64') 

theano_eigenvalues, theano_eigenvectors = get_Eigen(c_asymm_Matrix) #this one gives error because it expects float64 input and tries to downcast

print("---- asymm_Matrix :")
print(asymm_Matrix)
print()
print("---- np.linalg.eig")
print("eigenvalues :")
print(np_eigenvalues)
print("eigenvectors :")
print(np_eigenvectors)
print()
print("---- T.nlinalg.eig")
print("eigenvalues :")
print(theano_eigenvalues)
print("eigenvectors :")
print(theano_eigenvectors)

Output:

---- asymm_Matrix :
[[ 1.16270655  0.7266984  -0.87479655]
 [-3.04065245 -1.09469397 -1.57315266]
 [ 0.4850981   0.56072627  1.16248949]]
()
---- np.linalg.eig
eigenvalues :
[ 0.03696641+1.49515503j  0.03696641-1.49515503j  1.15656924+0.j        ]
eigenvectors :
[[ 0.22084674+0.37455923j  0.22084674-0.37455923j -0.66915708+0.j        ]
 [-0.87172957+0.j         -0.87172957-0.j          0.57392013+0.j        ]
 [ 0.20022446+0.10454574j  0.20022446-0.10454574j  0.47206407+0.j        ]]
()
---- T.nlinalg.eig
eigenvalues :
[ 0.03696636 +1.49515510e+00j  0.03696636 -1.49515510e+00j
  1.15656924 +1.52098596e-17j]
eigenvectors :
[[-0.22084673 -3.74559224e-01j -0.22084673 +3.74559224e-01j
   0.66915709 +0.00000000e+00j]
 [ 0.87172955 +0.00000000e+00j  0.87172955 +0.00000000e+00j
  -0.57392007 +5.55111512e-17j]
 [-0.20022446 -1.04545735e-01j -0.20022446 +1.04545735e-01j
  -0.47206411 +3.88578059e-16j]]

The output is slightly different but it's just because some very small imaginary parts which are displayed in the Theano output are zero in the numpy output.