SpectralNet (Moons 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.datasets import make_moons

n_samples = 1000

x_sample_train, labels_train = make_moons(n_samples = n_samples, noise = 0.05)
x_sample_test, labels_test = make_moons(n_samples = n_samples, noise = 0.05)

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

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

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

plt.show()
```
![image.png](https://qiita-image-store.s3.amazonaws.com/0/254604/b09117b5-654b-d964-5a60-dd3039aa8eb9.png)

##サンプルコード
```python
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 = 2
n_units_1 = 32
n_units_2 = 32
n_clusters = 2
n_nbrs = 5

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

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)
  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)
  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()
  
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/083bc016-f41b-fc7a-fedd-b599c1de7a7d.png)