共空间模式CSP

共空间模式（common spatial pattern，CSP）是脑-机接口领域常用的一类空间滤波算法，尤其在运动想象范式分类上具有较好的效果，是运动想象范式的基准算法之一。

目前，CSP及其改进算法的发展速度放缓，看似到达了算法的瓶颈期，近几年鲜少有较大的突破。尽管如此，CSP中的一些思想对脑-机接口算法设计仍然具有一定的启示作用。本文从CSP原始算法出发，讨论其变形和一系列改进算法，试图阐明其中的数学思想与神经科学的联系。

CSP数学原理

原始形式

2000年Graz的论文中提出的CSP是为2分类问题设计的，形式较为简单，然而如果你读CSP相关论文，就会发现CSP存在至少三种表述形式。这三种方式相互联系，又有所区分，很容易让初学者陷入混乱，不知道哪一种是正确形式。我接下来从2000年Graz的论文中的算法出发，讨论三种形式间的联系和不同。

假设我们采集脑电数据为$\{(\mathbf{X}^{(i)},y^{(i)})\}_{i=1}^{N_t}$，其中$\mathbf{X}^{(i)} \in \mathbb{R}^{N_c \times N_s}$是第$i$个样本，$N_c$是EEG导联个数，$N_s$是采样时间点的个数，$y^{(i)}$是$i$个样本的标签，$N_t$为总样本个数。第$i$个样本的协方差矩阵$\mathbf{C}^{(i)} \in \mathbb{R}^{N_c \times N_c}$为（所有样本均经过零均值处理）:

$$\begin{equation} \mathbf{C}^{(i)} = \frac{1}{N_t-1}\mathbf{X}^{(i)}\left(\mathbf{X}^{(i)}\right)^T \end{equation}$$

第$l$类的平均协方差矩阵$\overline{\mathbf{C}}\vphantom{C}^l$为：

$$\begin{equation} \overline{\mathbf{C}}\vphantom{C}^l = \frac{1}{|\mathcal{I}_l|} \sum_{i \in \mathcal{I}_l} \frac{\mathbf{C}^{(i)}}{\mathrm{tr}\left(\mathbf{C}^{(i)}\right)} \end{equation}$$

其中$\mathcal{I}_l$是标签为$l$的样本索引集合，$|\mathcal{I}_l|$则是集合中样本的个数，$\mathrm{tr}\left(\cdot\right)$求矩阵的迹。

为什么要使用$\mathrm{tr}\left(\cdot\right)$来对协方差矩阵实现迹归一化？

1990年Koles等人的文章中指出，迹归一化的目的是为了消除"被试间脑电信号幅值的变化"，注意到Koles等人的主要目的是区分健康人群和精神疾病人群，而个体的脑电幅值是有差异的。方差可以表征信号在时域上的能量高低，不同人群的协方差矩阵的绝对值不同。为了消除这种差异带来的影响，利用$\mathrm{tr}()$函数求得所有导联的总体能量，并对协方差矩阵迹归一化，从而安排除不同个体带来的干扰。Graz小组对同一个体不同试次的数据沿用了这种归一化方式，试图消除试次间的差异，发现也有一定的作用，这种迹归一化方式就一直流传下来。

然而，有些分析显示这种归一化方式会不利于最终的空间滤波器排序，建议不要使用迹归一化。实践中使不使用迹归一化还是要具体问题具体分析。我的感觉是没有必要在这里加入迹归一化，因为很多时候EEG预处理阶段已经使用了各种归一化手段来减弱噪声的影响。

接下来构建复合协方差矩阵$\overline{\mathbf{C}}\vphantom{C}_c$，并特征值分解，构建白化（whitening）矩阵$\mathbf{P}$：

$$\begin{equation} \begin{split} \overline{\mathbf{C}}\vphantom{C}^c &= \overline{\mathbf{C}}\vphantom{C}^1 + \overline{\mathbf{C}}\vphantom{C}^2\\ &= \mathbf{U}^c \mathbf{\Lambda}^c \left(\mathbf{U}^c\right)^T\\ \mathbf{P} &= \left(\mathbf{\Lambda}^c\right)^{-1/2}\left(\mathbf{U}^c\right)^T \end{split} \end{equation}$$

其中$\mathbf{U}^c$是特征向量矩阵（每一列是特征向量），$\mathbf{\Lambda}^c$是由特征值组成的对角矩阵。$\mathbf{P}$是白化矩阵，使得$\mathbf{P}\overline{\mathbf{C}}\vphantom{C}^c\left(\mathbf{P}\right)^T= \mathbf{I}$成立，注意到：

$$\begin{equation} \begin{split} \mathbf{I} &= \mathbf{P}\overline{\mathbf{C}}\vphantom{C}^c\left(\mathbf{P}\right)^T\\ &= \mathbf{P}\overline{\mathbf{C}}\vphantom{C}^1\left(\mathbf{P}\right)^T + \mathbf{P}\overline{\mathbf{C}}\vphantom{C}^2\left(\mathbf{P}\right)^T\\ &= \mathbf{S}^1 + \mathbf{S}^2\\ \mathbf{S}^1 &= \mathbf{P}\overline{\mathbf{C}}\vphantom{C}^1\left(\mathbf{P}\right)^T\\ \mathbf{S}^2 &= \mathbf{P}\overline{\mathbf{C}}\vphantom{C}^2\left(\mathbf{P}\right)^T \end{split} \end{equation}$$

对$\mathbf{S}^1$或$\mathbf{S}^2$做特征值分解，得到最终的空间滤波器$\mathbf{W}$：

$$\begin{equation} \begin{split} \mathbf{S}^1 &= \mathbf{U} \mathbf{\Lambda}^1 \left(\mathbf{U}\right)^T\\ \mathbf{S}^2 &= \mathbf{U} \mathbf{\Lambda}^2 \left(\mathbf{U}\right)^T\\ \mathbf{W} &= \mathbf{P}^T\mathbf{U} \end{split} \end{equation}$$

其中$\mathbf{S}^1$和$\mathbf{S}^2$具有相同的特征向量$\mathbf{U}$（这也是共空间模式名称的由来），这里假设$\mathbf{U}$的每一列是按照$\mathbf{\Lambda}^1$中的特征值从大到小排列的，可以看出$\mathbf{\Lambda}^2$中的特征值是从小到大排列的，满足$\mathbf{\Lambda}^1+\mathbf{\Lambda}^2=\mathbf{I}$的关系。

为什么$\mathbf{S}^1$和$\mathbf{S}^2$具有同样的特征向量和此消彼长的特征值关系？
这一点可以简单的证明如下：
假设$\mathbf{u}_j$和$\lambda_j^{1}$分别是$\mathbf{S}^1$的特征向量和特征值，即：

$$\begin{equation} \mathbf{S}^1\mathbf{u}_j=\lambda_j^{1}\mathbf{u}_j \end{equation}$$

注意到$\mathbf{S}^1+\mathbf{S}^2=\mathbf{I}$，把上式中的$\mathbf{S}^1$置换掉可得：

$$\begin{equation} \left(\mathbf{I}-\mathbf{S}^2\right)\mathbf{u}_j=\lambda_j^{1}\mathbf{u}_j \end{equation}$$

把上式变形一下可得：

$$\begin{equation} \mathbf{S}^2\mathbf{u}_j=(1-\lambda_j^{1})\mathbf{u}_j \end{equation}$$

显然$\mathbf{u}_j$也是$\mathbf{S}^2$的特征向量，只不过其特征值为$1-\lambda_j^{1}$。

脑-机接口中的空间滤波器是一组作用于EEG导联信号的向量，目的是为了加强空间分辨率或信噪比，可以简单理解为对导联信号的线性组合。事实上，不少空间滤波器本质上就是某些特征值分解问题的特征向量。

以上就是原始CSP算法的基本内容，在得到空间滤波器矩阵$\mathbf{W}$后（$\mathbf{W}$的每一列都是一个空间滤波器），选择前后各$m$个空间滤波器构建特征向量$\tilde{\mathbf{x}}$如下：

$$\begin{equation} \begin{split} \tilde{\mathbf{W}} &= \begin{bmatrix} \mathbf{w}_1, \cdots, \mathbf{w}_m, \mathbf{w}_{N_c-m+1}, \cdots, \mathbf{w}_{N_c} \end{bmatrix}\\ \tilde{\mathbf{x}} &= \mathrm{log}\left(\frac{\mathrm{diag}\left(\tilde{\mathbf{W}}^T\mathbf{X}\mathbf{X}^T\tilde{\mathbf{W}}\right)}{\mathrm{tr}\left(\tilde{\mathbf{W}}^T\mathbf{X}\mathbf{X}^T\tilde{\mathbf{W}}\right)}\right) \end{split} \end{equation}$$

其中$\mathbf{w}_m$表示$\mathbf{W}$的第$m$列，$\tilde{\mathbf{W}}$是最终选定的空间滤波器组，$\mathrm{diag}\left(\cdot\right)$是矩阵主对角线上的元素，$\mathrm{log}\left(\cdot\right)$对每个元素做对数变换，其主要目的是使数据近似正太分布。获得特征向量$\tilde{\mathbf{x}}$后，则可以使用线性判别分析（Linear Discriminant Analysis，LDA）、支持向量机（Support Vector Machine，SVM）等常见的机器学习模型构建分类器。

以上就是原始CSP算法的基本内容，简单回顾一下CSP算法，不难发现CSP实质求解的是这样一个问题，寻找正交矩阵$\mathbf{W}$对角化$\overline{\mathbf{C}}\vphantom{C}^1$和$\overline{\mathbf{C}}\vphantom{C}^2$，使得以下条件成立：

$$\begin{equation} \begin{split} \mathbf{W}^T\overline{\mathbf{C}}\vphantom{C}^1\mathbf{W} &= \mathbf{\Lambda}^1\\ \mathbf{W}^T\overline{\mathbf{C}}\vphantom{C}^2\mathbf{W} &= \mathbf{\Lambda}^2\\ \mathbf{\Lambda}^1 + \mathbf{\Lambda}^2 &= \mathbf{I}\\ \end{split} \end{equation}$$

让我们对以上的公式做一些变换，把第一个和第二个公式相加：

$$\begin{equation} \begin{split} \mathbf{W}^T\left(\overline{\mathbf{C}}\vphantom{C}^1+\overline{\mathbf{C}}\vphantom{C}^2\right)\mathbf{W} &= \mathbf{\Lambda}^1 + \mathbf{\Lambda}^2\\ &= \mathbf{I}\\ \end{split} \end{equation}$$

又因为$\mathbf{W}$是正交矩阵，故$\mathbf{W}^{-1}=\mathbf{W}^T$，从而：

$$\begin{equation} \left(\overline{\mathbf{C}}\vphantom{C}^1+\overline{\mathbf{C}}\vphantom{C}^2\right)\mathbf{W} = \mathbf{W} \end{equation}$$

把上式代入$\overline{\mathbf{C}}\vphantom{C}^1\mathbf{W}=\mathbf{W}\mathbf{\Lambda}^1$右边的$\mathbf{W}$，可得：

$$\begin{equation} \overline{\mathbf{C}}\vphantom{C}^1\mathbf{W} = \left(\overline{\mathbf{C}}\vphantom{C}^1+\overline{\mathbf{C}}\vphantom{C}^2\right)\mathbf{W}\mathbf{\Lambda}^1 \end{equation}$$

这个式子看起来很像特征向量定义的公式$\overline{\mathbf{C}}\vphantom{C}^1\mathbf{W}=\mathbf{W}\mathbf{\Lambda}^1$，只不过等式右边多了一个矩阵$\overline{\mathbf{C}}\vphantom{C}^1+\overline{\mathbf{C}}\vphantom{C}^2$。这类形式的特征值求解问题叫广义特征值问题，求解广义特征值问题是脑-机接口领域传统空间滤波方法的基础，大量的算法都可以转化为这一形式。

第二种形式

CSP的第二种形式与$\overline{\mathbf{C}}\vphantom{C}^1\mathbf{W} = \left(\overline{\mathbf{C}}\vphantom{C}^1+\overline{\mathbf{C}}\vphantom{C}^2\right)\mathbf{W}\mathbf{\Lambda}^1$密切相关，首先我们需要了解一个数学概念广义雷利商（generalized Rayleigh quotient）。

广义雷利商$\lambda$长这样：

$$\begin{equation} \begin{split} \lambda =\frac{\mathbf{w}^T\mathbf{A}\mathbf{w}}{\mathbf{w}^T\mathbf{B}\mathbf{w}}\\ \mathbf{A}, \mathbf{B} \succeq 0\\ \end{split} \end{equation}$$

其中$\mathbf{A}$和$\mathbf{B}$为半正定矩阵，$\mathbf{w}$是列向量。

如果我们求如下广义雷利商的优化问题，就会有一些有趣的结果：

$$\begin{equation} \begin{split} \max_{\mathbf{w}} \frac{\mathbf{w}^T\mathbf{A}\mathbf{w}}{\mathbf{w}^T\mathbf{B}\mathbf{w}} \end{split} \end{equation}$$

寻找$\mathbf{w}$使得$\lambda$最大，通常令$\mathbf{w}^T\mathbf{B}\mathbf{w}=1$，在数学上可以等价为求解下式：

$$\begin{equation} \mathbf{A}\mathbf{w} = \lambda\mathbf{B}\mathbf{w} \end{equation}$$

这个公式就是上一节提到的广义特征值问题，也就是说，寻找$\mathbf{w}$使广义雷利商最大的优化问题可以等价为求解$\mathbf{A}$和$\mathbf{B}$的广义特征值问题。如果我们继续寻找能够使$\lambda$第二大、第三大的$\mathbf{w}$，就会发现只要解出广义特征值问题的矩阵形式即可：

$$\begin{equation} \mathbf{A}\mathbf{W} = \mathbf{B}\mathbf{W}\mathbf{\Lambda} \end{equation}$$

不难发现，上一节中推导的CSP求解问题可以变形为求解广义雷利商问题：

$$\begin{equation} \begin{split} \overline{\mathbf{C}}\vphantom{C}^1\mathbf{W} = \left(\overline{\mathbf{C}}\vphantom{C}^1+\overline{\mathbf{C}}\vphantom{C}^2\right)\mathbf{W}\mathbf{\Lambda}^1 \ \ \Longleftrightarrow \ \ \argmax_{\mathbf{W}} \frac{\mathrm{tr}\left(\mathbf{W}^T\overline{\mathbf{C}}\vphantom{C}^1\mathbf{W}\right)}{\mathrm{tr}\left(\mathbf{W}^T\left(\overline{\mathbf{C}}\vphantom{C}^1+\overline{\mathbf{C}}\vphantom{C}^2\right)\mathbf{W}\right)} \end{split} \end{equation}$$

其中应满足$\mathbf{W}^T\left(\overline{\mathbf{C}}\vphantom{C}^1+\overline{\mathbf{C}}\vphantom{C}^2\right)\mathbf{W}=\mathbf{I}$的约束条件。

这就是CSP常见的第二种形式，它跟原始形式在数学上相互等价，由于分母在约束下是单位矩阵，也常写作如下优化问题：

$$\begin{equation} \begin{split} \argmax_{\mathbf{W}}\quad &\mathrm{tr}\left(\mathbf{W}^T\overline{\mathbf{C}}\vphantom{C}^1\mathbf{W}\right)\\ \textrm{s.t.}\quad &\mathbf{W}^T\left(\overline{\mathbf{C}}\vphantom{C}^1+\overline{\mathbf{C}}\vphantom{C}^2\right)\mathbf{W} = \mathbf{I}\\ \end{split} \end{equation}$$

第三种形式

CSP的第三种表述形式需要绕点弯路。首先还是从CSP的原始形式出发，即寻找正交矩阵$\mathbf{W}$使得以下条件成立：

$$\begin{equation} \begin{split} \mathbf{W}^T\overline{\mathbf{C}}\vphantom{C}^1\mathbf{W} &= \mathbf{\Lambda}^1\\ \mathbf{W}^T\overline{\mathbf{C}}\vphantom{C}^2\mathbf{W} &= \mathbf{\Lambda}^2\\ \mathbf{\Lambda}^1 + \mathbf{\Lambda}^2 &= \mathbf{I}\\ \end{split} \end{equation}$$

在第二个公式的左右两边同时右乘矩阵$\mathbf{W}^{-1}\left(\overline{\mathbf{C}}\vphantom{C}^2\right)^{-1}$，可以得到：

$$\begin{equation} \mathbf{W}^T=\mathbf{\Lambda}^2\mathbf{W}^{-1}\left(\overline{\mathbf{C}}\vphantom{C}^2\right)^{-1} \end{equation}$$

将该式带入$\mathbf{W}^T\overline{\mathbf{C}}\vphantom{C}^1\mathbf{W} = \mathbf{\Lambda}^1$，替换掉$\mathbf{W}^T$，可得：

$$\begin{equation} \mathbf{\Lambda}^2\mathbf{W}^{-1}\left(\overline{\mathbf{C}}\vphantom{C}^2\right)^{-1}\overline{\mathbf{C}}\vphantom{C}^1\mathbf{W} = \mathbf{\Lambda}^1 \end{equation}$$

左右两边左乘$\overline{\mathbf{C}}\vphantom{C}^2 \mathbf{W} \left(\mathbf{\Lambda}^2\right)^{-1}$，有：

$$\begin{equation} \begin{split} \overline{\mathbf{C}}\vphantom{C}^1\mathbf{W} &= \overline{\mathbf{C}}\vphantom{C}^2 \mathbf{W} \left(\mathbf{\Lambda}^2\right)^{-1} \mathbf{\Lambda}^1\\ &= \overline{\mathbf{C}}\vphantom{C}^2 \mathbf{W} \mathbf{\Lambda}\\ \mathbf{\Lambda} &= \left(\mathbf{\Lambda}^2\right)^{-1} \mathbf{\Lambda}^1\\ \end{split} \end{equation}$$

没错，我们又推出了熟悉的广义特征值问题，再考虑广义雷利商与之的联系，可以得到CSP的第三种形式：

$$\begin{equation} \begin{split} \argmax_{\mathbf{W}}\quad &\mathrm{tr}\left(\mathbf{W}^T\overline{\mathbf{C}}\vphantom{C}^1\mathbf{W}\right)\\ \textrm{s.t.}\quad &\mathbf{W}^T\overline{\mathbf{C}}\vphantom{C}^2\mathbf{W} = \mathbf{I}\\ \end{split} \end{equation}$$

相比CSP的原始形式和第二种形式，第三种形式更适合从直观上解释CSP在运动想象上有效的原因。运动想象会产生事件相关同步（ERS）和事件相关去同步（ERD）的现象，简单来说就是从电信号上看，某些脑区能量升高，某些脑区能量降低，故能量变化才是运动想象分类的关键特征。而方差可以看作一导信号能量的高低（协方差矩阵则是多导信号的综合反应），因此CSP的第三种形式实质体现的是这样一个问题：

寻找一种变换方式$\mathbf{w}$，使得变换后任务1的能量（$\mathbf{w}^T\overline{\mathbf{C}}\vphantom{C}^1\mathbf{w}$）和任务2的能量（$\mathbf{w}^T\overline{\mathbf{C}}\vphantom{C}^2\mathbf{w}$）差异最大化（其比值最大）。

CSP的这种特性恰好和运动想象产生的神经机制变化现象一致，CSP对能量特征做转换，从而强化了不同任务间能量的差异。

关于CSP的第三种形式，最后还需要注意的一点是其同CSP的第二种形式（或原始形式）并不完全等价，我们在推导第三种形式过程种始终没有用到这样一个约束条件$\mathbf{\Lambda}^1 + \mathbf{\Lambda}^2 = \mathbf{I}$。

这表明，第三种形式是CSP的一种泛化形式，其和CSP原始形式和第二种表述的差异仅在于特征值$\Lambda$不要求在0~1的范围内，具体来说，它们的特征值间存在这样一种关系：

$$\begin{equation} \begin{split} \mathbf{\Lambda} &= \left(\mathbf{\Lambda}^2\right)^{-1} \mathbf{\Lambda}^1\\ \mathbf{\Lambda}^1 &= (\mathbf{\Lambda} + I)^{-1}\mathbf{\Lambda}\\ \mathbf{\Lambda}^2 &= (\mathbf{\Lambda} + I)^{-1}\\ \end{split} \end{equation}$$

实现分析

CSP作为经典算法有各种实现，这里主要分析MNE的CSP源码，看看有啥可以学习的地方。

空间滤波器的选择

基本上，目前CSP算法中m的选择方案大多是根据经验选择（通常选择2～4）个。MNE的CSP选择了第二种形式的CSP算法，最后求解的特征值范围在0～1之间，因此可以对$\left| \lambda_i^1 - 0.5 \right|$先排序再取前M个成分的特征向量组成空间滤波器组$\tilde{\mathbf{W}}$（与前后各m个的做法有些许差别，但实践中很难有显著性上的差异，这种做法相对方便一些）。

实际上，针对2分类CSP算法，特征值与类平均协方差矩阵间黎曼距离在各个特征向量分量上的投影长度密切相关，具体的证明细节可以看Alexandre Barachant的这篇会议，对于$\overline{\mathbf{C}}\vphantom{C}^1$和$\overline{\mathbf{C}}\vphantom{C}^2$，有如下关系：

$$\begin{equation} \mathrm{\delta}_R\left(\overline{\mathbf{C}}\vphantom{C}^1,\overline{\mathbf{C}}\vphantom{C}^2\right) = \sqrt{\sum_{i=1}^{N_c} \mathrm{log}^2\left(\frac{\lambda_i^1}{1-\lambda_i^1}\right)} \end{equation}$$

也就是说，我们可以选定一个阈值$\epsilon$，将特征值及特征向量按$\mathrm{log}^2\left(\frac{\lambda}{1-\lambda}\right)$从大到小排序，选取最小的M使下式成立：

$$\begin{equation} \frac{\sqrt{\sum_{i=1}^{M} \mathrm{log}^2\left(\frac{\lambda_i^1}{1-\lambda_i^1}\right)}}{\mathrm{\delta}_R\left(\overline{\mathbf{C}}\vphantom{C}^1,\overline{\mathbf{C}}\vphantom{C}^2\right)} \ge \epsilon \end{equation}$$

再将前M个特征向量组成空间滤波器组$\tilde{\mathbf{W}}$。当$\epsilon=0.9$时，我们可以说选定的空间滤波器组可以贡献大约90%左右的2类之间的黎曼距离。

最后一种常用的空间滤波器选择方法是计算空间滤波后的特征同标签之间的互信息，再按互信息从大到小排列选择前M个空间滤波器。互信息常用于CSP的衍生算法FBCSP的特征筛选过程，也可以用于CSP（MNE中的CSP也提供了互信息的排序选项）。至于哪一种选择方法是最优的，目前似乎还没有定论（我感觉）。

协方差矩阵正则化

CSP中的正则化方法主要是对协方差矩阵做正则化处理，从十几年前开始，BCI研究者就在协方差矩阵的正则化处理上做了大量的工作，有些正则化方法与BCI的变异性问题也有着密切的联系，因此这一方面的正则化方法展开来讲就收不住啦。我们这里介绍的正则化方法的目的非常单纯，就是为了解决EEG中可能存在的协方差矩阵非正定的问题。

一般而言，本文的第一个公式$\mathbf{C}^{(i)}$在大多数情况下都是正定的。所谓矩阵$\mathbf{M}$是正定（definite-positive）的，是指对任意非0实向量$\mathbf{z}$，$\mathbf{z}^T\mathbf{M}\mathbf{z}$都是一个正数。不过在实践中，经常会遇到程序报类似b matrix is not definite positive这种的错误，这种情况来源于底层的特征值分解或广义特征值分解函数对于矩阵的正定性有较为严格的要求，但输入的协方差矩阵却不是正定的。

那么为啥协方差矩阵不是正定的呢？大概率可按照以下三种情况逐步排查：

• EEG采样信号$\mathbf{X}^{(i)}$中的$N_c \gt N_s$，即导联多于采样点
• 虽然导联多于采样点，但$\mathbf{X}^{(i)}$的秩小于$N_c$，即可能存在导联打串或做过共平均参考变换等导致矩阵不满秩的预处理操作
• 前两条都不满足，则可能是由于数值计算精度上出了问题，比如单个样本的协方差矩阵满足正定性，但平均协方差矩阵却不是正定的（我也不太懂数值计算方面的内容，这一条有待考证，但确实碰到过这样的现象）

总之，为了让计算进行下去，对协方差矩阵做正则化处理是很有必要的，协方差矩阵的正则化就是对协方差矩阵做以下变换：

$$\begin{equation} (1-\lambda) * \mathbf{C} + \lambda * \mu * \mathbf{I} \end{equation}$$

其中$\lambda$是待估计的正则化系数，$\mu=\frac{\mathrm{tr}(\mathbf{C})}{N_c}$是为了对单位矩阵的数值范围做限定。sklearn的covariance模块中列出了多种正则化处理方法，比如著名的ledoit-wolf正则化、oas正则化等方法，选个顺眼的用就行。