Separable_kernel_constituent_vectors.ipynb

import numpy as np
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
from matplotlib import cm

'''
    This code demonstrates that a Gaussian kernel M of rank 1 is also
    separable and can be obtained as the tensor product of two vectors
    X and Y
'''


def gcd(x, y):
    if x == 0:
        return y
    else:
        return gcd(y % x, x)


# Separate a 2D kernel m into x and y vectors
# Returns [bool valid, x, y] with valid set to True if it was possible
# to find two valid constituent x and y vectors
def separate(m):
    x = np.squeeze(np.asarray(m[0, :]))
    y = np.squeeze(np.asarray(m[:, 0]))
    nx = x.size
    ny = y.size
    gx = x[0]
    for i in range(1, nx):
        gx = gcd(gx, x[i])

    gy = y[0]
    for i in range(2, ny):
        gy = gcd(gy, y[i])

    x = x / gx
    y = y / gy
    scale = m[0, 0] / (x[0] * y[0])
    x = x * scale

    valid = np.allclose(np.outer(y, x), m)
    return valid, x, y


# Usage example
# |3|                 |6  24   9|
# |4| * |2  8  3| ==  |8  32  12|
# |1|                 |2   8   3|
mat = np.outer([[3], [4], [1]], [[2, 8, 3]])
separate(mat)


# Generate a symmetric separable Gaussian kernel of radius size (w==h==2*size+1)
def generate_gaussian_kernel(size, sigma = 3, central_value = 3.5, center=None):
    # Generate a size x size grid
    a = np.linspace(-size, size, 2 * size + 1)
    b = a
    x, y = np.meshgrid(a, b)

    return x, y, np.exp(-(x / sigma)**2 - (y / sigma)**2) * central_value

xgrid, ygrid, mat = generate_gaussian_kernel(6)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(xgrid, ygrid, mat, rstride=1, cstride=1, cmap='CMRmap')
plt.show()

if np.linalg.matrix_rank(mat) == 1:
    valid, x, y = separate(mat)
    if valid:
        print "Constituent vectors found: \nv1:", x, '\nv2:', y
    else:
        print "Constituent vectors not found"
else:
    print "Not a rank 1 matrix"