定制Theano Op进行数值积分

Question

我正在尝试编写一个自定义的Theano Op,它在数值上集成了两个值之间的函数.Op是PyMC3的自定义可能性,涉及一些积分的数值计算.我不能简单地使用@as_op装饰器,因为我需要使用HMC来执行MCMC步骤.任何帮助将不胜感激,因为这个问题似乎已经走到了好几次,但一直没有得到解决(如/sf/ask/2579711081/,Theano:实现一个积分函数).

显然,一种解决方案是在Theano中编写一个数值积分器,但是当非常好的积分器已经可用时,这似乎是浪费精力,例如通过scipy.integrate.

为了保持这个最小的例子,让我们尝试在Op中集成0到1之间的函数.以下在Op之外集成了Theano函数,并且在我的测试已经消失的情况下产生了正确的结果.

import theano
import theano.tensor as tt
from scipy.integrate import quad

x = tt.dscalar('x')
y = x**4 # integrand
f = theano.function([x], y)

print f(0)
print f(1)

ans = integrate.quad(f, 0, 1)[0]

print ans

但是,尝试在Op中进行集成似乎要困难得多.我目前的最大努力是:

import numpy as np
import theano
import theano.tensor as tt
from scipy import integrate

class IntOp(theano.Op):
    __props__ = ()

    def make_node(self, x):
        x = tt.as_tensor_variable(x)
        return theano.Apply(self, [x], [x.type()])

    def perform(self, node, inputs, output_storage):
        x = inputs[0]
        z = output_storage[0]

        f_to_int = theano.function([x], x)
        z[0] = tt.as_tensor_variable(integrate.quad(f_to_int, 0, 1)[0])

    def infer_shape(self, node, i0_shapes):
        return i0_shapes

    def grad(self, inputs, output_grads):
        ans = integrate.quad(output_grads[0], 0, 1)[0]
        return [ans]

intOp = IntOp()

x = tt.dmatrix('x')
y = intOp(x)

f = theano.function([x], y)

inp = np.asarray([[2, 4], [6, 8]], dtype=theano.config.floatX)
out = f(inp)

print inp
print out

这给出了以下错误:

Traceback (most recent call last):
  File "stackoverflow.py", line 35, in <module>
    out = f(inp)
  File "/usr/local/lib/python2.7/dist-packages/theano/compile/function_module.py", line 871, in __call__
    storage_map=getattr(self.fn, 'storage_map', None))
  File "/usr/local/lib/python2.7/dist-packages/theano/gof/link.py", line 314, in raise_with_op
    reraise(exc_type, exc_value, exc_trace)
  File "/usr/local/lib/python2.7/dist-packages/theano/compile/function_module.py", line 859, in __call__
    outputs = self.fn()
  File "/usr/local/lib/python2.7/dist-packages/theano/gof/op.py", line 912, in rval
    r = p(n, [x[0] for x in i], o)
  File "stackoverflow.py", line 17, in perform
    f_to_int = theano.function([x], x)
  File "/usr/local/lib/python2.7/dist-packages/theano/compile/function.py", line 320, in function
    output_keys=output_keys)
  File "/usr/local/lib/python2.7/dist-packages/theano/compile/pfunc.py", line 390, in pfunc
    for p in params]
  File "/usr/local/lib/python2.7/dist-packages/theano/compile/pfunc.py", line 489, in _pfunc_param_to_in
    raise TypeError('Unknown parameter type: %s' % type(param))
TypeError: Unknown parameter type: <type 'numpy.ndarray'>
Apply node that caused the error: IntOp(x)
Toposort index: 0
Inputs types: [TensorType(float64, matrix)]
Inputs shapes: [(2, 2)]
Inputs strides: [(16, 8)]
Inputs values: [array([[ 2.,  4.],
       [ 6.,  8.]])]
Outputs clients: [['output']]

Backtrace when the node is created(use Theano flag traceback.limit=N to make it longer):
  File "stackoverflow.py", line 30, in <module>
    y = intOp(x)
  File "/usr/local/lib/python2.7/dist-packages/theano/gof/op.py", line 611, in __call__
    node = self.make_node(*inputs, **kwargs)
  File "stackoverflow.py", line 11, in make_node
    return theano.Apply(self, [x], [x.type()])

HINT: Use the Theano flag 'exception_verbosity=high' for a debugprint and storage map footprint of this apply node.

我对此感到惊讶,尤其是TypeError,因为我认为我已经将output_storage变量转换为张量,但它似乎相信它仍然是一个ndarray.

Answer 1

我找到了你的问题,因为我试图在PyMC3中构建一个随机变量,代表一般的点过程(Hawkes,Cox,Poisson等),而似然函数有一个积分.我真的希望能够使用哈密尔顿蒙特卡罗或NUTS采样器,所以我需要在时间方面不可分割.

从你的尝试开始,我做了一个integOut theano Op,它似乎与我需要的行为一起正常工作.我已经在几个不同的输入上测试了它(不是我的统计模型,但它看起来很有希望!).我总是theano n00b,所以请原谅任何愚蠢.如果有人有任何反馈意见,我将非常感谢.不确定它正是您正在寻找的,但这是我的解决方案(示例位于底部和doc字符串中).*编辑:简化了一些残障的方法来解决这个问题.

import theano
import theano.tensor as T
from scipy.integrate import quad

class integrateOut(theano.Op):
    """
    Integrate out a variable from an expression, computing
    the definite integral w.r.t. the variable specified
    !!! Only implemented in this for scalars !!!


    Parameters
    ----------
    f : scalar
        input 'function' to integrate
    t : scalar
        the variable to integrate out
    t0: float
        lower integration limit
    tf: float
        upper integration limit

    Returns
    -------
    scalar
        a new scalar with the 't' integrated out

    Notes
    -----

    usage of this looks like:
    x = T.dscalar('x')
    y = T.dscalar('y')
    t = T.dscalar('t')

    z = (x**2 + y**2)*t

    # integrate z w.r.t. t as a function of (x,y)
    intZ = integrateOut(z,t,0.0,5.0)(x,y)
    gradIntZ = T.grad(intZ,[x,y])

    funcIntZ = theano.function([x,y],intZ)
    funcGradIntZ = theano.function([x,y],gradIntZ)

    """
    def __init__(self,f,t,t0,tf,*args,**kwargs):
        super(integrateOut,self).__init__()
        self.f = f
        self.t = t
        self.t0 = t0
        self.tf = tf

    def make_node(self,*inputs):
        self.fvars=list(inputs)
        # This will fail when taking the gradient... don't be concerned
        try:
            self.gradF = T.grad(self.f,self.fvars)
        except:
            self.gradF = None
        return theano.Apply(self,self.fvars,[T.dscalar().type()])

    def perform(self,node, inputs, output_storage):
        # Everything else is an argument to the quad function
        args = tuple(inputs)
        # create a function to evaluate the integral
        f = theano.function([self.t]+self.fvars,self.f)
        # actually compute the integral
        output_storage[0][0] = quad(f,self.t0,self.tf,args=args)[0]

    def grad(self,inputs,grads):
        return [integrateOut(g,self.t,self.t0,self.tf)(*inputs)*grads[0] \
            for g in self.gradF]

x = T.dscalar('x')
y = T.dscalar('y')
t = T.dscalar('t')

z = (x**2+y**2)*t

intZ = integrateOut(z,t,0,1)(x,y)
gradIntZ = T.grad(intZ,[x,y])
funcIntZ = theano.function([x,y],intZ)
funcGradIntZ = theano.function([x,y],gradIntZ)
print funcIntZ(2,2)
print funcGradIntZ(2,2)