Moving logdet to math

pymc3/distributions/multivariate.py

@@ -13,7 +13,7 @@
 
 import pymc3 as pm
 
-from pymc3.theanof import logdet
+from pymc3.math import logdet
 from . import transforms
 from .distribution import Continuous, Discrete, draw_values, generate_samples
 from ..model import Deterministic
@@ -124,10 +124,11 @@ def logp(self, value):
         delta = value - mu
         k = tau.shape[0]
 
+        result = k * tt.log(2 * np.pi)
         if self.gpu_compat:
-            result = k * tt.log(2 * np.pi) + tt.log(1. / det(tau))
+            result -= tt.log(det(tau))
         else:
-            result = k * tt.log(2 * np.pi) - logdet(tau)
+            result -= logdet(tau)
         result += (delta.dot(tau) * delta).sum(axis=delta.ndim - 1)
         return -1 / 2. * result
 

pymc3/math.py

@@ -10,6 +10,8 @@
     minimum, sgn, ceil, floor)
 from theano.tensor.nlinalg import det, matrix_inverse, extract_diag, matrix_dot, trace
 from theano.tensor.nnet import sigmoid
+from theano.gof import Op, Apply
+import numpy as np
 # pylint: enable=unused-import
 
 
@@ -25,3 +27,38 @@ def invlogit(x, eps=sys.float_info.epsilon):
 
 def logit(p):
     return tt.log(p / (1 - p))
+
+
+class LogDet(Op):
+    """Computes the logarithm of absolute determinant of a square
+    matrix M, log(abs(det(M))), on CPU. Avoids det(M) overflow/
+    underflow.
+
+    Note: Once PR #3959 (https://github.com/Theano/Theano/pull/3959/) by harpone is merged,
+    this must be removed. 
+    """
+    def make_node(self, x):
+        x = theano.tensor.as_tensor_variable(x)
+        o = theano.tensor.scalar(dtype=x.dtype)
+        return Apply(self, [x], [o])
+
+    def perform(self, node, inputs, outputs, params=None):
+        try:
+            (x,) = inputs
+            (z,) = outputs
+            s = np.linalg.svd(x, compute_uv=False)
+            log_det = np.sum(np.log(np.abs(s)))
+            z[0] = np.asarray(log_det, dtype=x.dtype)
+        except Exception:
+            print('Failed to compute logdet of {}.'.format(x))
+            raise
+
+    def grad(self, inputs, g_outputs):
+        [gz] = g_outputs
+        [x] = inputs
+        return [gz * matrix_inverse(x).T]
+
+    def __str__(self):
+        return "LogDet"
+
+logdet = LogDet()

pymc3/tests/test_math.py

@@ -0,0 +1,32 @@
+import numpy as np
+import theano
+from theano.tests import unittest_tools as utt
+from pymc3.math import LogDet, logdet
+from .helpers import SeededTest
+
+class TestLogDet(SeededTest):
+
+    def setUp(self):
+        utt.seed_rng()
+        self.op_class = LogDet
+        self.op = logdet
+
+    def validate(self, input_mat):
+        x = theano.tensor.matrix()
+        f = theano.function([x], self.op(x))
+        out = f(input_mat)
+        svd_diag = np.linalg.svd(input_mat, compute_uv=False)
+        numpy_out = np.sum(np.log(np.abs(svd_diag)))
+
+        # Compare the result computed to the expected value.
+        utt.assert_allclose(numpy_out, out)
+
+        # Test gradient:
+        utt.verify_grad(self.op, [input_mat])
+
+    def test_basic(self):
+        # Calls validate with different params
+        test_case_1 = np.random.randn(3, 3) / np.sqrt(3)
+        test_case_2 = np.random.randn(10, 10) / np.sqrt(10)
+        self.validate(test_case_1.astype(theano.config.floatX))
+        self.validate(test_case_2.astype(theano.config.floatX))

pymc3/tests/test_theanof.py

@@ -3,9 +3,6 @@
 import numpy as np
 import theano
 from ..theanof import DataGenerator, GeneratorOp, generator
-from theano.tests import unittest_tools as utt
-from pymc3.theanof import LogDet, logdet
-from .helpers import SeededTest
 
 def integers():
     i = 0
@@ -20,33 +17,6 @@ def integers_ndim(ndim):
         yield np.ones((2,) * ndim) * i
         i += 1
 
-class TestLogDet(SeededTest):
-
-    def setUp(self):
-        utt.seed_rng()
-        self.op_class = LogDet
-        self.op = logdet
-
-    def validate(self, input_mat):
-        x = theano.tensor.matrix()
-        f = theano.function([x], self.op(x))
-        out = f(input_mat)
-        svd_diag = np.linalg.svd(input_mat, compute_uv=False)
-        numpy_out = np.sum(np.log(np.abs(svd_diag)))
-
-        # Compare the result computed to the expected value.
-        utt.assert_allclose(numpy_out, out)
-
-        # Test gradient:
-        utt.verify_grad(self.op, [input_mat])
-
-    def test_basic(self):
-        # Calls validate with different params
-        test_case_1 = np.random.randn(3, 3) / np.sqrt(3)
-        test_case_2 = np.random.randn(10, 10) / np.sqrt(10)
-        self.validate(test_case_1.astype(theano.config.floatX))
-        self.validate(test_case_2.astype(theano.config.floatX))
-
 
 class TestGenerator(unittest.TestCase):
     def test_basic(self):

pymc3/theanof.py

@@ -4,17 +4,16 @@
 from .vartypes import typefilter, continuous_types
 from theano import theano, scalar, tensor as tt
 from theano.gof.graph import inputs
-from theano.gof import Op, Apply
+from theano.gof import Op
 from theano.configparser import change_flags
-from theano.tensor.nlinalg import matrix_inverse
 from .memoize import memoize
 from .blocking import ArrayOrdering
 from .data import DataGenerator
 
 __all__ = ['gradient', 'hessian', 'hessian_diag', 'inputvars',
            'cont_inputs', 'floatX', 'jacobian',
            'CallableTensor', 'join_nonshared_inputs',
-           'make_shared_replacements', 'generator', 'LogDet', 'logdet']
+           'make_shared_replacements', 'generator']
 
 
 def inputvars(a):
@@ -61,40 +60,6 @@ def floatX(X):
 Theano derivative functions
 """
 
-class LogDet(Op):
-    """Computes the logarithm of absolute determinant of a square
-    matrix M, log(abs(det(M))), on CPU. Avoids det(M) overflow/
-    underflow.
-
-    Note: Once PR #3959 (https://github.com/Theano/Theano/pull/3959/) by harpone is merged,
-    this must be removed. 
-    """
-    def make_node(self, x):
-        x = theano.tensor.as_tensor_variable(x)
-        o = theano.tensor.scalar(dtype=x.dtype)
-        return Apply(self, [x], [o])
-
-    def perform(self, node, inputs, outputs):
-        try:
-            (x,) = inputs
-            (z,) = outputs
-            s = np.linalg.svd(x, compute_uv=False)
-            log_det = np.sum(np.log(np.abs(s)))
-            z[0] = np.asarray(log_det, dtype=x.dtype)
-        except Exception:
-            print('Failed to compute logdet of {}.'.format(x))
-            raise
-
-    def grad(self, inputs, g_outputs):
-        [gz] = g_outputs
-        [x] = inputs
-        return [gz * matrix_inverse(x).T]
-
-    def __str__(self):
-        return "LogDet"
-
-logdet = LogDet()
-
 def gradient1(f, v):
     """flat gradient of f wrt v"""
     return tt.flatten(tt.grad(f, v, disconnected_inputs='warn'))