Commits

Jure Žbontar committed ea0b643

cca

  • Participants
  • Parent commits be90a2c

Comments (0)

Files changed (1)

 sys.argv.append('r')
 np.random.seed(42)
 
+
+
+samples = 100
+n = 101
+m = 35
+nm = n + m
+
+cov = np.random.random((nm, nm))
+cov = (cov + cov.T) / 2 + np.eye(nm)
+
+data = np.random.multivariate_normal(np.zeros(nm), cov, samples)
+data -= np.mean(data, axis=0)
+
+X = data[:,:n]
+Y = data[:,n:]
+
+#np.savetxt('data/X', X)
+#np.savetxt('data/Y', Y)
+
+Sxx = X.T.dot(X) / samples
+Sxy = X.T.dot(Y) / samples
+Syx = Y.T.dot(X) / samples
+Syy = Y.T.dot(Y) / samples
+
+#Sxx = np.array([[1, .4], [.4, 1]])
+#Sxy = np.array([[.5, .6], [.3, .4]])
+#Syx = Sxy.T
+#Syy = np.array([[1, .2], [.2, 1]])
+
 if sys.argv[1] in {'1', '2', '3', '4'}:
 
-    samples = 100
-    n = 1000
-    m = 35
-    nm = n + m
-
-    cov = np.random.random((nm, nm))
-    cov = (cov + cov.T) / 2 + np.eye(nm)
-
-    data = np.random.multivariate_normal(np.zeros(nm), cov, samples)
-    data -= np.mean(data, axis=0)
-
-    X = data[:,:n]
-    Y = data[:,n:]
-
-    #np.savetxt('data/X', X)
-    #np.savetxt('data/Y', Y)
-
-
-    Sxx = X.T.dot(X) / samples
-    Sxy = X.T.dot(Y) / samples
-    Syx = Y.T.dot(X) / samples
-    Syy = Y.T.dot(Y) / samples
-
-    #Sxx = np.array([[1, .4], [.4, 1]])
-    #Sxy = np.array([[.5, .6], [.3, .4]])
-    #Syx = Sxy.T
-    #Syy = np.array([[1, .2], [.2, 1]])
-
     Sxx_ = np.linalg.pinv(Sxx + 1e-5 * np.eye(Sxx.shape[0]))
     Syy_ = np.linalg.pinv(Syy + 1e-5 * np.eye(Syy.shape[0]))
 
+print('go')
+
 if sys.argv[1] == '1':
     Sxx_2 = scipy.linalg.sqrtm(Sxx_).real
 
     B = Ryy_.dot(U)
 
 if sys.argv[1] == 'r':
-    def center(x, xcenter):
-        if isinstance(xcenter, bool):
-            if xcenter:
-                xcenter = np.mean(x, axis=0)
-            else:
-                xcenter = np.zeros(ncx)
-        x -= xcenter
-        return xcenter
+    def cancor(x, y, xcenter=True, ycenter=True):
+        def center(x, xcenter):
+            if isinstance(xcenter, bool):
+                if xcenter:
+                    xcenter = np.mean(x, axis=0)
+                else:
+                    xcenter = np.zeros(ncx)
+            x -= xcenter
+            return xcenter
 
+        assert x.shape[0] == y.shape[0]
+        nr = x.shape[0]
 
-    x = np.array([[1,1.1,0.9], [-1,-1.2,-1.2], [2,2.3,2.4], [-2,-2.4,-2.3], [3,3,3], [4,4,4]])
-    y = np.array([[2,2.1],[-1,-1.2],[1,1.3],[-2,-2.4], [3,3], [4,4]])
-    xcenter = True
-    ycenter = True
+        ncx = x.shape[1]
+        ncy = y.shape[1]
 
-    assert x.shape[0] == y.shape[0]
-    nr = x.shape[0]
+        assert nr != 0 and ncx != 0 and ncy != 0
 
-    ncx = x.shape[1]
-    ncy = y.shape[1]
+        xcenter = center(x, xcenter)
+        ycenter = center(y, ycenter)
 
-    assert nr != 0 and ncx != 0 and ncy != 0
+        qx, rx = np.linalg.qr(x)
+        qy, ry = np.linalg.qr(y)
 
-    xcenter = center(x, xcenter)
-    ycenter = center(y, ycenter)
+        dx = np.linalg.matrix_rank(qx)
+        dy = np.linalg.matrix_rank(qy)
 
-    qx = np.linalg.qr(x)
-    qy = np.linalg.qr(y)
+        print(qx.shape)
 
-    dx = np.linalg.matrix_rank(qx)
-    dy = np.linalg.matrix_rank(qy)
+        assert dx > 0
+        assert dy > 0
 
-    assert dx > 0
-    assert dy > 0
+        u, s, v = np.linalg.svd(qx.T.dot(qy[:,:dy])[:dx])
+        xcoef = scipy.linalg.solve_triangular(rx[:dx,:dx], u)
+        ycoef = scipy.linalg.solve_triangular(ry[:dy,:dy], v.T)
 
-    d = np.zeros((nr, dy))
+        return xcoef, ycoef
 
+    A, B = cancor(X, Y)
 
-
-
+    print(X.shape, A.shape)
+    print(Y.shape, B.shape)
 
 if sys.argv[1] == 'sk':
     import sklearn.cross_decomposition