Playing with "deconvolution"

deconv.py

import numpy as np

i = np.array(list(range(1, 50)), dtype=np.float).reshape((7, 7))
k = np.array(list(range(1, 10)), dtype=np.float).reshape((3, 3))

print("Input:"); print(i)
print("Kernel:"); print(k)

# Forward convolution. We need to pad the input so that we can read from 
# its borders. (Not doing stride etc here.)
def conv2d(i, k):
    o = np.zeros_like(i)
    i = np.pad(i, (1, 1), mode="constant")
    for oy in range(7):
        for ox in range(7):
            for ky in range(3):
                for kx in range(3):
                    o[oy, ox] += k[ky, kx] * i[oy + ky, ox + kx]
            #o[oy, ox] += bias
    return o

o = conv2d(i, k)
print("Forward convolution:")
print(o)

# The backward pass of convolution. This computes the "gradient" of the 
# convolution.
# Each value in the input comes from 3x3 = 9 different values in the original
# array, so we copy the kernel into the output (weighted by the value in the
# input array), adding it to what was already there.
def conv2d_backward(i, k):
    o = np.zeros((i.shape[0] + 2, i.shape[1] + 2))
    for oy in range(7):
        for ox in range(7):
            for ky in range(3):
                for kx in range(3):
                    o[oy + ky, ox + kx] += k[ky, kx] * i[oy, ox]
    return o

print("Convolution backward pass:")
print(conv2d_backward(o, k))

print("Convolution backward pass (without padding):")
print(conv2d_backward(o, k)[1:-1,1:-1])

# Doing convolution the other way around, using the output from the original
# convolution as the input, but this time with the kernel flipped horizontally 
# and vertically.

k_flipped = np.flipud(np.fliplr(k))
print("Convolution with flipped kernel:")
print(conv2d(o, k_flipped))

# Note that the two arrays are identical, except for the extra border in the 
# backward pass. But we can ignore those values since they're just padding values.
#
# This proves that backward convolution is identical to convolution with a
# flipped kernel. (Note that k_flipped != k.T, since .T doesn't rotate the array
# by 180 degrees.)

# Note: I left out bias here. The bias is added after the convolution, so when
# you do backwards convolution or transposed convolution you'd first subtract 
# the bias. Therefore, it doesn't really have any influence.

# Just to make sure it _really_ works, do a gradient check. We need a loss
# function of some kind that uses the input in its calculation, so just do 
# the convolution and sum up all the numbers in the convolution output.

def compute_loss(i, k):
    o = conv2d(i, k)
    return np.sum(o)

print("Loss: %f" % compute_loss(i, k))

# Now compute the numerical gradient by perturbing each element of the input
# matrix. This is dLoss/dInput (not dLoss/dWeights).
numgrad = np.zeros_like(i)
perturb = np.zeros_like(i)
epsilon = 1e-6
for r in range(i.shape[0]):
    for c in range(i.shape[1]):
        perturb[r][c] = epsilon

        loss1 = compute_loss(i - perturb, k)
        loss2 = compute_loss(i + perturb, k)

        numgrad[r][c] = (loss2 - loss1) / (2*epsilon)
        perturb[r][c] = 0
print("Numerical gradient:")
print(numgrad)

# Also do the gradient by sending all ones back through the convolution.
# This should give the same answer as the numerical gradient (up to a certain
# amount of precision).
print("Computed gradient:")
print(conv2d_backward(np.ones_like(i), k)[1:-1,1:-1])