SpectralNet (MNIST data) の実装に関するメモ

##ポイント
* SpectralNet を実装し、具体的な数値例で確認。

##レファレンス
[1. SpectralNet: Spectral Clustering using Deep Neural Networks](https://arxiv.org/abs/1801.01587)

![image.png](https://qiita-image-store.s3.amazonaws.com/0/254604/d60ff1d0-89ed-8fb7-0d2a-712cbf2daf8e.png)

　　　　　　　　　　　　　（参照論文より引用）

##データ　
```python
from sklearn.manifold import TSNE

n_samples = 100
x_sample_train = mnist.train.images[:n_samples]
labels_train = np.argmax(mnist.train.labels[:n_samples], axis = 1)

x_embedded = TSNE(n_components=2).fit_transform(x_sample_train)

plt.scatter(x_embedded[:,0], x_embedded[:, 1], c = labels_train, cmap = 'Accent')
plt.colorbar()
plt.show()
```
![image.png](https://qiita-image-store.s3.amazonaws.com/0/254604/54e7b1be-55d7-b122-1c20-0b41686698e1.png)

##サンプルコード (K-means)
```python
tf.reset_default_graph()

# K-means clustering
def k_means(x, n_clusters):
  r_seed = np.random.randint(100)
  graph = tf.contrib.factorization.KMeans(x, num_clusters = n_clusters, \
                                         random_seed = r_seed)

  return graph.training_graph()


n_clusters = 10
x_input = x_sample_train.astype(np.float32)

all_scores, cluster_idx, scores, cluster_centers_initialized, init_op, training_op = k_means(x_input, n_clusters)
    
init_var = tf.global_variables_initializer()

with tf.Session() as sess:
  sess.run(init_var)
  sess.run(init_op)
  cluster_idx_init = sess.run(cluster_idx)
  
  for i in range(30):
    sess.run(training_op)
    
  cluster_idx_last = sess.run(cluster_idx)

x_embedded = TSNE(n_components=2).fit_transform(x_input)

plt.scatter(x_embedded[:,0], x_embedded[:, 1], c = cluster_idx_last[0], cmap = 'Accent')
plt.title('K-means clustering')
plt.colorbar()
plt.show()
```
![image.png](https://qiita-image-store.s3.amazonaws.com/0/254604/693dab06-56c4-1ad3-bdad-f82669c30494.png)

## サンプルコード (Spectral Clustering)
```python
tf.reset_default_graph()

# Spectral Clustering
def squared_distance(x, n_input):
  X = tf.expand_dims(x, 1)
  M = tf.tile(X, [1, n_input, 1])
  M_T = tf.transpose(M, [1, 0, 2])

  return tf.reduce_sum(tf.square(M - M_T), axis = -1)

def full_affinity(x, n_input, scale = 1.0):
  S = squared_distance(x, n_input)

  return tf.exp(-0.5 * S / scale**2)

def knn_affinity(x, n_input, n_nbrs):
  S = squared_distance(x, n_input)

  nn = tf.nn.top_k(-S, n_nbrs, sorted = True)
  nn_var = nn[0]
  nn_idx = nn[1]

  jj = tf.reshape(nn_idx, [-1, 1])
  ii = tf.range(n_input)
  ii = tf.reshape(ii, [-1, 1])
  ii = tf.tile(ii, [1, n_nbrs])
  ii = tf.reshape(ii, [-1, 1])
  indices = tf.concat([ii, jj], axis = 1)

  values = tf.ones(shape = [n_input * n_nbrs], dtype = tf.float32)
  mask = tf.SparseTensor(indices = tf.cast(indices, tf.int64), values = values, dense_shape = [n_input, n_input])
  mask = tf.sparse_reorder(mask)
  mask = tf.sparse_tensor_to_dense(mask)

  scale_squared = tf.reduce_mean(-nn_var[:, n_nbrs - 1])

  W = tf.exp(-0.5 * S / scale_squared) * mask
    
  W = 0.5 * (W + tf.transpose(W))
  
  return W

def laplacian(W):

  # unnormalized Laplacian
  diag = tf.reduce_sum(W, axis = -1)
  D = tf.matrix_diag(diagonal = diag)
  L = D - W 
  e, v = tf.self_adjoint_eig(L)

  # normalized Laplacian
  D_I = tf.matrix_diag(diagonal = 1.0 / diag)
  L_rw = tf.matmul(D_I, L)
  e_rw, v_rw = tf.self_adjoint_eig(L_rw)

  return v, v_rw

def k_means(x, n_clusters):
  r_seed = np.random.randint(100)
  graph = tf.contrib.factorization.KMeans(x, num_clusters = n_clusters, \
                                         random_seed = r_seed)

  return graph.training_graph()

n_clusters = 10
n_nbrs = 5
x_input = x_sample_train.astype(np.float32)

#W = full_affinity(x_input, n_sample, 1)
W = knn_affinity(x_input, n_samples, n_nbrs)
v, v_rw = laplacian(W)
y = v[:, :n_clusters]
all_scores, cluster_idx, scores, cluster_centers_initialized, init_op, training_op = k_means(y, n_clusters)
    
init_var = tf.global_variables_initializer()

with tf.Session() as sess:
  sess.run(init_var)
  sess.run(init_op)
  cluster_idx_init = sess.run(cluster_idx)
  
  for i in range(30):
    sess.run(training_op)

  y_last = sess.run(y)
  cluster_idx_last = sess.run(cluster_idx)
  
y_embedded = TSNE(n_components=2).fit_transform(y_last)

plt.scatter(y_embedded[:,0], y_embedded[:, 1], c = cluster_idx_last[0], cmap = 'Accent')
plt.title('Spectral clustering')
plt.colorbar()
plt.show()

plt.scatter(y_embedded[:,0], y_embedded[:, 1], c = labels_train, cmap = 'Accent')
plt.title('Labels')
plt.colorbar()
plt.show()
```
![image.png](https://qiita-image-store.s3.amazonaws.com/0/254604/6f6a4368-7584-7d5d-1581-3fc3964c01e6.png)

##サンプルコード (SpectralNet)
```python
from sklearn.manifold import TSNE

n_samples = 1000
x_sample_train = mnist.train.images[:n_samples]
labels_train = np.argmax(mnist.train.labels[:n_samples], axis = 1)
x_sample_test = mnist.test.images[:n_samples]
labels_test = np.argmax(mnist.test.labels[:n_samples], axis = 1)

#x_embedded_train = TSNE(n_components=2).fit_transform(x_sample_train)
#x_embedded_test = TSNE(n_components=2).fit_transform(x_sample_test)

#fig = plt.figure(figsize = (10, 3))

#ax1 = fig.add_subplot(1, 2, 1)
#ax1.scatter(x_embedded_train[:,0], x_embedded_train[:, 1], c = labels_train, cmap = 'Accent')
#ax1.set_title('Train')

#ax2 = fig.add_subplot(1, 2, 2)
#ax2.scatter(x_embedded_test[:,0], x_embedded_test[:, 1], c = labels_test, cmap = 'Accent')
#ax2.set_title('Test')

#plt.show()

tf.reset_default_graph()

# SpectralNet

# Utility functions
def squared_distance(x, n_input):
  X = tf.expand_dims(x, 1)
  M = tf.tile(X, [1, n_input, 1])
  M_T = tf.transpose(M, [1, 0, 2])
  return tf.reduce_sum(tf.square(M - M_T), axis = -1)

#def full_affinity(x, n_input, scale = 1.0):
#  S = squared_distance(x, n_input)
#  return tf.exp(-0.5 * S / scale**2)

def knn_affinity(x, n_input, n_nbrs):
  S = squared_distance(x, n_input)

  nn = tf.nn.top_k(-S, n_nbrs, sorted = True)
  nn_var = nn[0]
  nn_idx = nn[1]

  jj = tf.reshape(nn_idx, [-1, 1])
  ii = tf.range(n_input)
  ii = tf.reshape(ii, [-1, 1])
  ii = tf.tile(ii, [1, n_nbrs])
  ii = tf.reshape(ii, [-1, 1])
  indices = tf.concat([ii, jj], axis = 1)

  values = tf.ones(shape = [n_input * n_nbrs], dtype = tf.float32)
  n_input_int64 = tf.cast(n_input, tf.int64)
  mask = tf.SparseTensor(indices = tf.cast(indices, tf.int64), values = values, \
                         dense_shape = [n_input_int64, n_input_int64])
  mask = tf.sparse_reorder(mask)
  mask = tf.sparse_tensor_to_dense(mask)

  scale_squared = tf.reduce_mean(-nn_var[:, n_nbrs - 1])
  W = tf.exp(-0.5 * S / scale_squared) * mask
  return 0.5 * (W + tf.transpose(W))
  
#def laplacian(W):
#  # unnormalized Laplacian
#  diag = tf.reduce_sum(W, axis = -1)
#  D = tf.matrix_diag(diagonal = diag)
#  L = D - W 
#  e, v = tf.self_adjoint_eig(L)

#  # normalized Laplacian
#  D_I = tf.matrix_diag(diagonal = 1.0 / diag)
#  L_rw = tf.matmul(D_I, L)
#  e_rw, v_rw = tf.self_adjoint_eig(L_rw)

#  return v, v_rw

def k_means(x, n_clusters):
  r_seed = np.random.randint(100)
  graph = tf.contrib.factorization.KMeans(x, num_clusters = n_clusters, \
                                         random_seed = r_seed)
  return graph.training_graph()

def orthogonalization_weights(x, n_input):
  x_2 = tf.matmul(tf.transpose(x), x)
  L = tf.cholesky(x_2)
  return tf.transpose(tf.matrix_inverse(L)) * tf.sqrt(tf.cast(n_input, tf.float32))

# Neural Net
def weight_variable(name, shape):
  initializer = tf.truncated_normal_initializer(mean = 0.0, stddev = 0.01, dtype = tf.float32)
  return tf.get_variable(name, shape, initializer = initializer)

def bias_variable(name, shape):
  initializer = tf.constant_initializer(value = 0.0, dtype = tf.float32)
  return tf.get_variable(name, shape, initializer = initializer)

def beta_variable(name, shape):
  initializer = tf.constant_initializer(value = 0.0, dtype = tf.float32)
  return tf.get_variable(name, shape, initializer = initializer)

def gamma_variable(name, shape):
  initializer = tf.constant_initializer(value = 1.0, dtype = tf.float32)
  return tf.get_variable(name, shape, initializer = initializer)

def mlp(x, n_in, n_units_1, n_units_2, n_clusters):
  with tf.variable_scope('enc_dec'):
    w_1 = weight_variable('w_1', [n_in, n_units_1])
    b_1 = bias_variable('b_1', [n_units_1])

    h_1 = tf.add(tf.matmul(x, w_1), b_1)

    # batch normalization
    beta_1 = beta_variable('beta_1', [n_units_1])
    gamma_1 = gamma_variable('gamma_1', [n_units_1])

    mean_1, var_1 = tf.nn.moments(h_1, [0])
    h_1_n = gamma_1 * (h_1 - mean_1) / tf.sqrt(var_1 + 1e-5) + beta_1

    bn_1 = tf.nn.relu(h_1_n)

    w_2 = weight_variable('w_2', [n_units_1, n_units_2])
    b_2 = bias_variable('b_2', [n_units_2])

    h_2 = tf.nn.sigmoid(tf.add(tf.matmul(bn_1, w_2), b_2))

    # batch normalization
    beta_2 = beta_variable('beta_2', [n_units_1])
    gamma_2 = gamma_variable('gamma_2', [n_units_1])

    mean_2, var_2 = tf.nn.moments(h_2, [0])
    h_2_n = gamma_2 * (h_2 - mean_2) / tf.sqrt(var_2 + 1e-5) + beta_2

    bn_2 = tf.nn.relu(h_2_n)

    w_3 = weight_variable('w_3', [n_units_2, n_clusters])
    b_3 = bias_variable('b_3', [n_clusters])

    return tf.add(tf.matmul(bn_2, w_3), b_3)

def loss(y, n_input, W):
  S = squared_distance(y, n_input)
  n_input_float32 = tf.cast(n_input, tf.float32)
  return tf.reduce_mean(W * S) / n_input_float32**2
  
def training(loss, learning_rate):
  optimizer = tf.train.AdamOptimizer(learning_rate = learning_rate)
  train_step = optimizer.minimize(loss)
  return train_step

#def training_clipped(loss, learning_rate, clip_norm):
#  optimizer = tf.train.AdamOptimizer(learning_rate = learning_rate)

#  grads_and_vars = optimizer.compute_gradients(loss)
#  clipped_grads_and_vars = [(tf.clip_by_norm(grad, clip_norm = clip_norm), \
#                           var) for grad, var in grads_and_vars]
#  train_step = optimizer.apply_gradients(clipped_grads_and_vars)
#  return train_step

dim = 28*28
n_units_1 = 128
n_units_2 = 128
n_clusters = 10
n_nbrs = 5

batch_size = 50
learning_rate = 0.01
n_iter = 1000
show_step = 200
n_iter_km = 30

x = tf.placeholder(shape = [None, dim], dtype = tf.float32)
n = tf.placeholder(shape = (), dtype = tf.int32)

y_tilde = mlp(x, dim, n_units_1, n_units_2, n_clusters)
ortho_weights = orthogonalization_weights(y_tilde, n)
y = tf.matmul(y_tilde, ortho_weights)
W = knn_affinity(x, n, n_nbrs)

loss = loss(y, n, W)
train_step = training(loss, learning_rate)

# K-means
all_scores, cluster_idx, scores, cluster_centers_initialized, init_op, training_op = k_means(y, n_clusters)

init_var = tf.global_variables_initializer()

with tf.Session() as sess:
  sess.run(init_var)
  
  x_input_train = x_sample_train.astype(np.float32)
  x_input_test = x_sample_test.astype(np.float32)
  
  history_loss_train = []
  history_loss_test = []
  
  for i in range(n_iter):
    # Train
    rand_index = np.random.choice(len(x_input_train), size = batch_size)
    x_batch = x_input_train[rand_index]
    feed_dict = {x: x_batch, n: batch_size}
    sess.run(train_step, feed_dict = feed_dict)

    temp_loss = sess.run(loss, feed_dict = feed_dict)
    history_loss_train.append(temp_loss)
    if (i + 1) % show_step == 0:
      print ('--------------------')
      print ('Iteration: ' + str(i + 1) + '  Loss: ' + str(temp_loss))

    # Test
    rand_index = np.random.choice(len(x_input_test), size = batch_size)
    x_batch = x_input_test[rand_index]
    feed_dict = {x: x_batch, n: batch_size}
    temp_loss = sess.run(loss, feed_dict = feed_dict)
    history_loss_test.append(temp_loss)
      
  # K-means clustering
  # Train
  feed_dict_km = {x: x_sample_train, n: n_samples}
  sess.run(init_op, feed_dict = feed_dict_km)
  for i in range(n_iter_km):
    sess.run(training_op, feed_dict = feed_dict_km)
  y_last_train = sess.run(y, feed_dict = feed_dict_km)  
  cluster_idx_last_train = sess.run(cluster_idx, feed_dict = feed_dict_km)

  # Test
  feed_dict_km = {x: x_sample_test, n: n_samples}
  sess.run(init_op, feed_dict = feed_dict_km)
  for i in range(n_iter_km):
    sess.run(training_op, feed_dict = feed_dict_km)
  y_last_test = sess.run(y, feed_dict = feed_dict_km)  
  cluster_idx_last_test = sess.run(cluster_idx, feed_dict = feed_dict_km)

  
fig = plt.figure(figsize = (10, 3))
ax1 = fig.add_subplot(1, 2, 1)
ax1.plot(history_loss_train)
ax1.set_title('Loss (train)')
ax2 = fig.add_subplot(1, 2, 2)
ax2.plot(history_loss_test)
ax2.set_title('Loss (test)')
plt.show()

y_embedded_train = TSNE(n_components=2).fit_transform(y_last_train)
y_embedded_test = TSNE(n_components=2).fit_transform(y_last_test)

plt.scatter(y_embedded_train[:,0], y_embedded_train[:, 1], c = cluster_idx_last_train[0], cmap = 'Accent')
plt.title('Spectral clustering (train)')
plt.colorbar()
plt.show()

plt.scatter(y_embedded_train[:,0], y_embedded_train[:, 1], c = labels_train, cmap = 'Accent')
plt.title('Labels (train)')
plt.colorbar()
plt.show()

plt.scatter(y_embedded_test[:,0], y_embedded_test[:, 1], c = cluster_idx_last_test[0], cmap = 'Accent')
plt.title('Spectral clustering (test)')
plt.colorbar()
plt.show()

plt.scatter(y_embedded_test[:,0], y_embedded_test[:, 1], c = labels_train, cmap = 'Accent')
plt.title('Labels (test)')
plt.colorbar()
plt.show()
  

#fig = plt.figure(figsize = (10, 3))
#ax1 = fig.add_subplot(1, 2, 1)
#ax1.scatter(x_sample_train[:,0], x_sample_train[:, 1], c = cluster_idx_last_train[0], cmap = 'RdBu')
#ax1.set_title('Clustering (train)')
#ax2 = fig.add_subplot(1, 2, 2)
#ax2.scatter(x_sample_test[:,0], x_sample_test[:, 1], c = cluster_idx_last_test[0], cmap = 'RdBu')
#ax2.set_title('Clustering (test)')
#plt.show()
```
![image.png](https://qiita-image-store.s3.amazonaws.com/0/254604/73729a37-a3fb-cfd8-bf25-0595d425b8b5.png)

![image.png](https://qiita-image-store.s3.amazonaws.com/0/254604/566c2ead-c173-7490-05fa-61326d26caf8.png)

![image.png](https://qiita-image-store.s3.amazonaws.com/0/254604/7a08fa93-61f0-f95f-69e9-bd6b8266d971.png)