From 95b2046104edbcacd07d273b049a784a2c8b7fb8 Mon Sep 17 00:00:00 2001 From: Bhargav Srinivasa Date: Mon, 13 Feb 2017 19:29:48 +0530 Subject: [PATCH 1/9] Adding logdet --- pymc3/theanof.py | 32 ++++++++++++++++++++++++++++++++ 1 file changed, 32 insertions(+) diff --git a/pymc3/theanof.py b/pymc3/theanof.py index e22b066754..8711138374 100644 --- a/pymc3/theanof.py +++ b/pymc3/theanof.py @@ -1,4 +1,5 @@ import numpy as np +import theano from .vartypes import typefilter, continuous_types from theano import theano, scalar, tensor as tt @@ -59,6 +60,37 @@ 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 in the Theano repo 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 = numpy.linalg.svd(x, compute_uv=False) + log_abs_det = numpy.sum(numpy.log(numpy.abs(s))) + z[0] = numpy.asarray(log_abs_det, dtype=x.dtype) + except Exception: + print('Failed to compute logabsdet 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 "LogAbsDet" + def gradient1(f, v): """flat gradient of f wrt v""" return tt.flatten(tt.grad(f, v, disconnected_inputs='warn')) From 9c034c7c9ff1f651f038937fabbf5c22391714cd Mon Sep 17 00:00:00 2001 From: Bhargav Srinivasa Date: Mon, 13 Feb 2017 21:15:10 +0530 Subject: [PATCH 2/9] Typos, added tests --- pymc3/tests/test_theanof.py | 30 +++++++++++++++++++++++++++++- pymc3/theanof.py | 21 ++++++++++++--------- 2 files changed, 41 insertions(+), 10 deletions(-) diff --git a/pymc3/tests/test_theanof.py b/pymc3/tests/test_theanof.py index 35d8c5752b..e67bb1310d 100644 --- a/pymc3/tests/test_theanof.py +++ b/pymc3/tests/test_theanof.py @@ -3,7 +3,8 @@ 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 def integers(): i = 0 @@ -18,6 +19,33 @@ def integers_ndim(ndim): yield np.ones((2,) * ndim) * i i += 1 +class TestLogDet(unittest.TestCase): + + 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 8711138374..228b02c6ef 100644 --- a/pymc3/theanof.py +++ b/pymc3/theanof.py @@ -4,8 +4,9 @@ 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 +from theano.gof import Op, Apply 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 @@ -60,12 +61,13 @@ def floatX(X): Theano derivative functions """ -class logdet(Op): +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 in the Theano repo is merged, this must be removed. + 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) @@ -76,11 +78,11 @@ def perform(self, node, inputs, outputs): try: (x,) = inputs (z,) = outputs - s = numpy.linalg.svd(x, compute_uv=False) - log_abs_det = numpy.sum(numpy.log(numpy.abs(s))) - z[0] = numpy.asarray(log_abs_det, dtype=x.dtype) + 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 logabsdet of {}.'.format(x)) + print('Failed to compute logdet of {}.'.format(x)) raise def grad(self, inputs, g_outputs): @@ -89,13 +91,14 @@ def grad(self, inputs, g_outputs): return [gz * matrix_inverse(x).T] def __str__(self): - return "LogAbsDet" + return "LogDet" + +logdet = LogDet() def gradient1(f, v): """flat gradient of f wrt v""" return tt.flatten(tt.grad(f, v, disconnected_inputs='warn')) - empty_gradient = tt.zeros(0, dtype='float32') From 5b0efc907b9daff8b0a4cb6886669cefe0de3c15 Mon Sep 17 00:00:00 2001 From: Bhargav Srinivasa Date: Mon, 13 Feb 2017 21:29:09 +0530 Subject: [PATCH 3/9] Seeded test --- pymc3/tests/test_theanof.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pymc3/tests/test_theanof.py b/pymc3/tests/test_theanof.py index e67bb1310d..72d0d8be29 100644 --- a/pymc3/tests/test_theanof.py +++ b/pymc3/tests/test_theanof.py @@ -5,6 +5,7 @@ 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 @@ -19,7 +20,7 @@ def integers_ndim(ndim): yield np.ones((2,) * ndim) * i i += 1 -class TestLogDet(unittest.TestCase): +class TestLogDet(SeededTest): def setUp(self): utt.seed_rng() From f2729f2fec8f0d51993e5ba1e4d5d845977a7eae Mon Sep 17 00:00:00 2001 From: Bhargav Srinivasa Date: Wed, 15 Feb 2017 14:42:44 +0530 Subject: [PATCH 4/9] Added to __all__, changed distribution --- pymc3/distributions/multivariate.py | 4 +++- pymc3/theanof.py | 2 +- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/pymc3/distributions/multivariate.py b/pymc3/distributions/multivariate.py index ac1a251297..c7a8f96c7c 100755 --- a/pymc3/distributions/multivariate.py +++ b/pymc3/distributions/multivariate.py @@ -12,6 +12,8 @@ from theano.tensor.nlinalg import det, matrix_inverse, trace import pymc3 as pm + +from pycm3.theanof import logdet from . import transforms from .distribution import Continuous, Discrete, draw_values, generate_samples from ..model import Deterministic @@ -113,7 +115,7 @@ def logp(self, value): delta = value - mu k = tau.shape[0] - result = k * tt.log(2 * np.pi) + tt.log(1. / det(tau)) + result = k * tt.log(2 * np.pi) - logdet(tau) result += (delta.dot(tau) * delta).sum(axis=delta.ndim - 1) return -1 / 2. * result diff --git a/pymc3/theanof.py b/pymc3/theanof.py index 228b02c6ef..c7d53cdbf9 100644 --- a/pymc3/theanof.py +++ b/pymc3/theanof.py @@ -14,7 +14,7 @@ __all__ = ['gradient', 'hessian', 'hessian_diag', 'inputvars', 'cont_inputs', 'floatX', 'jacobian', 'CallableTensor', 'join_nonshared_inputs', - 'make_shared_replacements', 'generator'] + 'make_shared_replacements', 'generator', 'LogDet', 'logdet'] def inputvars(a): From a3602bb92e0b567a42781af9a5f158a0ea9ee15f Mon Sep 17 00:00:00 2001 From: Bhargav Srinivasa Date: Thu, 16 Feb 2017 00:15:41 +0530 Subject: [PATCH 5/9] Fixed typo --- pymc3/distributions/multivariate.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pymc3/distributions/multivariate.py b/pymc3/distributions/multivariate.py index c7a8f96c7c..91a376d2d2 100755 --- a/pymc3/distributions/multivariate.py +++ b/pymc3/distributions/multivariate.py @@ -13,7 +13,7 @@ import pymc3 as pm -from pycm3.theanof import logdet +from pymc3.theanof import logdet from . import transforms from .distribution import Continuous, Discrete, draw_values, generate_samples from ..model import Deterministic From 3145e975c74fa8d30280886d192bc799c1650896 Mon Sep 17 00:00:00 2001 From: Bhargav Srinivasa Date: Fri, 17 Feb 2017 01:19:00 +0530 Subject: [PATCH 6/9] Added GPU warning --- pymc3/distributions/multivariate.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/pymc3/distributions/multivariate.py b/pymc3/distributions/multivariate.py index 91a376d2d2..e2aadc1510 100755 --- a/pymc3/distributions/multivariate.py +++ b/pymc3/distributions/multivariate.py @@ -85,14 +85,22 @@ class MvNormal(Continuous): Covariance matrix. tau : array, optional Precision matrix. + + Flags + ---------- + gpu_compat : False, because LogDet is not GPU compatible yet. + """ - def __init__(self, mu, cov=None, tau=None, *args, **kwargs): + def __init__(self, mu, cov=None, tau=None, gpu_compat=False, *args, **kwargs): super(MvNormal, self).__init__(*args, **kwargs) self.mean = self.median = self.mode = self.mu = mu = tt.as_tensor_variable(mu) tau, cov = get_tau_cov(mu, tau=tau, cov=cov) self.tau = tt.as_tensor_variable(tau) self.cov = tt.as_tensor_variable(cov) + self.gpu_compat = gpu_compat + if gpu_compat is False and theano.config.device == 'gpu': + warnings.warn("This function is not GPU compatible. Please check the gpu_compat flag") def random(self, point=None, size=None): mu, cov = draw_values([self.mu, self.cov], point=point) From ffa85bef13e0729cbdb44ad06b72ae3acaa89884 Mon Sep 17 00:00:00 2001 From: Bhargav Srinivasa Date: Fri, 17 Feb 2017 12:34:23 +0530 Subject: [PATCH 7/9] gpu_compat flag --- pymc3/distributions/multivariate.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/pymc3/distributions/multivariate.py b/pymc3/distributions/multivariate.py index e2aadc1510..ddede46764 100755 --- a/pymc3/distributions/multivariate.py +++ b/pymc3/distributions/multivariate.py @@ -89,6 +89,7 @@ class MvNormal(Continuous): Flags ---------- gpu_compat : False, because LogDet is not GPU compatible yet. + If this is set as true, the GPU compatible (but numerically unstable) log(det) is used. """ @@ -100,7 +101,7 @@ def __init__(self, mu, cov=None, tau=None, gpu_compat=False, *args, **kwargs): self.cov = tt.as_tensor_variable(cov) self.gpu_compat = gpu_compat if gpu_compat is False and theano.config.device == 'gpu': - warnings.warn("This function is not GPU compatible. Please check the gpu_compat flag") + warnings.warn("The function used is not GPU compatible. Please check the gpu_compat flag") def random(self, point=None, size=None): mu, cov = draw_values([self.mu, self.cov], point=point) @@ -123,7 +124,10 @@ def logp(self, value): delta = value - mu k = tau.shape[0] - result = k * tt.log(2 * np.pi) - logdet(tau) + if gpu_compat is False: + result = k * tt.log(2 * np.pi) - logdet(tau) + elif gpu_compat is True: + result = k * tt.log(2 * np.pi) + tt.log(1. / det(tau)) result += (delta.dot(tau) * delta).sum(axis=delta.ndim - 1) return -1 / 2. * result From 11b699071184a0b7cfc298990c41f6240e03b998 Mon Sep 17 00:00:00 2001 From: Bhargav Srinivasa Date: Fri, 17 Feb 2017 13:28:01 +0530 Subject: [PATCH 8/9] Fixed typo --- pymc3/distributions/multivariate.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pymc3/distributions/multivariate.py b/pymc3/distributions/multivariate.py index ddede46764..43389cfdb7 100755 --- a/pymc3/distributions/multivariate.py +++ b/pymc3/distributions/multivariate.py @@ -124,9 +124,9 @@ def logp(self, value): delta = value - mu k = tau.shape[0] - if gpu_compat is False: + if self.gpu_compat is False: result = k * tt.log(2 * np.pi) - logdet(tau) - elif gpu_compat is True: + elif self.gpu_compat is True: result = k * tt.log(2 * np.pi) + tt.log(1. / det(tau)) result += (delta.dot(tau) * delta).sum(axis=delta.ndim - 1) return -1 / 2. * result From 1e0d9b73513bb11934ce6b061125a67ff8742667 Mon Sep 17 00:00:00 2001 From: Bhargav Srinivasa Date: Fri, 17 Feb 2017 14:41:13 +0530 Subject: [PATCH 9/9] Changed if condition --- pymc3/distributions/multivariate.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pymc3/distributions/multivariate.py b/pymc3/distributions/multivariate.py index 43389cfdb7..d1876ca2bf 100755 --- a/pymc3/distributions/multivariate.py +++ b/pymc3/distributions/multivariate.py @@ -124,10 +124,10 @@ def logp(self, value): delta = value - mu k = tau.shape[0] - if self.gpu_compat is False: - result = k * tt.log(2 * np.pi) - logdet(tau) - elif self.gpu_compat is True: + if self.gpu_compat: result = k * tt.log(2 * np.pi) + tt.log(1. / det(tau)) + else: + result = k * tt.log(2 * np.pi) - logdet(tau) result += (delta.dot(tau) * delta).sum(axis=delta.ndim - 1) return -1 / 2. * result