Hermetian args

theano/tensor/slinalg.py

@@ -356,7 +356,8 @@ class Expm(Op):
 
 
 class ExpmGrad(Op):
-    """Gradient of the matrix exponential of a square array
+    """Gradient of the matrix exponential of a square array. Currently
+    implemented only for Hermetian arguments.
     """
 
     def make_node(self, A, gw):
@@ -370,22 +371,21 @@ class ExpmGrad(Op):
     def infer_shape(self, node, shapes):
         return [shapes[0]]
 
-    def perform(self, node, (A, gw), (out,)):
-        # Kalbfleisch and Lawless, J. Am. Stat. Assoc. 80 (1985) Equation 3.4
-        w, UL, UR = scipy.linalg.eig(A, left=True, right=True)
-
-        if numpy.linalg.norm(w - numpy.real(w)) > numpy.finfo(A.dtype).eps:
-            warnings.warn("ExpmGrad not correct for matrices with complex "
-                          "eigenvalues")
+    def perform(self, node, (A, gA), (out,)):
+        if not numpy.allclose(A, A.T):
+            raise NotImplementedError(
+                "ExpmGrad is only implemented for symmetric matrices")
 
-        G = (UR.conj().T).dot(gw).dot(UL)
+        # Kalbfleisch and Lawless, J. Am. Stat. Assoc. 80 (1985) Equation 3.4
+        w, U = scipy.linalg.eigh(A)
 
+        G = (U.T).dot(gA).dot(U)
         exp_w = numpy.exp(w)
         V = numpy.subtract.outer(exp_w, exp_w) / numpy.subtract.outer(w, w)
-        V[numpy.diag_indices_from(V)] = exp_w
+        numpy.fill_diagonal(V, exp_w)
         V = numpy.multiply(V, G, V)
 
-        out[0] = numpy.real(UL.dot(V).dot(UR.conj().T)).astype(A.dtype)
+        out[0] = (U.dot(V).dot(U.T)).astype(A.dtype)
 
 
 expm = Expm()

theano/tensor/tests/test_slinalg.py

@@ -217,36 +217,3 @@ def test_expm_grad_1():
     A = A + A.T
 
     tensor.verify_grad(expm, [A,], rng=rng)
-
-def test_expm_grad_2():
-    # with non-symmetric matrix
-    if not imported_scipy:
-        raise SkipTest("Scipy needed for the expm op.")
-    rng = numpy.random.RandomState(utt.fetch_seed())
-    A = rng.randn(3, 3).astype(config.floatX)
-
-    tensor.verify_grad(expm, [A,], rng=rng)
-
-# def test_expm_grad_3():
-#     if not imported_scipy:
-#         raise SkipTest("Scipy needed for the expm op.")
-#     from theano.gradient import grad
-#     rng = numpy.random.RandomState(utt.fetch_seed())
-#     A = rng.randn(3, 3).astype(config.floatX)
-#
-#     h = 1e-7
-#     def e(i,j):
-#         v = numpy.zeros((3, 3))
-#         v[i, j] = 1
-#         return v
-#
-#     x = tensor.matrix()
-#     grad_expm_f = function([x], grad(expm(x)[0,0], x))
-#     expm_f = function([x], expm(x)[0,0])
-#
-#     g = lambda i, j:  (expm_f(A + h*e(i,j)) - expm_f(A)) / h
-#     numgrad = numpy.array([[g(i,j) for i in range(3)] for j in range(3)])
-#
-#     print(grad_expm_f(A))
-#     print(numgrad)
-#