Kernelized Correlation Filters KCF算法原理详解（阅读笔记）

参考博文：

• 【KCF算法解析】High-Speed Tracking with Kernelized Correlation Filters笔记
• KCF论文理解与源码解析
• KCF算法公式推导
• KCF(High-Speed Tracking with Kernelized Correlation Filters)论文详解
• KCF目标跟踪方法分析与总结

预备知识

1. 岭回归

首先温习一下最小二乘法.

在矩阵形式下，最小二乘法的解$w=(X^{T}X{)}^{-1}{X}^{T}Y$，$X$是一个$m\times n$的矩阵，其中每一行代表样本，列代表样本的某一个属性。为了求得$w$，则必须要求矩阵$(X^{T}X)$满秩（满秩是可逆的充要条件）。然而在解决现实问题时，当$X$中的某些属性相关性较大时，$X$不一定满秩（线性相关$\iff$行列式为0$\iff$矩阵不满秩）。除此之外，当$n$的数量大于$m$即样本矩阵的列数大于行数时，$X$也不满秩，因为：

此时，由于$X^{T}X$不满秩，则不可求其逆，因此无法算出最优解。那是否可以求一个近似最优解出来呢？

一个比较直观的想法是，在$X^{T}X$中加上一个单位阵（即只有对角线有元素，其它处全为0的矩阵），该对角阵值的大小由一个参数$\lambda$控制。

这样做的好处是，使得$X^{T}X$变成满秩矩阵，这样是不是就可以求解了？而由于这个满秩矩阵是在原矩阵的基础上添加了扰动项（对角阵）得来的，所以最后求得的$w$一般不是最优解，但是也是比较接近最优解了。（当最优解很难很难求解出时，可以较为简单地求出接近最优解的解一般是一种更好的选择）

最优解$w=(X^{T}X{)}^{-1}{X}^{T}Y$，在加上由参数$\lambda$控制的对角阵后，求得的次优解为$w=(X^{T}X+\lambda I_n{)}^{-1}{X}^{T}Y$ ，其中：

• $I_n$是单位矩阵，对角线全是1，类似于“山岭”。

• $\lambda$是岭系数，改变其数值可以改变单位矩阵对角线的值。

再反推目标函数，可以得到目标函数为：

$\underset{w}{\mathrm{argmin}}||X\cdot w-Y|{|}^2+\lambda ||w|{|}^2$

也就是说，岭回归其实就是在最小二乘法的基础上，加入了对$w$的惩罚项，当$\lambda$值越大时，惩罚强度就越大，生成的模型就越简单，$||w|{|}^2$也被称为L2范数。从整体上来看，岭回归是对最小二乘回归的一种补充，其损失了无偏性，但是获得了较高计算精度和数值稳定性。

这篇岭回归详解 从零开始 从理论到实践讲得可以，还有这篇岭回归原理简单分析与理解。

2. 循环移位和循环矩阵

参考【KCF算法解析】High-Speed Tracking with Kernelized Correlation Filters笔记 理解即可，注意生成循环矩阵的过程用$X=C(x)$表示，其中$x$是原向量，$X$是生成的循环矩阵。

3. 傅里叶对角化

All circulant matrices are made diagonal by the Discrete Fourier Transform (DFT), regardless of the generating vector x.
任意循环矩阵可以被傅里叶变换矩阵对角化。

参照这篇博客循环矩阵傅里叶对角化理解即可，我也没太懂。

4. 方向梯度直方图（HOG）

图像特征：方向梯度直方图 HOG

正文

全文内容概览参考这篇目标跟踪算法——KCF入门详解即可。

公式推导参考这篇KCF算法公式推导即可。

1. 线性回归

1.1. 岭回归

KCF的基本模型是岭回归器。

线性分类器：
$$f(x)=w^{T}x\text{\hspace{1em}(1)}$$
岭回归的目标函数：
$$\underset{w}{\mathrm{argmin}}||X\cdot w-Y|{|}^2+\lambda ||w|{|}^2\text{\hspace{1em}(2)}$$
基于岭回归求得的闭式解（到底什么是“闭合解（Closed-form Solution）”）：
$$w=(X^{T}X+\lambda I_n{)}^{-1}{X}^{T}Y\text{\hspace{1em}(3)}$$
变换到复数域的闭式解：
$$w=(X^{H}X+\lambda I_n{)}^{-1}{X}^{H}Y\text{\hspace{1em}(4)}$$
其中，$X$表示样本矩阵，其中每一行代表每一个样本，每一列代表所有样本的某个特征值，$X^H$表示其复数的共轭转置。

1.2. 基于循环矩阵获取正负样本

假设感兴趣区域的图像patch由一个长度为$n$的向量$x$表示，$x$也被称作基样本。这个很好理解吧，目标跟踪一般是给定一组连续帧，以及需要跟踪的目标在起始帧的位置信息，然后在接下来的连续帧中执行跟踪。那么我们最开始能够获得的目标特征，就用$x$表示。

之后，原文中的一句话提到：

Our goal is to train a classifier with both the base sample (a positive example) and several virtual samples obtained by translating it (which serve as negative examples).

即，作者计划通过基样本来生成正负样本，然后送给分类器来训练，以获得区分目标区域和其它区域的能力。这句话也很好理解吧，就是基于初始帧以及给定的目标信息，获得目标的特征，然后在此基础上生成正负样本，正样本代表需要跟踪的目标，负样本代表背景或者其它干扰信息。这是目标跟踪领域的一个常见思路，典型的如TLD、MOOSE都采用了这种方法。

如何基于基样本生成正负样本呢？这里得用前面提到的循环矩阵。

对于基样本$x$，生成的循环矩阵$X=C(x)$，文中给了个示意图：

OK，此时，基于循环矩阵得到的正负样本，采用岭回归训练线性模型，得到在复数域下的解析解和前面的公式相同，但是公式中字母的意义不同：
$$w=(X^{H}X+\lambda I_n{)}^{-1}{X}^{H}Y\text{\hspace{1em}(5)}$$
此时，$X$表示由初始向量$x$（基样本）生成的循环矩阵，$X^H$表示其复数的共轭转置，即$X^H=(X^{\ast }{)}^T$，$X^{\ast }$是$X$的复数共轭形式。

1.3. 基于傅里叶对角化的求解

还记得前面提到的任意循环矩阵可以被傅里叶变换矩阵对角化吧，那么，用傅里叶对角化，我们就可以更快速地求解。

对于循环矩阵$X$，其可以被傅里叶对角化：
$$X=C(x)=F\mathrm{diag}(\hat{x})F^H\text{\hspace{1em}(6)}$$
其中，$\hat{x}$是$x$的离散傅里叶变换，即$\hat{x}=\mathcal{F}(\vec{x})$。$F$是离散傅里叶变换矩阵，是一个不依赖于$x$的常数矩阵。注意，离散傅里叶变换矩阵是一个酉矩阵，即其满足$F^{H}F=F{F}^H=I$。

之后，经过公式推导，参考KCF算法公式推导，可以得到线性模型参数$w$的求解结果：
$$\hat{w}=\frac{{\hat{x}}^{\ast }\odot \hat{y}}{{\hat{x}}^{\ast }\odot \hat{x}+\lambda }\text{\hspace{1em}(7)}$$
注意，这里求出的$\hat{w}$是$w$的傅里叶变换形式，使用逆向傅里叶变换即可以轻松地在空间域中恢复$w$，这样我们就能在样本线性可分的前提下，基于训练样本，快速地求解出较优解$w$（也就是能够快速地训练好分类模型）。

结果对比：

• 原闭式解$w=(X^{H}X+\lambda I_n{)}^{-1}{X}^{H}Y$的计算复杂度为$O(n^3)$

• 简化后的解$\hat{w}=\frac{{\hat{x}}^{\ast }\odot \hat{y}}{{\hat{x}}^{\ast }\odot \hat{x}+\lambda }$的复杂度为$O(nlogn)$

所以，引入循环矩阵，不仅增强了样本，提高了准确率，而且简化了运算，提高了运行速度。

2. 使用非线性回归对模型进行训练

2.1. 应用kernel-trick的非线性模型的求解结果

上面只考虑了样本矩阵是线性可分的情况，而当数据在当前的特征空间中不能线性可分时，我们可以考虑把数据特征从当前特征空间映射到更高维的空间，从而使得在低维空间不可分的情况到高维空间之后就能变的线性可分了。核函数主要的目的就是把一个线性问题映射到一个非线性核空间，这样把在低维空间不可分的到核空间之后就能够可分了。

• 因此，需要找到一个非线性映射函数$\phi (x)$，使得映射后的样本在高维空间中线性可分。那么在这个高维空间中，就可以使用岭回归来寻找一个分类器$f(x_i)=w^T\phi (x_i)$，其中，$\phi (x_i)$表示对样本$x_i$通过非线性映射函数$\phi$变换到高维空间后的样本向量。

• 对于分类模型$f(x_i)=w^T\phi (x_i)$，$f(x_i)$和$x_i$都已知，$\phi (\cdot )$取决于我们选取的映射函数，因此主要是想求得对$w$的解。在使用Kernel-trick之后，$w$可以用训练样本的线性组合来表示：
  $$w=\sum _i{\alpha }_i\phi (x_i)\text{\hspace{1em}(8)}$$
  则最优化问题不再是求解$w$，而是$\alpha$。至于为什么可以将$w$表示为样本的线性组合，可以参考这篇高维映射 与 核方法（Kernel Methods），里面有详细的解释，建议精读。

• 线性条件下的回归问题，经过非线性变换后，对于新样本$z$的预测值为：
  $$\begin{array}{rl}f(z)& =w^{T}z\\ & =(\sum _i^n{\alpha }_i\phi (x_i){)}^T\cdot \phi (z)\\ & =\sum _i^n{{\alpha }_i\phi (x_i)}^T\phi (z)\\ & =\sum _i^n{\alpha }_i\kappa (x_i,z)\end{array}\text{\hspace{1em}(9)}$$
  其中，$x_i$代表训练分类器时的训练样本，$z$代表新出现的测试样本。即，预测其实是新样本与所有训练样本在高维空间中内积的加权平均。

对于$\alpha$的求解，原文写得比较省略，原文中非线性情况下的最终求解结果（公式16）为
$$\alpha =(K+\lambda I{)}^{-1}y\text{\hspace{1em}(10)}$$
注意，这里的优化目标不是用$w$而是用字母$\alpha$表示，$K$代表核矩阵。

推导过程还是推荐仔细阅读这篇高维映射 与 核方法（Kernel Methods），阅读完之后也能了解什么是核技巧，以及为什么优化目标从线性可分情况下的$w$换成了$\alpha$。

2.2. 快速核方法

当采用基于Kernel trick的非线性模型训练时，最终求解目标$\alpha =(K+\lambda I{)}^{-1}y$包含了核矩阵$K$。

而从核矩阵$K$的公式中可以显而易见地发现，当数据样本量过大时，会导致核矩阵$K$过大，从而使得计算的时间复杂度和空间复杂度都很高。

核方法并非没有缺点。原来的线性模型显然是一种参数化方法，好处是训练完成的模型只需要存储权重向量 。然而在使用了核方法以后，由于我们用训练集替换掉了权重项，因此相当于转化成了非参数化的方法，显然增加了需要存储的数据量，同时每一次做预测都需要用到全部训练数据集。

因此，本节中提到了一种快速的核回归方法。即怎么去计算核函数，怎么快速得到最优解。

原文的思路是：如果能证明核矩阵$K$对于数据集来说是循环矩阵，就可以对其进行对角化，以加速求解过程。

通过施加一个条件来允许K循环。

定理一：给定循环矩阵$C(x)$，如果核函数$\kappa$对于任何置换矩阵$M$都满足$\kappa (x,x^{\prime })=\kappa (Mx,M{x}^{\prime })$，则$C(x)$在核函数$\kappa$下对应的核矩阵$K$是循环矩阵。

证明过程：假定有一个置换矩阵$P$，如：
$$P=\left[\begin{array}{ccccc}0& 0& 0& \cdots & 1\\ 1& 0& 0& \cdots & 0\\ 0& 1& 0& \cdots & 0\\ 0& 0& 1& \cdots & 0\\ \cdots & \cdots & \cdots & \cdots & \cdots \end{array}\right]\text{\hspace{1em}(11)}$$
即对于基样本$x$，有：
$$\begin{gathered} x_1=P^{0}x=x=[x_1,x_2,\cdots ,x_n{]}^T \\ {x}_2=P^{1}x=[x_n,x_1,\cdots ,x_{n-1}{]}^T \\ \cdots \\ {x}_n=P^{n-1}x=[x_2,x_3,\cdots ,x_1{]}^T\text{\hspace{1em}(12)} \end{gathered}$$
则:
$$\begin{array}{rl}{K}_{ij}& =\kappa (x_i,x_j)\\ & =\kappa (P^{i}x,P^{j}x)\end{array}\text{\hspace{1em}(13)}$$
从核矩阵$K$的公式中，可以发现：

核矩阵中第$i$行第$j$列的值就是第$i$和第$j$个样本经过核函数之后得到的值再逐元素相乘计算出的内积。

而我们的样本矩阵是在基向量$x$的基础上经过循环变换得到的循环矩阵，也就是说，这个式子，刚好满足我们之前的所有假设和前提。这里可能有点绕，得细细思考一会。

现在，由于定理中提到的是任意变换$M$，那我们就假设$M$这个置换矩阵为$P^{-i}$。

由于$\kappa (x,x^{\prime })=\kappa (Mx,M{x}^{\prime })$，所以有：
$$\begin{array}{rl}\kappa (M{x}_i,M{x}_j)& =\kappa (x_i,x_j)\\ & =\kappa (P^{-i}{P}^{i}x,P^{-i}{P}^{j}x)\\ & =\kappa (x,P^{j-i}x)\\ & =K_{ij}\end{array}\text{\hspace{1em}(14)}$$
由于$P^n=P^0$，也就是当$j-i$的值每到$n$时为一个循环，所以上式可以改写为$K_{ij}=\kappa (x,P^{(j-i)modn}x)$，即：
$$K=\left[\begin{array}{cccc}\kappa (x,P^{0}x)& \kappa (x,P^{1}x)& \cdots & \kappa (x,P^{n-1}x)\\ \kappa (x,P^{1}x)& \kappa (x,P^{2}x)& \cdots & \kappa (x,P^{0}x)\\ ⋮& ⋮& \ddots & ⋮\\ \kappa (x,P^{n-1}x)& \kappa (x,P^{0}x)& \cdots & \kappa (x,P^{n-2}x)\end{array}\right]\text{\hspace{1em}(15)}$$
可以看到，此时核矩阵$K$是一个循环矩阵。

原文中提到，下面的核函数都满足定理一：

1. 径向基函数核，如：高斯核；
2. 点积核，如：线性核、多项式核；
3. 加权核，如：交集、${\chi }^2$；
4. 指数加权核

当$K$为循环矩阵时，由前面关于循环矩阵的知识，可以得知核矩阵$K$中的所有元素都是基于其第一行这个向量循环移位得到的。令$k^{xx}$为核矩阵$K$中的第一行向量，则$K=C(k^{xx})$。

原文中非线性情况下的最终求解结果为
$$\alpha =(K+\lambda I{)}^{-1}y\text{\hspace{1em}(10)}$$
由于$K$为循环矩阵，则可以应用傅里叶对角化求解：
$$\begin{array}{rl}\alpha & =(K+\lambda I{)}^{-1}y\\ & =(C(k^{xx})+\lambda I{)}^{-1}y\\ & =(F\mathrm{diag}({\hat{k}}^{xx})F^H+\lambda I{)}^{-1}y\\ & =(F\mathrm{diag}({\hat{k}}^{xx}+\lambda )F^H{)}^{-1}y\\ & =(F^H{)}^{-1}(\mathrm{diag}({\hat{k}}^{xx}+\lambda {)}^{-1})F^{-1}y\\ & =F\mathrm{diag}({\hat{k}}^{xx}+\lambda {)}^{-1}{F}^{H}y\\ & =F\mathrm{diag}(\frac{1}{{\hat{k}}^{xx}+\lambda })F^{H}y\end{array}\text{\hspace{1em}(16)}$$
这里用到了离散傅里叶变换矩阵是酉矩阵的性质，对公式两边进行离散傅里叶变换得到：
$$\hat{\alpha }=\frac{\hat{y}}{{\hat{k}}^{xx}+\lambda }\text{\hspace{1em}(17)}$$
从上面这个公式，也就是原文中的公式(17)可以发现，当保证核矩阵$K$为循环矩阵时，通过傅里叶对角化，可以使得对$\alpha$的求解过程只需要涉及到核矩阵$K$中的第一行向量${\hat{k}}^{xx}$（$1\times n$大小），而${\hat{k}}^{xx}$是随着样本数量线性变化的。这意味着什么？传统的核方法需要计算一个$n\times n$的核矩阵，与样本数量呈平方增长。

这意味着快速核方法大大节省了空间资源和时间资源。

2.3. 快速检测

在2.2节中，得到了系数$\alpha$的简化求解形式（训练过程），而为了进一步提高算法的运行速度，文章对检测的过程也进行了简化加速。

首先由训练样本和其对应的标签训练检测器，其中训练集是由目标区域和由其移位得到的若干样本组成，对应的标签是根据距离越近正样本可能性越大的准则赋值的，然后可以得到$\alpha$。

由式(9)可得，当对新样本$z$（待检测的图像patch）进行检测时，有：
$$f(z)=K\alpha \text{\hspace{1em}(18)}$$
注意，岭回归中包含的正则化只在训练时使用，对于预测，直接用求解好的$\alpha$和核矩阵$K$即可。

现在，用$K^z$替换$K$，考虑到实际情况，式(18)变为：
$$f(z)=(K^z{)}^T\alpha \text{\hspace{1em}(19)}$$
其中，$K^z=C(k^{xz})$，表示训练样本和待检测样本之间的核矩阵，是一个非对称矩阵。$k^{xz}$表示由基样本$x$和基础待检测图像patch$z$计算得到的核相关向量，$C$仍代表循环位移操作。

由前面的内容易知，$K^z$中的元素可以被表示为$\kappa (P^{i-1}z,P^{j-1}x)$，$P$的定义如式(11)所示。

基于定理一，同样可以证明训练样本和检测样本之间的核矩阵$K^z$为循环矩阵，因此可以对其进行傅里叶对角化，有：
$$\begin{array}{rl}f(z)& =(K^z{)}^T\alpha \\ & =(C(k^{xz}){)}^T\alpha \\ & =(F\mathrm{diag}({\hat{k}}^{xz})F^H{)}^T\alpha \\ & =F^H\mathrm{diag}({\hat{k}}^{xz})F\alpha \end{array}\text{\hspace{1em}(20)}$$
考虑到离散傅里叶变换矩阵$F$的性质，则有：
$$Ff(z)=\mathrm{diag}({\hat{k}}^{xz})F\alpha \text{\hspace{1em}(21)}$$
舍掉$F$，再进行傅里叶变换，则可以得到最终的快速检测公式：
$$\hat{f}(z)={\hat{k}}^{xz}\odot \hat{\alpha }\text{\hspace{1em}(22)}$$
该公式可以理解为用$\alpha$对核矩阵$K^z$进行空间滤波操作（或者是线性加权），生成的$f(z)$是一个向量，代表所有可能位置的置信度值，最大的置信度值对应的patch即作为检测的输出。

3. 快速核相关

尽管前面已经提出了快速训练和快速检测的方法，但它们仍依赖于计算两个图像块的核相关性这一过程（如$k^{xx}$和$k^{xz}$），根据前面的内容，核相关操作包含了计算两个输入（及所有相对位移）的高维特征点乘值，这是最后一个需要解决的计算瓶颈。因此本节对核相关过程进行了加速。

这一节主要就是介绍几种常见的核函数的形式，以及每种核函数所对应的快速计算方式。

• 点积核：
  $$k^{x{x}^{\prime }}=g({\mathcal{F}}^{-1}({\hat{x}}^{\ast }\odot {\hat{x}}^{\prime }))\text{\hspace{1em}(23)}$$

• 多项式核：
  $$k^{x{x}^{\prime }}=({\mathcal{F}}^{-1}({\hat{x}}^{\ast }\odot {\hat{x}}^{\prime })+a{)}^b\text{\hspace{1em}(24)}$$

• 对于某种径向基函数$h$，其核函数可表示为：
  $$k^{x{x}^{\prime }}=h(||x^{\prime }|{|}^2+||x|{|}^2-2{\mathcal{F}}^{-1}({\hat{x}}^{\ast }\odot {\hat{x}}^{\prime }))\text{\hspace{1em}(25)}$$

• 高斯核：
  $$k^{x{x}^{\prime }}=exp(-\frac{1}{{\sigma }^2}(||x^{\prime }|{|}^2+||x|{|}^2-2{\mathcal{F}}^{-1}({\hat{x}}^{\ast }\odot {\hat{x}}^{\prime })))\text{\hspace{1em}(26)}$$

4. KCF算法流程

1. 在初始帧中基于初始的目标bounding box，扩大2.5倍（跟踪区域是目标大小的2.5倍可以提供上下文信息），得到大小为$M\times N$的ROI。
2. 计算初始图像patch的HOG特征，得到31维的HOG特征图，即HOG特征图的大小为$m\times n\times 31$（$m$和$n$的大小取决于$M$、$N$以及HOG的cell size参数）。
3. 对HOG特征图进行余弦窗口加权，从而减轻由于边界移位导致的图像不光滑。
4. 利用二维高斯函数得到长宽和HOG特征图一样的标签矩阵$y$（即大小也为$m\times n$）。越靠近中心的像素点具有更高的权重，使得训练出的模型更关心目标的中心区域。
5. 基于快速核相关公式得到$k^{xx}$。
6. 针对输入的HOG特征图，训练出非线性回归器：
   $$\hat{\alpha }=\frac{\hat{y}}{{\hat{k}}^{xx}+\lambda }\text{\hspace{1em}(17)}$$
7. 在接下来的一帧中，基于上一帧的目标中心，框选出大小为$M\times N$的ROI，求得HOG特征图，进行余弦窗加权，之后，同样基于快速核相关求得$k^{xz}$。
8. 基于下述公式计算傅域下的响应矩阵，然后计算逆傅里叶变换得到$f(z)$，即，最终的响应矩阵：
   $$\hat{f}(z)={\hat{k}}^{xz}\odot \hat{\alpha }\text{\hspace{1em}(22)}$$
9. 找到响应矩阵$f(z)$中响应值最大的位置，若该位置的响应值大于给定阈值，则输出该位置作为当前全局图像中的目标位置。若最大响应值仍小于给定阈值，则采取补救方法（这个还不太清楚是什么，好像是全局图像搜索，然后重新初始化）。
10. 模型更新， 基于当前输出的目标位置，选取样本，类似于步骤1-6计算得到新的${\alpha }^{\prime }$，然后采用线性插值来更新模型：
    $${\alpha }_{new}=m{\alpha }_{old}+(1-m){\alpha }^{\prime }\text{\hspace{1em}(27)}$$
    $m$为学习率，是一个0-1之间的值。$m$越大说明以前的模型影响越大（太大容易跟丢），越小说明以前的模型影响越小（太小鲁棒性低）。
11. 重复7-10的检测过程。
    

5. 总结

相关滤波算法中最有代表性的算法之一，原理复杂，但是其实现比较简单。
神中神，伟大，无需多言。