Theano中的循环相关性

Question

我正在尝试用Theano计算两个信号的循环交叉相关，这将用于后续的损失计算和优化。但我不太确定该怎么做。

它的定义如下：

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

• “周期性”求和，或者说是“对于每个第k个分量”。

我可以先做普通的相关性计算，然后再进行周期性求和，但我不太清楚怎么在符号上实现这个（可能需要用到scan？）

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（以及所有与复数有关的函数，按我所知）梯度还没有实现，我想（会报错：Elemwise{real,no_inplace}.grad illegally returned an integer-valued variable. (Input index 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]