随着对CCA的深入研究,是时候对CCA进行一下总结了。
本菜鸡主要研究方向为故障诊断,故会带着从应用角度进行理解。
![典型相关分析(Canonical Correlation Analysis,CCA)原理及Python、MATLAB实现](https://img-blog.csdnimg.cn/d8d94472576048039f19bde2ca53791f.gif)
典型相关分析
基本原理
从字面意义上理解CCA,我们可以知道,简单说来就是对不同变量之间做相关分析。较为专业的说就是,一种度量两组变量之间相关程度的多元统计方法。
关于相似性度量距离问题,在这里有一篇Blog可以参考参考。
首先,从基本的入手。
当我们需要对两个变量进行相关关系分析时,则常常会用到相关系数来反映。学过概率统计的小伙伴应该都知道的吧。还是解释一下。
相关系数:是一种用以反映变量之间相关关系密切程度的统计指标。相关系数是按积差方法计算,同样以两变量与各自平均值的离差为基础,通过两个离差相乘来反映两变量之间相关程度;着重研究线性的单相关系数。
其中,表示
的协方差矩阵,
为
的方差,
为
的方差.
复习了一下大学本科概率统计知识,那么,如果我们需要分析的对象是两组或者多组向量,又该怎么做呢?
CCA的数学表达
这里举例两组变量,那么我们的公式会是这样:
我们会得到一个这样的矩阵:
这样的话,我们把每个变量的相关系数都求了出来,不知道会不会和我一样觉得这样很繁琐呢。如果我们能找到两组变量之间的各自的线性组合,那么我们就只分析讨论线性组合之间的相关分析。
典型相关系数:是先对原来各组变量进行主成分分析,得到新的线性关系的综合指标,再通过综合指标之间的线性相关系数来研究原各组变量间相关关系。
现在我们利用主成分分析(PCA)的思想,可以把多个变量与多个变量之间的相关转化成两个变量之间的相关。
先得到两组变量的协方差矩阵
其中,.
把上面两组变量分别组合成两个变量U、V,则用线性表示
然后,找出最大可能的相关系数,
使得,,这样,就得到了典型相关系数;而其中的
为典型相关变量。
典型相关分析最朴素的思想:首先分别在每组变量中找出第一对典型变量,使其具有最大相关性,然后在每组变量中找出第二对典型变量,使其分别与本组内的第一对典型变量不相关,第二对本身具有次大的相关性。如此下去,直到进行到K步,两组变量的相关系被提取完为止,可以得到K组变量。
So,
注意:此时的若不能反映两组变量之间的相关关系,我们需要继续构造下一组关系变量来表示,具体可构造
对这样的关系
直到为止
其中,我们需要一个约束条件满足,使得
典型相关系数公式
在此约束条件下,系数得到最大,则使得
最大。
典型相关系数及变量的求法
下面一起来求(这里只例举第一典型相关系数)
(一起来复习高数–拉格朗日乘数法)
前提条件,我们有个计算公式,约束条件也有了,故这是一个求解条件极值题呀!!!
列出我们的拉格朗日函数:
其中,表示拉格朗日乘子参数。
由上述典型相关系数公式 可知,我们只需求其系数
的最大值,即可
对做一阶偏导运算:
也就是
将上式分别左乘得,
由于约束条件可知,解得,
此时,我们来对比一下上面列出的求解公式,有没有,是不是,一模一样,我们的拉格朗日乘子
,也就是说,
即为我们要求解的典型相关系数。
我们由式,可得
由于,再将此代入
中,可得
上面的式子是不是很熟悉,
故,的特征值为
,特征向量为
。
到此,我们可求得相应的及
。
典型相关分析应用
基于 CCA 的故障检测方法
对于CCA应用在故障检测中,基于 CCA 的故障检测方法可以视为基于 PCA 和基于 PLS 故障检测方法的一种扩展。
基本思想:是利用典型相关关系构建一个残差发生器, 通过对残差信号的评价做出故障检测的相应决策。该方法中提出了 4 个统计量, 将输入空间分为两个部分, 即与输出空间相关的子空间和与输出空间不相关的子空间;同理,将输出空间分为两个部分, 即与输入空间相关的子空间和与输入空间不相关的子空间。
设 和
表示测量的过程输入和输出向量,
和
分别表示相应的数据维数。对两个向量进行去均值, 可得
式中: 和
分别表示输入变量均值和输出变量均值。假设采样得到 N 个过程数据, 并组成如下输入输出数据集
式中:是按式(1)(2)中心化后的输入输出向量, 相应的平均值
并且, 协方差矩阵 和输入输出的互协方差矩阵
分别为:
结合 CCA 方法, 可得:
式中: κ 为主元个数, 为典型相关系数值。
令,
由 CCA 方法可知,和
具有密切的相关性。
但是在实际系统中, 测量变量难免受到噪声影响, 两者之间的相关性可表示为:
式中: 为噪声项, 并且与
弱相关。基于此, 残差向量
其中的
为系数矩阵,上面介绍了CCA的数学运算,这里应该能理解,只不过多加了一个噪声向量。
假设输入输出数据服从高斯分布。已知线性变换不改变随机变量的分布, 所以残差信号也服
从高斯分布, 其协方差矩阵:
同理, 还可以得到另一残差向量
其协方差矩阵
由式(9)(11) 可以看出, 残差 r1和 r2的协方差相同。 对于故障检测, 可构造如下两个统计量:
注意到统计量 用于检测发生在输出子空间且输入相关的那部分故障, 为了检测与输入不相关的那部分故障, 可构造一个统计量
同理, 为了检测发生在输入空间且与输出不相关的那部分故障, 可构造另一统计量
由以上分析可知, 通过确定主元个数 κ, 可以得到4 个统计量、
、
和
进行故障检测。
关于过程故障监控的统计量
,在深度学习、机器学习、故障诊断领域用的较多,这里可参考
的相关内容。
应用部分参考自一篇Paper[1]. CHEN Zhiw en,DING S X,ZHANG Kai,et al.Canonical correlation analysis- based fault detection methods with application to alumina evaporation process[J].Control Engineering Practice,2016,46:51- 58.
Python代码:
## 通过sklearn工具包内置的CCA实现
import numpy as np
from sklearn.cross_decomposition import CCA
from icecream import ic # ic用于显示,类似于print
A = [[3, 4, 5, 6, 7] for i in range(2000)]
B = [[8, 9, 10, 11, 12] for i in range(2000)]
# 注意在A、B中的数为输入变量及输出变量参数
# 建模
cca = CCA(n_components=1) # 若想计算第二主成分对应的相关系数,则令cca = CCA(n_components=2)
# 训练数据
cca.fit(X, Y)
# 降维操作
X_train_r, Y_train_r = cca.transform(X, Y)
#输出相关系数
ic(np.corrcoef(X_train_r[:, 0], Y_train_r[:, 0])[0, 1]) #如果想计算第二主成分对应的相关系数 print(np.corrcoef(X_train_r[:, 1], Y_train_r[:, 1])[0, 1])
另有一个包含可视化CCA的Python代码在 这里。
Matlab代码:
function[ccaEigvector1, ccaEigvector2] = CCA(data1, data2)
dataLen1 = size(data1, 2);
dataLen2 = size(data2, 2);
% Construct the scatter of each view and the scatter between them
data = [data1 data2];
covariance = cov(data);
% Sxx = covariance(1 : dataLen1, 1 : dataLen1) + eye(dataLen1) * 10^(-7);
Sxx = covariance(1 : dataLen1, 1 : dataLen1);
% Syy = covariance(dataLen1 + 1 : size(covariance, 2), dataLen1 + 1 : size(covariance, 2)) ...
% + eye(dataLen2) * 10^(-7);
Syy = covariance(dataLen1 + 1 : size(covariance, 2), dataLen1 + 1 : size(covariance, 2));
Sxy = covariance(1 : dataLen1, dataLen1 + 1 : size(covariance, 2));
% Syx = Sxy';
% using SVD to compute the projection
Hx = (Sxx)^(-1/2);
Hy = (Syy)^(-1/2);
H = Hx * Sxy * Hy;
[U, D, V] = svd(H, 'econ');
ccaEigvector1 = Hx * U;
ccaEigvector2 = Hy * V;
% make the canonical correlation variable has unit variance
ccaEigvector1 = ccaEigvector1 * diag(diag((eye(size(ccaEigvector1, 2)) ./ sqrt(ccaEigvector1' * Sxx * ccaEigvector1))));
ccaEigvector2 = ccaEigvector2 * diag(diag((eye(size(ccaEigvector2, 2)) ./ sqrt(ccaEigvector2' * Syy * ccaEigvector2))));
end
坚持读Paper,坚持做笔记!!!
To Be No.1
哈哈哈哈过路能❤关注、收藏、点个赞❤三连就最好不过了
ღ( ´・ᴗ・` )
❤
『
对自己的爱好保持热情,不要太功利!
』
文章出处登录后可见!