From b9a51a2ab85459be3e7abebd8648802b296a1f67 Mon Sep 17 00:00:00 2001
From: James Hensman <james.hensman@gmail.com>
Date: Fri, 24 Jan 2014 10:19:03 +0000
Subject: [PATCH] added PCA to linalg

---
 GPy/util/linalg.py | 24 ++++++++++++++++++++++++
 1 file changed, 24 insertions(+)

diff --git a/GPy/util/linalg.py b/GPy/util/linalg.py
index 4cafe044..4cc2d7e3 100644
--- a/GPy/util/linalg.py
+++ b/GPy/util/linalg.py
@@ -547,3 +547,27 @@ def backsub_both_sides(L, X, transpose='left'):
     else:
         tmp, _ = lapack.dtrtrs(L, np.asfortranarray(X), lower=1, trans=0)
         return lapack.dtrtrs(L, np.asfortranarray(tmp.T), lower=1, trans=0)[0].T
+
+def PCA(Y, input_dim):
+    """
+Principal component analysis: maximum likelihood solution by SVD
+
+:param Y: NxD np.array of data
+:param input_dim: int, dimension of projection
+
+
+:rval X: - Nxinput_dim np.array of dimensionality reduced data
+:rval W: - input_dimxD mapping from X to Y
+
+"""
+    if not np.allclose(Y.mean(axis=0), 0.0):
+        print "Y is not zero mean, centering it locally (GPy.util.linalg.PCA)"
+
+        # Y -= Y.mean(axis=0)
+
+    Z = linalg.svd(Y - Y.mean(axis=0), full_matrices=False)
+    [X, W] = [Z[0][:, 0:input_dim], np.dot(np.diag(Z[1]), Z[2]).T[:, 0:input_dim]]
+    v = X.std(axis=0)
+    X /= v;
+    W *= v;
+    return X, W.T