目录
前言
典型相关分析（Canonical Correlation Analysis，CCA）是一种用于探索两组变量之间相关性的多元统计方法。它的核心目标是找到两组变量的线性组合，使得这些组合之间的相关性最大化，从而揭示两组变量之间的潜在联系。在BCI脑机接口领域，SSVEP稳态视觉诱发电位中，经常使用的一种方法。
理论学习
关于cov协方差矩阵，补充下图，可以看到组成矩阵的各个元素，清楚各个位置上元素的求解。
随后回到CCA，现在研究x,y，两对线性组合之间的关系。
此处将（8.6）带入（8.5）相关系数p公式中，得到式(8.7)。 其中用到了如下的运算性质
得到了假设情况下，相关系数p的表达式。要求解该假设下，p的最大值。使用拉格朗日乘数法。
这里矩阵求导，补充说明，使用了以下两个矩阵求导公式。
因为得到了 μ = λ = p 的结论，所以将（8.9)中的 μ 换为 λ 。
总结思路
所以如下图可知，知晓理论求解的步骤后，只要用Matlab求 x, y 的各种cov协方差，配合inv求逆矩阵，即可得到M1和M2矩阵。然后用matlab求出M1和M2的特征值与对应的特征向量。然后找特征值最大的特征向量a,b即可。（因为 μ=λ=p）特征值最大时，相关系数p也是最大的。
由此，实现上面图片的运行逻辑，得到以下matlab代码。文件命名为CCA.m
function [A,B]=CCA(X,Y)
%CCA canonical correlation analysis
% [A,B]=cca(X,Y)
dimx=size(X,2); dimy=size(Y,2); %列数(变量数)
dim=min(dimx,dimy);
Sall = cov([X Y]); %cov求协方差矩阵。原阵X大小为M*N，则cov(X)大小为N*N的矩阵
Sxx=Sall(1:dim,1:dim);
Syy=Sall(dim+1:end,dim+1:end);
Sxy=Sall(1:dim,dim+1:end);
Syx=Sall(dim+1:end,1:dim);
iSxx=inv(Sxx); %矩阵求逆
iSyy=inv(Syy);
M1 = iSxx*Sxy*iSyy*Syx;
M2 = iSyy*Syx*iSxx*Sxy;
[vec1,val1]=eig(M1); %val对角阵，元素是特征值。 vec各列为对应特征值的特征向量
[vec2,val2]=eig(M2);
% A......找到每列的最大值和对应的行索引
[maxVals1, rowIndices1] = max(val1);
% 找到最大值的列标
[~, colIndex1] = max(maxVals1);
A = vec1(:, colIndex1);
% B......找到每列的最大值和对应的行索引
[maxVals2, rowIndices2] = max(val2);
% 找到最大值的列标
[~, colIndex2] = max(maxVals2);
B = vec2(:, colIndex2);matlab实践
测试用数据，从书本198页，8.3 R语言上机实现中，例8.1拿的数据。
课本表明其计算结果如下
所以把课本答案摘出来，如果算法可行，应得到如下A,B矩阵为结果
A=[0.0314; B=[0.0661;
-0.4932; 0.0168;
0.0082]; -0.0140];测试代码test.m
测试用代码, 使用书本数据，做计算验证。 文件名 test.m
baseA = [191, 36, 50;
193, 38, 58;
189, 35, 46;
211, 38, 56;
176, 31, 74;
169, 34, 50;
154, 34, 64;
193, 36, 46;
176, 37, 54;
156, 33, 54;
189, 37, 52;
162, 35, 62;
182, 36, 56;
167, 34, 60;
154, 33, 56;
166, 33, 52;
247, 46, 50;
202, 37, 62;
157, 32, 52;
138, 33, 68];
baseB = [5, 162, 60;
12, 101, 101;
13, 155, 58;
8, 101, 38;
15, 200, 40;
17, 120, 38;
14, 210, 105;
6, 70, 31;
4, 60, 25;
15, 225, 73;
2, 110, 60;
12, 105, 37;
4, 101, 42;
6, 125, 40;
17, 251, 250;
13, 210, 115;
1, 50, 50;
12, 210, 120;
11, 230, 80;
2, 110, 43];
%% cca
x = baseA;
y = baseB;
[A, B] = CCA(x, y);
A
B将CCA.m和test.m放在同一文件夹内，运行tset.m得到如下运行结果
发现与课本结果，只需乘一个常数倍数即可得到。乘任意常数不会改变相关关系，所以求解应该是对的。这个倍数可能是不同语言函数求解导致，一些函数运算会自动进行规范化。
结果规范化
我们把A,B结果进行规范化处理，修改如下
function [A,B]=CCA(X,Y)
%CCA canonical correlation analysis
% [A,B]=cca(X,Y)
dimx=size(X,2); dimy=size(Y,2); %列数(变量数)
dim=min(dimx,dimy);
Sall = cov([X Y]); %cov求协方差矩阵。原阵X大小为M*N，则cov(X)大小为N*N的矩阵
Sxx=Sall(1:dim,1:dim);
Syy=Sall(dim+1:end,dim+1:end);
Sxy=Sall(1:dim,dim+1:end);
Syx=Sall(dim+1:end,1:dim);
iSxx=inv(Sxx); %矩阵求逆
iSyy=inv(Syy);
M1 = iSxx*Sxy*iSyy*Syx;
M2 = iSyy*Syx*iSxx*Sxy;
[vec1,val1]=eig(M1); %val对角阵，元素是特征值。 vec各列为对应特征值的特征向量
[vec2,val2]=eig(M2);
% A......找到每列的最大值和对应的行索引
[maxVals1, rowIndices1] = max(val1);
% 找到最大值的列标
[~, colIndex1] = max(maxVals1);
A = vec1(:, colIndex1);
A = A* diag(diag((eye(size(A, 2))./ sqrt(A’ * Sxx * A))))* -1:
% B......找到每列的最大值和对应的行索引
[maxVals2, rowIndices2] = max(val2);
% 找到最大值的列标
[~, colIndex2] = max(maxVals2);
B = vec2(:, colIndex2);
B = B * diag(diag((eye(size(B, 2))./ sqrt(B’ * Syy * B))))* -1:运行结果如下图
课本结果 A=[0.0314; -0.4932; 0.0082]; B=[0.0661; 0.0168; -0.0140];
其运行结果与书本解答，仅有千分位的精度差别。鉴于课本使用R语言求解，我们是matlab，不同语言软件有差异，会有些许影响结果精度，所以该代码应该是准确的。
改用SVD求解特征值特征向量
将求解步骤钟，解特征值和特征向量改为使用svd分解，即可得到上一篇文章的CCA.m
%******** CCA.m *********%
function [A,B,C1]=CCA(X,Y)
%CCA canonical correlation analysis
% [A,B]=cca(X,Y)
dimx=size(X,2); dimy=size(Y,2); %列数(变量数)
Sall = cov([X Y]); %cov求协方差矩阵。原阵X大小为M*N，则cov(X)大小为N*N的矩阵
dim=min(dimx,dimy);
Sxx=Sall(1:dim,1:dim);
Syy=Sall(dim+1:end,dim+1:end);
Sxy=Sall(1:dim,dim+1:end);
Syx=Sall(dim+1:end,1:dim);
iSxx=inv(Sxx); %矩阵求逆
iSyy=inv(Syy);
M1 = iSxx*Sxy*iSyy*Syx;
M2 = iSyy*Syx*iSxx*Sxy;
% 使用 SVD 分解
[U1, S1, V1] = svd(M1); % M1 的奇异值分解
[U2, S2, V2] = svd(M2); % M2 的奇异值分解
% 手动单位化 U1 和 U2：使得每列（奇异向量）的模长为 1
for i = 1:size(U1, 2)
U1(:, i) = U1(:, i) / sqrt(U1(:, i)' * U1(:, i)); % 每列单位化
end
for i = 1:size(U2, 2)
U2(:, i) = U2(:, i) / sqrt(U2(:, i)' * U2(:, i)); % 每列单位化
end
% 手动归一化 S1 和 S2，使得最大奇异值为 1
S1 = S1 / max(diag(S1)); % 将 S1 的最大奇异值归一化为 1
S2 = S2 / max(diag(S2)); % 将 S2 的最大奇异值归一化为 1
% A ...... 选择最大奇异值对应的奇异向量
[maxVals1, colIndex1] = max(diag(S1)); % 找到最大奇异值所在的列
A = U1(:, colIndex1); % 选择对应的奇异向量
C1 = maxVals1; % 最大奇异值
A = A * diag(1 ./ sqrt(A' * Sxx * A)); % 归一化，注意矩阵的大小匹配
% B ...... 选择最大奇异值对应的奇异向量
[maxVals2, colIndex2] = max(diag(S2)); % 找到最大奇异值所在的列
B = U2(:, colIndex2); % 选择对应的奇异向量
B = B * diag(1 ./ sqrt(B' * Syy * B)); % 归一化，注意矩阵的大小匹配
参考书籍
《多元统计分析（第5版）》何晓群