cca and regularized kernel cca

kcca.py

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)