# Principal Component Analysis in Tensorflow

To learn the low-level API of Tensorflow I am trying to implement some traditional machine learning algorithms. The following Python script implements Principal Component Analysis using gradient descent. It finds component weights that maximize the variance of each component. The weights are constrained to be orthonormal, as required by the PCA definition.

The results are consistent with Scikit-Learn's PCA implementation, so I assume the code works correctly.

Since these are my first baby steps with Tensorflow, I doubt that I used the API in an idiomatic way. How can my use of Tensorflow be improved?

import numpy as np
import tensorflow as tf
import matplotlib.pyplot as plt

x_data = np.random.randn(100, 3) @ [[0, 0, 1], [1, -2, 0], [0, 1, 0]]  # data set
n_components = 3  # number of desired components

def proj(v, u):
'''project vector v on vector u'''
return u * (tf.reduce_sum(u * v) / tf.reduce_sum(u**2))

def orthonormalize(w):
'''orthonormalize vectors basis w'''
n = w.shape
ortho_vecs = []
for i in range(n):
next_vector = w[:, i]
for vec in ortho_vecs:
next_vector -= proj(w[:, i], vec)
ortho_vecs.append(next_vector)
u = tf.stack(ortho_vecs, axis=-1)
return u / tf.reduce_sum(u**2, axis=-2)**0.5  # normalize basis vectors

# reset tensorflow (because I run this in a persisting Jupyter kernel)
tf.reset_default_graph()

# placeholder tensor for input data
x = tf.placeholder(tf.float32, shape=[None, x_data.shape[-1]])

# PCA is one densely connected layer with no biases and linear activation function.
# Each unit corresponds to one component, and components are constrained to be
# orthonormal.
y = tf.layers.dense(x, units=n_components, use_bias=False,
kernel_constraint=orthonormalize)

# mean and variance of the resulting components
mean_y = tf.reduce_mean(y, axis=-2)
var_y = tf.reduce_mean((y - mean_y)**2, axis=-2)

# minimize the sum of component variances
# This seems not only to converge to the correct components, but even sorts
# the components from largest to smallest. I don't know why... possibly has to
# do with the orthogonalization where the first component is unconstrained and
# each successive component is more and more constrained...
train = optimizer.minimize(-tf.reduce_sum(var_y))

with tf.Session() as sess:
# initialize weights
init = tf.global_variables_initializer()
sess.run(init)

# iterate training
train_var = []
train_weights = []
for _ in range(1000):
# stochastic gradient descent (need at least two samples for computing
# the variance... call it micro-batch?)
x_batch = x_data[np.random.randint(x_data.shape, size=2)]
#x_batch = x_data  # full data set

_, v, w = sess.run([train, var_y] + tf.trainable_variables(),
feed_dict={x: x_batch})

train_weights.append(w.ravel())

train_var.append(np.var(x_data @ w, axis=0))

# plot convergence of weights and variances
plt.subplot(2, 1, 1)
plt.plot(train_var)
plt.subplot(2, 1, 2)
plt.plot(train_weights)

writer = tf.summary.FileWriter("test", sess.graph)

# print final component weights
print(sess.run(tf.trainable_variables()))