From 1c922f0a3176c708831a57eb0f5fadea3a82f789 Mon Sep 17 00:00:00 2001 From: Bhargav Srinivasa Date: Sun, 19 Feb 2017 12:53:43 +0530 Subject: [PATCH 1/3] Moved logdet to math --- pymc3/distributions/multivariate.py | 2 +- pymc3/math.py | 37 +++++++++++++++++++++++++++ pymc3/tests/test_math.py | 33 ++++++++++++++++++++++++ pymc3/tests/test_theanof.py | 28 --------------------- pymc3/theanof.py | 39 ++--------------------------- 5 files changed, 73 insertions(+), 66 deletions(-) create mode 100644 pymc3/tests/test_math.py diff --git a/pymc3/distributions/multivariate.py b/pymc3/distributions/multivariate.py index d1876ca2bf..0089324b34 100755 --- a/pymc3/distributions/multivariate.py +++ b/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 diff --git a/pymc3/math.py b/pymc3/math.py index fa020148b2..eb7ff7bf12 100644 --- a/pymc3/math.py +++ b/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): + 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() diff --git a/pymc3/tests/test_math.py b/pymc3/tests/test_math.py new file mode 100644 index 0000000000..f2730b8f7e --- /dev/null +++ b/pymc3/tests/test_math.py @@ -0,0 +1,33 @@ +import unittest +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)) diff --git a/pymc3/tests/test_theanof.py b/pymc3/tests/test_theanof.py index 72d0d8be29..60f2bb7a67 100644 --- a/pymc3/tests/test_theanof.py +++ b/pymc3/tests/test_theanof.py @@ -4,7 +4,6 @@ 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(): @@ -20,33 +19,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): diff --git a/pymc3/theanof.py b/pymc3/theanof.py index c7d53cdbf9..a9b7c0f494 100644 --- a/pymc3/theanof.py +++ b/pymc3/theanof.py @@ -4,9 +4,8 @@ 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 @@ -14,7 +13,7 @@ __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')) From 604b5063f8e531b4b6ad11ddc1942cc1b37a10f8 Mon Sep 17 00:00:00 2001 From: Bhargav Srinivasa Date: Sun, 19 Feb 2017 14:37:22 +0530 Subject: [PATCH 2/3] Removed unused imports, added params --- pymc3/math.py | 2 +- pymc3/tests/test_math.py | 1 - pymc3/tests/test_theanof.py | 2 -- 3 files changed, 1 insertion(+), 4 deletions(-) diff --git a/pymc3/math.py b/pymc3/math.py index eb7ff7bf12..be0904da85 100644 --- a/pymc3/math.py +++ b/pymc3/math.py @@ -42,7 +42,7 @@ def make_node(self, x): o = theano.tensor.scalar(dtype=x.dtype) return Apply(self, [x], [o]) - def perform(self, node, inputs, outputs): + def perform(self, node, inputs, outputs, params=None): try: (x,) = inputs (z,) = outputs diff --git a/pymc3/tests/test_math.py b/pymc3/tests/test_math.py index f2730b8f7e..49dd235d66 100644 --- a/pymc3/tests/test_math.py +++ b/pymc3/tests/test_math.py @@ -1,4 +1,3 @@ -import unittest import numpy as np import theano from theano.tests import unittest_tools as utt diff --git a/pymc3/tests/test_theanof.py b/pymc3/tests/test_theanof.py index 60f2bb7a67..6cbd521442 100644 --- a/pymc3/tests/test_theanof.py +++ b/pymc3/tests/test_theanof.py @@ -3,8 +3,6 @@ import numpy as np import theano from ..theanof import DataGenerator, GeneratorOp, generator -from theano.tests import unittest_tools as utt -from .helpers import SeededTest def integers(): i = 0 From d3ba445714459ca355c1be6943f5909b3393b408 Mon Sep 17 00:00:00 2001 From: Bhargav Srinivasa Date: Sun, 19 Feb 2017 19:31:52 +0530 Subject: [PATCH 3/3] Changed math style --- pymc3/distributions/multivariate.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/pymc3/distributions/multivariate.py b/pymc3/distributions/multivariate.py index 0089324b34..d2ec711174 100755 --- a/pymc3/distributions/multivariate.py +++ b/pymc3/distributions/multivariate.py @@ -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