Example of Self-Organizing Map (Kohonen Network) based on the Olivetti faces dataset (Cupy-based)

som-olivetti-cupy.py

import cupy as cp
import matplotlib.pyplot as plt
import numpy as np

from sklearn.datasets import fetch_olivetti_faces

# Set random seed for reproducibility
np.random.seed(1000)
cp.random.seed(1000)

# Retrieve the dataset
faces = fetch_olivetti_faces()

nb_iterations = 5000
nb_startup_iterations = 100
nb_patterns = 100
pattern_length = faces['data'].shape[1]
pattern_width = pattern_height = int(np.sqrt(pattern_length))
eta0 = 1.0
sigma0 = 10
tau = 400.0
matrix_side = 20
precomputed_distances = np.zeros((matrix_side, matrix_side, matrix_side, matrix_side))

# Create the training set
X = faces['data']
X_train = cp.array(X[0:nb_patterns])

# Initialize the weights
W = cp.random.uniform(cp.min(X_train), cp.max(X_train), size=(matrix_side, matrix_side, pattern_length))

def winning_unit(xt):
    distances = cp.linalg.norm(W - xt, ord=2, axis=2)
    max_activation_unit = cp.argmax(distances)
    return int(cp.floor(max_activation_unit / matrix_side)), max_activation_unit % matrix_side
  
def eta(t):
    return eta0 * cp.exp(-float(t) / tau)
  
def sigma(t):
    return float(sigma0) * cp.exp(-float(t) / tau)

# Precompute the distances
for i in range(matrix_side):
    for j in range(matrix_side):
        for k in range(matrix_side):
            for t in range(matrix_side):
                precomputed_distances[i, j, k, t] = np.power(float(i) - float(k), 2) + np.power(float(j) - float(t), 2)
                
precomputed_distances = cp.array(precomputed_distances)

def distance_matrix(xt, yt, sigmat):
    dm = precomputed_distances[xt, yt, :, :]
    de = 2.0 * cp.power(sigmat, 2)
    return cp.exp(-dm / de)
  
  
if __name__ == '__main__':
  # Main cycle
  sequence = np.arange(0, nb_patterns)

  for t in range(nb_iterations):
      np.random.shuffle(sequence)

      if t < nb_startup_iterations:
          etat = eta(t)
          sigmat = sigma(t)
      else:
          etat = 0.1
          sigmat = 0.01

      for n in sequence:
          x_sample = X_train[n]

          xw, yw = winning_unit(x_sample)
          dm = distance_matrix(xw, yw, sigmat)

          dW = etat * cp.expand_dims(dm, axis=2) * (x_sample - W)
          W += dW
          
  # Show the weight matrix
  matrix_w = cp.zeros((matrix_side * pattern_height, matrix_side * pattern_width))

  for i in range(matrix_side):
      for j in range(matrix_side):
          matrix_w[i * pattern_height:i * pattern_height + pattern_height, 
                   j * pattern_height:j * pattern_height + pattern_width] = W[i, j].reshape((pattern_height, pattern_width)) * 255.0
          
  fig, ax = plt.subplots(figsize=(12, 12))

  ax.matshow(matrix_w.tolist(), cmap='gray')
  ax.set_xticks([])
  ax.set_yticks([])
  
  plt.show()