Циркулярная корреляция в Theano

Question

Я пытаюсь вычислить circular cross-correlation двух сигналов с помощью Theano, чтобы использовать его для дальнейшего расчета потерь, которые я бы оптимизировал. Но я не совсем уверен, как это сделать.

Она определяется следующим образом:

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

• "периодическим" суммирование или как "для каждого к-го компонента".

я мог бы сделать обычную корреляцию, а затем выполнять периодическое суммирование, но это не совсем понятно, как это сделать (периодическое суммирование) символически (с помощью сканирования, возможно?)

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

или использовать круговой крест декорреляции теоремы и вычислить как Иффт (FFT (х) * FFT (у)):

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

, но в этом случае я не могу на самом деле вычислить градиент, поскольку градиент для ОБПФА (и все функции, которые есть что-то делать с комплексными числами в гене аль, AFAIK) еще не реализована, я думаю (прерывается с ошибкой: 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.]] 

---

здесь также несколько вариантов (путем прямого циркулянтного расчета) и время Comparation:

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]