add numeric.eig, update function.Eig

nutils/function.py

@@ -1790,7 +1790,7 @@ def __init__( self, func, symmetric=False, sort=False ):
     self.symmetric = symmetric
     self.func = func
     self.shape = func.shape
-    self.eig = numpy.linalg.eigh if symmetric else numpy.linalg.eig 
+    self.eig = numpy.linalg.eigh if symmetric else numeric.eig
 
   def evalf( self, arr ):
     assert arr.ndim == len(self.shape)+1

nutils/numeric.py

@@ -295,6 +295,34 @@ def diagonalize( arg ):
   diag[:] = arg
   return diagonalized
 
+def eig( A ):
+  '''If A has repeated eigenvalues, numpy.linalg.eig sometimes fails to produce
+  the complete eigenbasis. This function aims to fix that by identifying the
+  problem and completing the basis where necessary.'''
+
+  L, V = numpy.linalg.eig( A )
+
+  # check repeated eigenvalues
+  for index in numpy.ndindex( A.shape[:-2] ):
+    unique, inverse = numpy.unique( L[index], return_inverse=True )
+    if len(unique) < len(inverse): # have repeated eigenvalues
+      repeated, = numpy.where( numpy.bincount(inverse) > 1 )
+      vectors = V[index].T
+      for i in repeated: # indices pointing into unique corresponding to repeated eigenvalues
+        where, = numpy.where( inverse == i ) # corresponding eigenvectors
+        for j, n in enumerate(where):
+          W = vectors[where[:j]]
+          vectors[n] -= numpy.dot( numpy.dot( W, vectors[n] ), W ) # gram schmidt orthonormalization
+          scale = numpy.linalg.norm(vectors[n])
+          if scale < 1e-8: # vectors are near linearly dependent
+            u, s, vh = numpy.linalg.svd( A - unique[i] * numpy.eye(len(A)) )
+            nnz = numpy.argsort( abs(s) )[:len(where)]
+            vectors[where] = vh[nnz].conj()
+            break
+          vectors[n] /= scale
+
+  return L, V
+
 def isbool( a ):
   return isboolarray( a ) and a.ndim == 0 or numpy.issubdtype( type(a), numpy.bool )