Clean stuff up, no matrix inverse

theano/tensor/slinalg.py

 import logging
+import warnings
 
 logger = logging.getLogger(__name__)
 import numpy
@@ -370,17 +371,21 @@ class ExpmGrad(Op):
         return [shapes[0]]
 
     def perform(self, node, (A, gw), (out,)):
-        w, M = scipy.linalg.eig(A)
+        # Kalbfleisch and Lawless, J. Am. Stat. Assoc. 80 (1985) Equation 3.4
+        w, UL, UR = scipy.linalg.eig(A, left=True, right=True)
 
-        G = scipy.linalg.solve(M, gw).dot(M)
+        if numpy.linalg.norm(w - numpy.real(w)) > numpy.finfo(A.dtype).eps:
+            warnings.warn("ExpmGrad not correct for matrices with complex "
+                          "eigenvalues")
+
+        G = (UR.conj().T).dot(gw).dot(UL)
 
         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.multiply(V, G, V)
+        V = numpy.multiply(V, G, V)
 
-        Mi = scipy.linalg.inv(M)
-        out[0] = numpy.real(M.dot(V).dot(Mi))
+        out[0] = numpy.real(UL.dot(V).dot(UR.conj().T)).astype(A.dtype)
 
 
 def expm(A):

theano/tensor/tests/test_slinalg.py

@@ -224,31 +224,29 @@ def test_expm_grad_2():
         raise SkipTest("Scipy needed for the expm op.")
     rng = numpy.random.RandomState(utt.fetch_seed())
     A = rng.randn(3, 3).astype(config.floatX)
-    A = A
 
     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,1], x))
-    expm_f = function([x], expm(x)[0,1])
-
-    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)
-
+# 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)
+#