Implemented linalg.eigh.

theano/sandbox/linalg/__init__.py

 
 from ops import (cholesky, matrix_inverse, solve,
         diag, extract_diag, alloc_diag,
-        det, psd, eig,
+        det, psd, eig, eigh,
         trace, spectral_radius_bound)

theano/sandbox/linalg/ops.py

@@ -882,8 +882,8 @@ class Eig(Op):
 
     """
 
-    def __init__(self):
-        pass
+    def __init__(self, numop):
+        self._numop = numop
 
     def props(self):
         """Function exposing different properties of each instance of the
@@ -907,9 +907,10 @@ class Eig(Op):
 
     def perform(self, node, (x,), (w, v)):
         try:
-            w[0], v[0] = [z.astype(x.dtype) for z in numpy.linalg.eig(x)]
+            w[0], v[0] = [z.astype(x.dtype) for z in self._numop(x)]
         except numpy.linalg.LinAlgError:
-            logger.debug('Failed to find eig of %s' % str(node.inputs[0]))
+            logger.debug('Failed to find %s of %s' % (node.inputs[0], 
+                                                      self._numop.__name__))
             raise
 
     def infer_shape(self, node, shapes):
@@ -917,7 +918,7 @@ class Eig(Op):
         return [(n,), (n,n)]
 
     def __str__(self):
-        return "Eig"
+        return self._numop.__name__.capitalize()
 
     def grad(self, inputs, g_outputs):
         r"""The gradient function should return
@@ -943,7 +944,8 @@ class Eig(Op):
         gw, gv = g_outputs
         return [EigGrad()(x, w, v, gw, gv)]
 
-eig = Eig()
+eig = Eig(numpy.linalg.eig)
+eigh = Eig(numpy.linalg.eigh)
 
 class EigGrad(Op):
     """Gradient of an eigensystem.
@@ -968,6 +970,16 @@ class EigGrad(Op):
 
     def perform(self, node, inputs, outputs):
         """
+        Implements the "reverse-mode" gradient for the eigensystem of
+        a square matrix.
+
+        Let
+
+            .. math:: w, v = \mbox{eig}(x).
+
+        By definition of the eigensystem,
+
+            .. math:: \sum_j x_{ij}\,v_{jn} = w_n\,v_{in}.
         """
         x, w, v, gw, gv = inputs
         N = x.shape[0]

theano/sandbox/linalg/tests/test_linalg.py

@@ -29,7 +29,7 @@ from theano.sandbox.linalg.ops import (cholesky,
                                        imported_scipy,
                                        Eig,
                                        )
-from theano.sandbox.linalg import eig
+from theano.sandbox.linalg import eig, eigh
 from nose.plugins.skip import SkipTest
 
 
@@ -471,10 +471,10 @@ class test_Solve(utt.InferShapeTester):
                                 self.op_class)
 
 class test_Eig(utt.InferShapeTester):
+    op_class = Eig
+    op = eig
     def setUp(self):
         super(test_Eig, self).setUp()
-        self.op_class = Eig
-        self.op = Eig()
         self.rng = numpy.random.RandomState(utt.fetch_seed())
         self.A = theano.tensor.matrix()
         X = numpy.asarray(self.rng.rand(5, 5),
@@ -494,16 +494,15 @@ class test_Eig(utt.InferShapeTester):
         A = theano.tensor.matrix()
         self.assertEquals([e.eval({A: [[1]]}) for e in self.op(A)],
                           [[1.0], [[1.0]]])
+        x = [[0, 1], [1, 0]] 
+        w, v = [e.eval({A: x}) for e in self.op(A)]
+        assert_array_almost_equal(numpy.dot(x,v), w * v)
 
-        w, v = [e.eval({A: [[0, 1], [1, 0]]}) 
-                for e in self.op(A)]
-        assert_array_almost_equal(w, [1, -1])
-        x = math.sqrt(2)/2
-        assert_array_almost_equal(v, [[x, -x], [x, x]])
-
+class test_Eigh(test_Eig):
+    op = eigh
     def test_grad(self):
         S = self.S
         def fun(x):
-            w, v = eig(x)
+            w, v = self.op(x)
             return w.sum() + v.sum()*0
         utt.verify_grad(fun, [S], rng=self.rng)