通过三种方法计算CCA

Question

我最近研究了CCA的概念，想在MATLAB中实现它。然而，有一个现有的 matlab 命令 canoncorr 存在。我想编写自己的代码。我进行了广泛的研究并发现了三种方法：

1：Hardoon： 该方法使用拉格朗日乘子将问题分解为广义特征值问题。代码可以在这里找到：cca_hardoon 为了理智起见，我还在这里给出了代码：数据必须事先居中。

function [Wx, Wy, r] = cca(X,Y)

% CCA calculate canonical correlations
%
% [Wx Wy r] = cca(X,Y) where Wx and Wy contains the canonical correlation
% vectors as columns and r is a vector with corresponding canonical
% correlations.
%
% Update 31/01/05 added bug handling.

if (nargin ~= 2)
  disp('Inocorrect number of inputs');
  help cca;
  Wx = 0; Wy = 0; r = 0;
  return;
end


% calculating the covariance matrices
z = [X; Y];
C = cov(z.');
sx = size(X,1);
sy = size(Y,1);
Cxx = C(1:sx, 1:sx) + 10^(-8)*eye(sx);
Cxy = C(1:sx, sx+1:sx+sy);
Cyx = Cxy';
Cyy = C(sx+1:sx+sy,sx+1:sx+sy) + 10^(-8)*eye(sy);

%calculating the Wx cca matrix
Rx = chol(Cxx);
invRx = inv(Rx);
Z = invRx'*Cxy*(Cyy\Cyx)*invRx;
Z = 0.5*(Z' + Z);  % making sure that Z is a symmetric matrix
[Wx,r] = eig(Z);   % basis in h (X)
r = sqrt(real(r)); % as the original r we get is lamda^2
Wx = invRx * Wx;   % actual Wx values

% calculating Wy
Wy = (Cyy\Cyx) * Wx; 

% by dividing it by lamda
Wy = Wy./repmat(diag(r)',sy,1);

2。 MATLAB 方法 请注意，数据居中是在代码本身内完成的。

3.仅通过正态 SVD 进行 CCA： 此方法不需要 qr 分解，仅使用 svd 分解。我在这里引用了这篇文章：cca by svd。请参阅下面的文字文章，这些文章摘自所引用的文章。

我尝试自己编写这个程序，但没有成功。

function [A,B,r,U,V] = cca_by_svd(x,y)
% computing the means
N = size(x,1); mu_x = mean(x,1); mu_y = mean(y,1);
% substracting the means
x = x - repmat(mu_x,N,1); y = y - repmat(mu_y,N,1);

x = x.'; y = y.';
% computing the covariance matrices
Cxx = (1/N)*x*(x.');  Cyy = (1/N)*y*(y.'); Cxy = (1/N)*x*(y.');

%dimension
m = min(rank(x),rank(y));
%m = min(size(x,1),size(y,1));

% computing the quare root inverse of the matrix
[V,D]=eig(Cxx); d = diag(D);
% Making all the eigen values positive
d = (d+abs(d))/2; d2 = 1./sqrt(d); d2(d==0)=0; Cxx_iv=V*diag(d2)*inv(V);

% computing the quare root inverse of the matrix
[V,D]=eig(Cyy); d = diag(D);
% Making all the eigen values positive
d = (d+abs(d))/2; d2 = 1./sqrt(d); d2(d==0)=0; Cyy_iv=V*diag(d2)*inv(V);

Omega = Cxx_iv*Cxy*Cyy_iv;
[C,Sigma,D] = svd(Omega);
A = Cxx_iv*C; A = A(:,1:m);
B = Cyy_iv*D.'; B = B(:,1:m);
A = real(A); B = real(B);
U = A.'*x; V = B.'*y;
r = Sigma(1:m,1:m);

我正在运行此代码片段：

clc;clear all;close all;
load carbig;
X = [Displacement Horsepower Weight Acceleration MPG];
nans = sum(isnan(X),2) > 0;
x = X(~nans,1:3);
y = X(~nans,4:5);
[A1, B1, r1, U1, V1] = canoncorr(x, y);

[A2, B2, r2, U2, V2] = cca_by_svd(x, y);

[A3, B3, r3] = cca(x.',y.',1);

投影向量如下：

>> A1
A1 =
    0.0025    0.0048
    0.0202    0.0409
   -0.0000   -0.0027
>> A2
A2 =
    0.0025    0.0048
    0.0202    0.0410
   -0.0000   -0.0027
>> A3
A3 =
   -0.0302   -0.0050   -0.0022
    0.0385   -0.0420   -0.0176
    0.0020    0.0027   -0.0001
>> B1
B1 =
   -0.1666   -0.3637
   -0.0916    0.1078
>> B2
B2 =
   -0.1668   -0.3642
   -0.0917    0.1079
>> B3
B3 =
   0.0000 + 0.0000i   0.3460 + 0.0000i   0.1336 + 0.0000i
   0.0000 + 0.0000i  -0.0967 + 0.0000i   0.0989 + 0.0000i

问题：有人可以告诉我哪里出了问题吗？我提到的三种方法都解决相同的问题，理想情况下它们的解决方案应该收敛。我承认我的代码'cca_by_svd'可能是错误的，但hardoon的代码和matlab的输出应该是相同的。请指出我哪里出错了。 编辑 我已经重新检查并更正了我的代码。现在对于这个数据集，方法 2 和 3 收敛。

Answer 1

有一些事情是

cca(X,Y)

做不到的，而

canoncorr

却能做到：

一是对数据进行标准化。如果将

X = normc(X')'

（也适用于

Y

）添加到

cca(X,Y)

函数中，则输出

r

将与

canoncorr

的输出相匹配。如果您查看

canoncorr

的代码，您会发现它是从

X

和

Y

的 QR 分解开始的。

另一个区别是

eig

按升序对特征值进行排序，因此

cca(X,Y)

应该翻转

eig(Z)

的输出。

注意：尽管纠正了这些差异，但我无法完全恢复 Wx 和 Wy 以匹配

canoncorr

的输出。理想情况下，Wx'*Wx 在

cca

和

canoncorr

之间看起来应该完全相同。

Answer 2

使用 matlab 的默认 canoncorr 函数能得到相同的结果吗？