I'm trying to compute circular cross-correlation of two signals with Theano to use it in further calculation of loss that I would optimize over. But I'm not quite sure how to do that.
It is defined as following:
(f * g)[n] = sum_k f[k]g[k+n]
ccc[n] = \sum_k (f*g)[n-kN]
I could do an ordinary correlation and then perform periodic summation, but it's not quite clear how to do that (periodic summation) symbolically (using scan, probably?)
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]]) )
or use circular cross-correlation theorem and compute as a 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]) )
but in this case I can't actually compute a gradient because gradient for ifft (and all functions that has something to do with complex numbers in general, afaik) is not implemented yet, I guess (aborts with an error: Elemwise{real,no_inplace}.grad illegally returned an integer-valued variable. (Input index 0, dtype complex128))
Here's a working solution I came up with (definitely not optimal as soon as FFT is not used):
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]]) )
returns
[[ 94. 108. 108.]
[ 149. 157. 159.]]
as expected.
And can be analytically differentiated:
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.]]
here's also few more options (through direct circulant calculation) and time-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
which gives
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]