Created
September 21, 2017 17:03
-
-
Save eickenberg/926850ddfa2a5800048efa8ae02f2ae9 to your computer and use it in GitHub Desktop.
cca and regularized kernel cca
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
def kcca_dual_coef(K1, K2, gamma, n_coefs=None, reduced_problem=False): | |
eigen_shape = tuple(np.array(K1.shape) * 2) | |
eigen_matrix1 = np.zeros(eigen_shape) | |
em_view = eigen_matrix1.view() | |
em_view.shape = 2, eigen_shape[0] // 2, 2, eigen_shape[1] // 2 | |
em_view = em_view.transpose(0, 2, 1, 3) | |
eigen_matrix2 = np.zeros(eigen_shape) | |
em_view2 = eigen_matrix2.view() | |
em_view2.shape = 2, eigen_shape[0] // 2, 2, eigen_shape[1] // 2 | |
em_view2 = em_view2.transpose(0, 2, 1, 3) | |
eye = np.eye(eigen_shape[0] // 2) | |
if n_coefs is None: | |
n_coefs = K1.shape[0] | |
eigvals = (eigen_shape[0] - n_coefs, eigen_shape[0] - 1) | |
if reduced_problem: | |
em_view[0, 1] = K2 | |
em_view[1, 0] = K1 | |
em_view2[0, 0] = K1 + gamma * eye | |
em_view2[1, 1] = K2 + gamma * eye | |
eigen_matrix1 += eigen_matrix1.T | |
eigen_matrix2 += eigen_matrix2.T | |
return linalg.eig(eigen_matrix1, eigen_matrix2, eigvals=eigvals) | |
else: | |
K1K2 = K1.dot(K2) | |
em_view[0, 1] = K1K2 | |
em_view[1, 0] = K1K2.T | |
em_view2[0, 0] = K1.dot(K1 + gamma * eye) # becomes K1 ** 2 + gamma * K1 | |
em_view2[1, 1] = K2.dot(K2 + gamma * eye) # becomes K2 ** 2 + gamma * K2 | |
eigen_matrix1 += eigen_matrix1.T | |
eigen_matrix2 += eigen_matrix2.T | |
corr, components = linalg.eigh(eigen_matrix1, eigen_matrix2, eigvals=eigvals) | |
return corr, components | |
def cca_primal_coef(X1, X2, gamma=0., center_data=True, n_coefs=None): | |
if center_data == True: | |
center1, center2 = X1.mean(axis=0), X2.mean(axis=0) | |
elif center_data == False: | |
center1, center2 = 0, 0 | |
else: | |
center1, center2 = center_data | |
X1 = X1 - center1 | |
X2 = X2 - center2 | |
S11 = X1.T.dot(X1) | |
S22 = X2.T.dot(X2) | |
S12 = X1.T.dot(X2) | |
S21 = S12.T | |
eigen_shape = (X1.shape[1] + X2.shape[1],) * 2 | |
eigen_matrix1 = np.zeros(eigen_shape) | |
em_view = eigen_matrix1.view() | |
em_view.shape = 2, eigen_shape[0] // 2, 2, eigen_shape[1] // 2 | |
em_view = em_view.transpose(0, 2, 1, 3) | |
eigen_matrix2 = np.zeros(eigen_shape) | |
em_view2 = eigen_matrix2.view() | |
em_view2.shape = 2, eigen_shape[0] // 2, 2, eigen_shape[1] // 2 | |
em_view2 = em_view2.transpose(0, 2, 1, 3) | |
eye1 = np.eye(X1.shape[1]) | |
eye2 = np.eye(X2.shape[1]) | |
em_view[0, 1] = S12 | |
em_view[1, 0] = S21 | |
em_view2[0, 0] = S11 + gamma * eye1 | |
em_view2[1, 1] = S22 + gamma * eye2 | |
half = X1.shape[1] | |
eigen_matrix1 += eigen_matrix1.T | |
eigen_matrix2 += eigen_matrix2.T | |
if n_coefs is None: | |
n_coefs = min(X1.shape[1], X2.shape[1], X1.shape[0], X2.shape[0]) | |
eigvals = (eigen_matrix1.shape[1] - n_coefs, eigen_matrix1.shape[1] - 1) | |
corr, components = linalg.eigh(eigen_matrix1, eigen_matrix2, eigvals=eigvals) | |
half = X1.shape[1] | |
components = components[:, ::-1] | |
corr = corr[::-1] | |
return corr, components, components[:half], components[half:] | |
def cca(X1, X2, gamma=0., center_data=True, n_coefs=None, use_dual='auto'): | |
n_samples, n_features1 = X1.shape | |
n_samples, n_features2 = X2.shape | |
if use_dual == 'auto': | |
if n_features1 + n_features2 > 2 * n_samples: | |
use_dual = True | |
else: | |
use_dual = False | |
if not use_dual: | |
return cca_primal_coef(X1, X2, gamma, center_data, n_coefs=n_coefs) | |
else: | |
# use linear kernel | |
if center_data is True: | |
center1, center2 = X1.mean(axis=0), X2.mean(axis=0) | |
elif center_data is False: | |
center1, center2 = 0., 0. | |
X1 = X1 - center1 | |
X2 = X2 - center2 | |
K1 = X1.dot(X1.T) | |
K2 = X2.dot(X2.T) | |
c, C = kcca_dual_coef(K1, K2, gamma, n_coefs=n_coefs) | |
half = C.shape[0] // 2 | |
primal1 = X1.T.dot(C[:half]) | |
primal2 = X2.T.dot(C[half:]) | |
return c, C, primal1, primal2 | |
if __name__ == "__main__": | |
X1 = np.random.randn(100, 10) | |
p = np.random.permutation(X1.shape[1]) | |
X2 = X1[:, p] + 1. * np.random.randn(*X1.shape) | |
X1_ = np.random.randn(50, 10) | |
X2_ = X1_[:, p] | |
r = cca(X1, X2, gamma=10000., use_dual=False, center_data=False, n_coefs=10) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment