Initial commit

theano/tensor/slinalg.py

@@ -329,3 +329,58 @@ def kron(a, b):
                          o.shape[1] * o.shape[3]) +
                         tuple([o.shape[i] for i in range(4, o.ndim)]))
     return o
+
+
+class Expm(Op):
+    """Compute the matrix exponential of a square array
+    """
+
+    def make_node(self, A):
+        assert imported_scipy, (
+            "Scipy not available. Scipy is needed for the Expm op")
+
+        A = as_tensor_variable(A)
+        assert A.ndim == 2
+        expm = theano.tensor.matrix(dtype=A.dtype)
+        return Apply(self, [A,], [expm,])
+
+    def perform(self, node, (A,), (expm,)):
+        expm[0] = scipy.linalg.expm(A)
+
+    def grad(self, (A,), (g_out,)):
+        return [ExpmGrad()(A, g_out)]
+
+
+class ExpmGrad(Op):
+    """Gradient of the matrix exponential of a square array
+    """
+
+    def make_node(self, A, gw):
+        assert imported_scipy, (
+            "Scipy not available. Scipy is needed for the Expm op")
+        A = as_tensor_variable(A)
+        assert A.ndim == 2
+        print(gw.shape)
+        out = theano.tensor.matrix(dtype=A.dtype)
+        return Apply(self, [A, gw], [out,])
+
+    def perform(self, node, (A, gw), outputs):
+         w, M = scipy.linalg.eig(A)
+         outputs[0][0] = numpy.zeros_like(A)
+
+        # Mi = scipy.linalg.inv(M)
+        # G = Mi.dot(gw).dot(M)
+
+        # V = np.zeros_like(x)
+        # for i in range(V.shape[0]):
+        #     for j in range(V.shape[1]):
+        #         if i != j:
+        #             V[i,j] = G[i,j] * (exp(w[i]) - exp(w[j])) / (w[i] - w[j])
+        #         else:
+        #             V[i,j] = G[i,i] * exp(w[i])
+
+        # outputs[0] = M.dot(V).dot(Mi)
+
+
+def expm(A):
+    return Expm()(A)

theano/tensor/tests/test_slinalg.py

@@ -20,7 +20,8 @@ from theano.tensor.slinalg import ( Cholesky,
                                     solve,
                                     Eigvalsh,
                                     EigvalshGrad,
-                                    eigvalsh
+                                    eigvalsh,
+                                    expm
                                     )
 
 from nose.plugins.skip import SkipTest
@@ -189,3 +190,29 @@ class test_Solve(utt.InferShapeTester):
                                                dtype=config.floatX)],
                                 self.op_class,
                                 warn=False)
+
+
+def test_expm():
+    if not imported_scipy:
+        raise SkipTest("Scipy needed for the expm op.")
+    rng = numpy.random.RandomState(utt.fetch_seed())
+    A = rng.randn(5, 5).astype(config.floatX)
+
+    ref = scipy.linalg.expm(A)
+
+    x = tensor.matrix()
+    m = expm(x)
+    expm_f = function([x], m)
+
+    val = expm_f(A)
+    numpy.testing.assert_array_almost_equal(val, ref)
+
+
+def test_expm_grad():
+    if not imported_scipy:
+        raise SkipTest("Scipy needed for the expm op.")
+    rng = numpy.random.RandomState(utt.fetch_seed())
+    A = rng.randn(5, 5).astype(config.floatX)
+
+    eps = None
+    yield (lambda: utt.verify_grad(expm, [A], 3, rng, eps=eps))