python - Theano 中的循环相关

Question

我正在尝试计算 circular cross-correlation两个信号与 Theano 一起使用它来进一步计算我将优化的损失。但我不太确定该怎么做。

定义如下:

(f * g)[n] = sum_k f[k]g[k+n]
ccc[n] = \sum_k (f*g)[n-kN]

• “周期性”求和或类似“针对每个第 k 个分量”。

我可以做一个普通的相关，然后执行周期求和，但不太清楚如何以符号方式(可能使用扫描？)做那个(周期求和)

conv2d = T.signal.conv.conv2d

x = T.dmatrix()
y = T.dmatrix()
veclen = x.shape[1]

corr_expr = conv2d(x, y[:, ::-1], image_shape=(1, veclen), border_mode='full')

# circ_corr = T.sum([corr_expr[k::veclen] for k in T.arange(veclen)])

corr = theano.function([x, y], outputs=circ_corr)

corr( np.array([[2, 3, 5]]), np.array([[7, 11, 13]]) )

或使用循环互相关定理并计算为 iFFT(FFT(x)*FFT(y)):

import theano.sandbox.fourier as dft
x = T.dmatrix()
y = T.dvector()
veclen = x.shape[1]

exp = T.real( 
        dft.ifft( 
            dft.fft(x, veclen, axis=1) 
            * dft.fft(y[::-1], y.shape[0], axis=1).reshape((1, -1)), 
            veclen, axis=1
        ) 
    )[:, ::-1]
f = theano.function([x, y], outputs=exp)

f(np.array([[2, 3, 5], [3, 4, 4], [5, 6, 7]]), np.array([7, 11, 13]) )

但在这种情况下，我实际上无法计算梯度，因为 ifft 的梯度(以及所有与一般复数有关的函数，afaik)尚未实现，我猜(因错误而中止:Elemwise{real,no_inplace}.grad 非法返回一个整数值变量。(输入索引 0，dtype complex128))

Answer 1

这是我想出的一个可行的解决方案(如果不使用 FFT 就绝对不是最佳解决方案):

def circular_crosscorelation(X, y):
    """ 
    Input:
        symbols for X [n, m]
        and y[m,]

    Returns: 
        symbol for circular cross corelation of each of row in X with 
        cc[n, m]
    """
    n, m = X.shape
    corr_expr = T.signal.conv.conv2d(X, y[::-1].reshape((1, -1)), image_shape=(1, m), border_mode='full')
    corr_len = corr_expr.shape[1]
    pad = m - corr_len%m
    v_padded = T.concatenate([corr_expr, T.zeros((n, pad))], axis=1)
    circ_corr_exp = T.sum(v_padded.reshape((n, v_padded.shape[1] / m, m)), axis=1)
    return circ_corr_exp[:, ::-1]

X = T.dmatrix()
y = T.dmatrix()
cc = theano.function([X, y], circular_crosscorelation(X, y))
print cc( np.array([[2, 3, 5], [4, 5, 6]]), np.array([[7, 11, 13]]) )

返回

[[  94.  108.  108.]
 [ 149.  157.  159.]]

正如预期的那样。

并且可以解析微分:

score = T.sum(circ_corr_exp**2)
grad = T.grad(score, x)
g = theano.function([x, y], outputs=grad)
print g( np.array([[2, 3, 5], [4, 5, 6]]), np.array([[7, 11, 13]]) )

>> [[ 6332.  6388.  6500.]
>>  [ 9554.  9610.  9666.]]

---

这里还有更多选项(通过直接循环计算)和时间比较:

def circulant_np(v):
    row = np.arange(len(v))
    col = -np.arange(len(v))
    idx = (row[:, np.newaxis] + col)%len(v)
    return v[idx]

print circulant_np(np.array([1, 2, 3, 5]))

def c_corr_np(a, b):
    return circulant_np(a).dot(b[::-1])

def circulant_t(v):
    row = T.arange(v.shape[0])
    col = -T.arange(v.shape[0])
    idx = (row.reshape((-1, 1)) + col)%v.shape[0]
    return v[idx]

def c_corr_t_f(a, b):
    """ 1d correlation using circulant matrix """
    return circulant_t(a).dot(b[::-1])

a = T.dvector('a')
b = T.dvector('b')
c_corr_t = theano.function([a, b], c_corr_t_f(a, b))

print c_corr_np(np.array([2, 3, 5]), np.array([7, 11, 13]))
print c_corr_t(np.array([2, 3, 5]), np.array([7, 11, 13]))
print c_corr( np.array([[2, 3, 5]]), np.array([[7, 11, 13]]) )

%timeit c_corr_np(np.array([2, 3, 5]), np.array([7, 11, 13]))
%timeit c_corr_t(np.array([2, 3, 5]), np.array([7, 11, 13]))
%timeit c_corr( np.array([[2, 3, 5]]), np.array([[7, 11, 13]]) ) # = circular_crosscorelation

给出

10000 loops, best of 3: 30.6 µs per loop
10000 loops, best of 3: 132 µs per loop
10000 loops, best of 3: 149 µs per loop

和逆交叉校正:

def inverse_circular_crosscorelation(y):
    """ 
    Input:
        symbol for y[1, m]

    Returns: 
        symbol for y_inv s.t. 
        cc( y, y_inv ) = (1, 0 ... 0)
    """

    A = circulant_t(y.reshape((-1, )))
    b = T.concatenate([T.zeros((y.shape[1] - 1, )), T.ones((1, ))]).reshape((-1, 1))
    return T.nlinalg.matrix_inverse(A).dot(b).reshape((1, -1))[:, ::-1]