Theano中的圓形相關

Question

我試圖用Theano計算兩個信號的circular cross-correlation，以便進一步計算我將優化的損耗。但我不太確定如何做到這一點。

它被定義爲以下：

(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]])) 

或使用圓形截面-correlation定理和計算的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 

and inverse cross-corr：

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]