Как поместить «произвольную» операцию в скользящее окно с помощью Theano?

Question

Я хочу определить некоторую функцию на матрице X. Например, mean(pow(X - X0, 2)), где X0 - другая матрица (X0 фиксирована/постоянна). Чтобы сделать его более конкретным, допустим, что и X, и X0 - это матрицы 10 x 10. Результатом операции является реальное число.

Теперь у меня большая матрица (допустим, 500 x 500). Я хочу применить описанную выше операцию ко всем подматрицам 10 x 10 «большой» матрицы. Другими словами, я хочу сдвинуть окно 10 x 10 над «большой» матрицей. Для каждого местоположения окна я должен получить реальное число. Итак, в качестве конечного результата мне нужно получить вещественную матрицу (или 2D-тензор) (ее форма должна быть 491 x 491).

То, что я хочу иметь, близко к сверточному слою, но не совсем то же самое, потому что я хочу использовать среднее квадратическое отклонение вместо линейной функции, представленной нейроном.

Answer 1

Это всего лишь решение для работы с цифрами, надеясь, что этого достаточно. Я предполагаю, что ваша функция состоит из операции над матричными элементами и среднего значения, то есть масштабированной суммы. Следовательно, достаточно посмотреть на Y как в

Y = np.power(X-X0, 2) 

Так что нам нужно иметь дело только с определением окошками среднего. Заметим, что для 1-го случая матричный продукт с соответствующим вектором единиц может быть определен для вычисления среднего значения, например.

h = np.array([0, 1, 1, 0]) # same dimension as y 
m1 = np.dot(h, y)/2 
m2 = (y[1] + y[2])/2 
print(m1 == m2) # True 

Двухмерный корпус аналогичен, но имеет два матричных умножения, один для строк и один для столбцов. Например.

m_2 = np.dot(np.dot(h, Y), h)/2**2 

Для построения скользящего окна нам необходимо построить матрицу сдвинутых окон, например.

H = [[1, 1, 1, 0, 0, ..., 0], 
    [0, 1, 1, 1, 0, ..., 0], 
      . 
      . 
      . 
    [0, ..., 0, 0, 1, 1, 1]] 

вычислить все суммы

S = np.dot(np.dot(H, Y), H.T) 

Полный примером для матрицы с (m, m) окна будет

import numpy as np 

n, m = 500, 10 
X0 = np.ones((n, n)) 
X = np.random.rand(n, n) 
Y = np.power(X-X0, 2) 

h = np.concatenate((np.ones(m), np.zeros(n-m))) # window at position 0 
H = np.vstack((np.roll(h, k) for k in range(n+1-m))) # slide the window 
M = np.dot(np.dot(H,Y), H.T)/m**2 # calculate the mean 
print(M.shape) # (491, 491) 

Альтернативными, но, вероятно, несколько менее эффективным способом для построения H is

H = np.sum(np.diag(np.ones(n-k), k)[:-m+1, :] for k in range(m)) 

Обновление

Вычисление среднего квадрата отклонения также возможно с этим подходом. Для этого, мы обобщаем вектор идентичность |x-x0|^2 = (x-x0).T (x-x0) = x.T x - 2 x0.T x + x0.T x0 (пространство обозначает скалярную или матричное умножение и .T транспонированной вектор) в матричном случае:

Мы предполагаем W является (m,n) матрицей, содержащей единичную матрицы блока (m.m), который способен для извлечения подматрицы (k0,k1)-ой (m,m) на Y = W Z W.T, где Z представляет собой матрицу (n,n), содержащую данные.Расчет разницы

D = Y - X0 = Y = W Z W.T - X0 

просто, где X0 и D является (m,m) матрицей. Квадратный корень квадратной суммы элементов называется Frobenius norm. На основании этих identities, можно записать как квадрат суммы

s = sum_{i,j} D_{i,j}^2 = trace(D.T D) = trace((W Z W.T - X0).T (H Z H.T - X0)) 
    = trace(W Z.T W.T W Z W.T) - 2 trace(X0.T W Z W.T) + trace(X0.T X0) 
    =: Y0 + Y1 + Y2 

Термин Y0 можно интерпретировать как H Z H.T от метода от above.The термина Y1 может быть интерпретирован как взвешенное среднее значение по Z и Y2 является константой , который необходимо определить только один раз. Таким образом, возможная реализация будет:

import numpy as np 

n, m = 500, 10 
x0 = np.ones(m) 
Z = np.random.rand(n, n) 

Y0 = Z**2 
h0 = np.concatenate((np.ones(m), np.zeros(n-m))) 
H0 = np.vstack((np.roll(h0, k) for k in range(n+1-m))) 
M0 = np.dot(np.dot(H0, Y0), H0.T) 

h1 = np.concatenate((-2*x0, np.zeros(n-m))) 
H1 = np.vstack((np.roll(h1, k) for k in range(n+1-m))) 
M1 = np.dot(np.dot(H1, Z), H0.T) 

Y2 = np.dot(x0, x0) 
M = (M0 + M1)/m**2 + Y2