CSP(Common spatial patterns)共空间模式算法简介

CSP算法基础知识

在了解CSP算法之前的相关知识基础

空间滤波器(spatial filters)

滤波是指通过操作接受或拒绝一定的频率。具体到空间滤波器，指的是通过空间滤波器对矩阵的元素进行修改，达到预期的效果。

空间滤波器由一下两部分组成：

一个邻域（通常为规模较小矩阵）
对该邻域的覆盖的矩阵元素执行的预定义操作

下图为使用3 × \times × 3的邻域的线性空间滤波举例，过程就是对原矩阵的元素逐个按照预选定义的滤波器的操作进行运算

对于矩阵中的任意一点 ( x , y ) (x,y) (x,y)，滤波器的响应 g ( x , y ) g(x,y) g(x,y)是滤波器系数和滤波器覆盖的矩阵系数乘积之和：
g ( x , y ) = ω ( − 1 , − 1 ) f ( x − 1 , y − 1 ) + ω ( − 1 , 0 ) f ( x − 1 , y ) + . . . + ω ( 1 , 1 ) f ( x + 1 , y + 1 ) g(x,y)=\omega(-1,-1)f(x-1,y-1)+\omega(-1,0)f(x-1,y)+...+\omega(1,1)f(x+1,y+1) g(x,y)=ω(−1,−1)f(x−1,y−1)+ω(−1,0)f(x−1,y)+...+ω(1,1)f(x+1,y+1)
将其推广至 m × n m\times n m×n的线性滤波器，设 m = 2 a + 1 m=2a+1 m=2a+1， n = 2 b + 1 n=2b+1 n=2b+1，则该滤波器响应函数 g ( x , y ) g(x,y) g(x,y)：
g ( x , y ) = ∑ s = − a a ∑ t = − b b ω ( s , t ) f ( x + s , y + t ) g(x,y)=\sum_{s=-a}^a\sum_{t=-b}^b\omega(s,t)f(x+s,y+t) g(x,y)=s=−a∑at=−b∑bω(s,t)f(x+s,y+t)

协方差与协方差矩阵

此处的协方差指的是样本的协方差（区别于随机变量的协方差）

样本的协方差

对于含有多个属性的样本，当分析两个属性之间的线性相关程度时，将每个多属性样本看做一个多维度向量，需要计算两个维度之间的协方差，协方差越大，两维度越相关。

设样本对应的多维随机变量为 X = [ X 1 , X 2 , . . . , X n ] T X=[X_1,X_2,...,X_n]^T X=[X1,X2,...,Xn]T，样本集合为 { X . j = [ x 1 j , x 2 j , . . . , x n j ] T ∣ 1 ≤ j ≤ m } \{X_{.j}=[x_{1j},x_{2j},...,x_{nj}]^T|1\le j\le m\} {X.j=[x1j,x2j,...,xnj]T∣1≤j≤m}，m为样本数量。对于变量中a属性与b属性的线性相关性，其中 1 ≤ a ≤ n 1\le a\le n 1≤a≤n， 1 ≤ b ≤ n 1\le b\le n 1≤b≤n，n为样本维度，其协方差公式为
q a b = ∑ j = 1 m ( x a j − x a ˉ ) ( x b j − x b ˉ ) m − 1 q_{ab}=\frac{\sum_{j=1}^m(x_{aj}-\bar{x_a})(x_{bj}-\bar{x_b})}{m-1} qab=m−1∑j=1m(xaj−xaˉ)(xbj−xbˉ)
其中 m − 1 m-1 m−1代表随机变量，自由度减1。

协方差矩阵

对于多属性变量，我们需要计算各个属性间的线性相关性。因此形成正方形矩阵

样本集合为 { X . j = [ x 1 j , x 2 j , . . . , x n j ] T ∣ 1 ≤ j ≤ m } \{X_{.j}=[x_{1j},x_{2j},...,x_{nj}]^T|1\le j\le m\} {X.j=[x1j,x2j,...,xnj]T∣1≤j≤m}可依据生成 n × n n\times n n×n的矩阵，将每一项按照协方差公式计算，使用 ∑ ^ \hat{\sum} ∑^ 表示样本的协方差矩阵
∑ ^ = [ q 11 q 12 ⋯ q 1 n q 21 q 22 ⋯ q 2 n ⋮ ⋮ ⋱ ⋮ q n 1 q n 2 ⋯ q n n ] = 1 m − 1 × [ ∑ j = 1 m ( x 1 j − x 1 ˉ ) ( x 1 j − x 1 ˉ ) ∑ j = 1 m ( x 1 j − x 1 ˉ ) ( x 2 j − x 2 ˉ ) ⋯ ∑ j = 1 m ( x 1 j − x 1 ˉ ) ( x n j − x n ˉ ) ∑ j = 1 m ( x 2 j − x 2 ˉ ) ( x 1 j − x 1 ˉ ) ∑ j = 1 m ( x 2 j − x 2 ˉ ) ( x 2 j − x 2 ˉ ) ⋯ ∑ j = 1 m ( x 2 j − x 2 ˉ ) ( x n j − x n ˉ ) ⋮ ⋮ ⋱ ⋮ ∑ j = 1 m ( x n j − x n ˉ ) ( x 1 j − x 1 ˉ ) ∑ j = 1 m ( x n j − x n ˉ ) ( x 2 j − x 2 ˉ ) ⋯ ∑ j = 1 m ( x n j − x n ˉ ) ( x n j − x n ˉ ) ] = 1 m − 1 ∑ j = 1 m ( x . j − x . ˉ ) ( x . j − x . ˉ ) T \hat{\sum}=\left[ \begin{matrix} q_{11} & q_{12} & \cdots & q_{1n} \\ q_{21} & q_{22} & \cdots & q_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ q_{n1} & q_{n2} & \cdots & q_{nn} \\ \end{matrix} \right] \\{}\\ =\frac1{m-1}\times \left[ \begin{matrix} {\sum_{j=1}^m(x_{1j}-\bar{x_1})(x_{1j}-\bar{x_1})} & {\sum_{j=1}^m(x_{1j}-\bar{x_1})(x_{2j}-\bar{x_2})} & \cdots & {\sum_{j=1}^m(x_{1j}-\bar{x_1})(x_{nj}-\bar{x_n})} \\{}\\ {\sum_{j=1}^m(x_{2j}-\bar{x_2})(x_{1j}-\bar{x_1})} & {\sum_{j=1}^m(x_{2j}-\bar{x_2})(x_{2j}-\bar{x_2})} & \cdots & {\sum_{j=1}^m(x_{2j}-\bar{x_2})(x_{nj}-\bar{x_n})} \\ \vdots & \vdots & \ddots & \vdots \\ {\sum_{j=1}^m(x_{nj}-\bar{x_n})(x_{1j}-\bar{x_1})} & {\sum_{j=1}^m(x_{nj}-\bar{x_n})(x_{2j}-\bar{x_2})} & \cdots & {\sum_{j=1}^m(x_{nj}-\bar{x_n})(x_{nj}-\bar{x_n})} \\ \end{matrix} \right] \\{}\\ =\frac1{m-1}\sum_{j=1}^m(x_{.j}-\bar{x_.})(x_{.j}-\bar{x_.})^T ∑^=⎣⎢⎢⎢⎡q11q21⋮qn1q12q22⋮qn2⋯⋯⋱⋯q1nq2n⋮qnn⎦⎥⎥⎥⎤=m−11×⎣⎢⎢⎢⎢⎢⎡∑j=1m(x1j−x1ˉ)(x1j−x1ˉ)∑j=1m(x2j−x2ˉ)(x1j−x1ˉ)⋮∑j=1m(xnj−xnˉ)(x1j−x1ˉ)∑j=1m(x1j−x1ˉ)(x2j−x2ˉ)∑j=1m(x2j−x2ˉ)(x2j−x2ˉ)⋮∑j=1m(xnj−xnˉ)(x2j−x2ˉ)⋯⋯⋱⋯∑j=1m(x1j−x1ˉ)(xnj−xnˉ)∑j=1m(x2j−x2ˉ)(xnj−xnˉ)⋮∑j=1m(xnj−xnˉ)(xnj−xnˉ)⎦⎥⎥⎥⎥⎥⎤=m−11j=1∑m(x.j−x.ˉ)(x.j−x.ˉ)T
在对样本进行归一化
y . j = x . j − x . ˉ z i . = y i . / σ i y_{.j}=x_{.j}-\bar{x_.} \\ z_{i.}=y_{i.}/\sigma_i y.j=x.j−x.ˉzi.=yi./σi
协方差矩阵 ∑ ^ \hat{\sum} ∑^可表示为
∑ ^ = 1 m − 1 ∑ j = 1 m z . j z . j T \hat{\sum}=\frac1{m-1}\sum_{j=1}^mz_{.j}z_{.j}^T ∑^=m−11j=1∑mz.jz.jT

方阵特征值与特征向量

A是n阶矩阵，如果数λ和n维非零列向量x使关系式
A x = λ x Ax=\lambda x Ax=λx
则称数λ为矩阵A的特征值，向量x为矩阵A的特征向量

矩阵对角化及同时对角化

矩阵对角化

若n阶矩阵 A A A与对角矩阵 Λ \Lambda Λ相似，即 P − 1 A P = Λ P^{-1}AP=\Lambda P−1AP=Λ
Λ = [ λ 1 λ 2 ⋱ λ n ] \Lambda=\left[ \begin{matrix} \lambda_1 & \quad & \quad & \quad\\ \quad & \lambda_2 & \quad & \quad\\ \quad & \quad & \ddots & \quad \\ \quad & \quad & \quad & \lambda_n \\ \end{matrix} \right] Λ=⎣⎢⎢⎡λ1λ2⋱λn⎦⎥⎥⎤
则 λ 1 , λ 2 , . . . , λ n \lambda_1,\lambda_2,...,\lambda_n λ1,λ2,...,λn为n个特征值， P P P中每个列向量 p i p_i pi是 A A A对应于特征值 λ i \lambda_i λi的特征向量。

该过程被称为矩阵对角化，也是特征值分解过程

PS：对于对称矩阵，其 P P P为正交矩阵，即 P − 1 = P T P^{-1}=P^T P−1=PT，故 A = P Λ P T A=P\Lambda P^T A=PΛPT

矩阵同时对角化

设 A , B A,B A,B是数域F上两个n阶矩阵,若存在n阶可逆矩阵P,使 P − 1 A P P^{-1}AP P−1AP, P − 1 B P P^{-1}BP P−1BP同时为对角矩阵,则称 A , B A,B A,B可同时相似对角化.

矩阵白化

给出任意矩阵，一般情况下，其协方差矩阵不是对角矩阵，矩阵白化是指用一个白化矩阵 A A A，使得 Y = A X Y=AX Y=AX的协方差矩阵转化为对角矩阵，省去过程推导，对于矩阵 X X X，其协方差矩阵:
c o v ( X ) = 1 m − 1 X T X cov(X)=\frac1{m-1}X^TX cov(X)=m−11XTX
对矩阵进行特征值分解:
c o v ( X ) = U λ U T cov(X)=U\lambda U^T cov(X)=UλUT
U和 λ \lambda λ分别表示cov(X)的特征向量和特征值对角矩阵
则白化矩阵为：
P = λ − 1 U T P=\sqrt{\lambda^{-1}}U^T P=λ−1 UT

CSP算法

共空间模式（CSP）是一种对两分类任务下的空域滤波特征提取算法，能够从多通道的脑机接口数据里面提取出每一类的空间分布成分。

算法简介

CSP算法旨在设计空间滤波器使得两组脑电时空信号矩阵滤波后，方差值差异最大化，从而得到具有较高区分度的特征向量。用于下一步将特征向量送入分类器进行分类。

上图中，脑电时空信号矩阵的维数 N × T N\times T N×T，N代表脑电通道数（导联数目），T代表时间层面的采样点个数。一个矩阵的单位为一个trail，表示一次测试。滤波器的维数是 2 m × N 2m\times N 2m×N，N同样代表导联数目，2m为生成该滤波器时，人为设定的特征选取个数。

上面说明了CSP算法的结果的用途，下面说一下其生成的过程。其利用的基本原理是协方差矩阵对角化。输入为训练集的全部trail矩阵，每个矩阵对应的任务标签（比如左手或右手运动想象），规定的特征数。

下面详细介绍一下CSP算法的步骤

构造空间滤波器

两类数据 X 1 X_1 X1， X 2 X_2 X2归一化后协方差矩阵 R 1 R_1 R1， R 2 R_2 R2分别为，其中 t r a c e ( X ) trace(X) trace(X)表示矩阵对角线上元素和，矩阵有两个下标，1，2表示标签，i表示第i个trail：
R 1 , i = X 1 , i X 1 , i T t r a c e ( X 1 , i X 1 , i T ) R 2 , i = X 2 , i X 2 , i T t r a c e ( X 2 , i X 2 , i T ) R_{1,i}=\frac{X_{1,i}X_{1,i}^T}{trace(X_{1,i}X_{1,i}^T)}\\{}\\ R_{2,i}=\frac{X_{2,i}X_{2,i}^T}{trace(X_{2,i}X_{2,i}^T)} R1,i=trace(X1,iX1,iT)X1,iX1,iTR2,i=trace(X2,iX2,iT)X2,iX2,iT

将所有矩阵的协方差矩阵求出，求各自标签的均值空间协方差矩阵，其中，K代表同一标签的时空矩阵trail数量
R c = ∑ i R c , i R c ˉ = 1 K R c ( c = 1 , 2 ) R_c=\sum_iR_{c,i}\\ \bar{R_c}=\frac1KR_c\quad(c=1,2) Rc=i∑Rc,iRcˉ=K1Rc(c=1,2)

之后，求混合空间协方差矩阵 R R R，此处 R R R为对角矩阵:
R = R 1 ˉ + R 2 ˉ R=\bar{R_1}+\bar{R_2} R=R1ˉ+R2ˉ

利用混合空间协方差矩阵，求白化特征值矩阵P，其中U是R的特征向量矩阵，λ是对应的特征值构成的对角阵
R = U λ U T P = λ − 1 U T R=U\lambda U^T\\ P=\sqrt{\lambda^{-1}}U^T R=UλUTP=λ−1 UT

构造空间滤波器

使用白化矩阵对 R 1 R_1 R1， R 2 R_2 R2进行如下变换，注意此处 R c R_c Rc为相同标签矩阵和：
S 1 = P R 1 P T , S 2 = P R 2 P T S 1 = B 1 λ 1 B 1 T , S 2 = B 2 λ 2 B 2 T λ 1 + λ 2 = I B 1 = B 2 = B S_1=PR_1P^T,\quad S_2=PR_2P^T\\ S_1=B_1\lambda_1B_1^T,\quad S_2=B_2\lambda_2B_2^T\\ \lambda_1+\lambda_2=I\\ B_1=B_2=B S1=PR1PT,S2=PR2PTS1=B1λ1B1T,S2=B2λ2B2Tλ1+λ2=IB1=B2=B
其中 λ 1 \lambda_1 λ1， λ 2 \lambda_2 λ2均为对应特征值对角矩阵，I为单位矩阵

从上面推导中得出，两类矩阵的特征值相加总为1， S 1 S_1 S1的最大特征值对应 S 2 S_2 S2的最小特征值，且两个矩阵的特征向量 B 1 B_1 B1， B 2 B_2 B2相等，两矩阵可同时对角化。

S 1 S_1 S1的最大特征值所对应的特征向量使 S 2 S_2 S2有最小的特征值，反之亦然，故根据最开始设定的特征选取个数，在 S 1 S_1 S1中提取最大最小特征值对应的特征向量。

把λ1中的特征值按照降序排列，则λ2中对应的特征值按升序排列，在λ1中的特征值最大最小各选择m个特征值，将对应的特征向量并整合作为 B ^ \hat B B^

最终，所求的空间滤波器为：
ω = B ^ T P \omega=\hat B^TP ω=B^TP

特征提取

在得到空间滤波器之后，使用该滤波器得到所需要的特征向量，以备后续分类使用。

将原始数据 X N × T X_{N\times T} XN×T使用 ω \omega ω空间滤波
Z 2 m × T = ω 2 m × N X N × T Z_{2m\times T}=\omega_{2m\times N}X_{N\times T} Z2m×T=ω2m×NXN×T

选取特征向量 f p i f_p^i fpi
f p i = log ⁡ ( v a r p i ∑ p = 1 2 m v a r p i ) f_p^i=\log (\frac{var_p^i}{\sum_{p=1}^{2m}var_p^i}) fpi=log(∑p=12mvarpivarpi)
其中 v a r p i var_p^i varpi是 Z 2 m × T Z_{2m\times T} Z2m×T中第p行的方差

将所有 f p i f_p^i fpi整合为 f i f^i fi即为每个trial EEG数据提取的特征

参考文献

[1]Zheng, X., Li, J., Ji, H., Duan, L., Li, M., Pang, Z., … & Tianhao, G. (2020). Task transfer learning for EEG classification in motor imagery-based BCI system. Computational and mathematical methods in medicine, 2020.

[2]Yilu Xu, Qingguo Wei, Hua Zhang, Ronghua Hu, Jizhong Liu, Jing Hua and Fumin Guo.Transfer Learning Based on Regularized Common Spatial Patterns Using Cosine Similarities of Spatial Filters for Motor-Imagery BCI *[J].Journal of Circuits, Systems, and Computers,:World Scientific Publishing Company,2019.28(7).

[3]特征提取丨共空间模式 Common Spatial Pattern (CSP)

[4]协方差与协方差矩阵

[5][线性代数] 矩阵白化

[6]数字图像处理3之空间滤波（平滑空间滤波器）