【目标跟踪】KCF高速跟踪详解

Henriques, João F., et al. “High-speed tracking with kernelized
correlation filters.” Pattern Analysis and Machine Intelligence, IEEE
Transactions on 37.3 (2015): 583-596.

本文的跟踪方法效果甚好，速度奇高，思想和实现均十分简洁。其中利用循环矩阵进行快速计算的方法尤其值得学习。另外，作者在主页上十分慷慨地给出了各种语言的实现代码。
本文详细推导论文中的一系列步骤，包括论文中未能阐明的部分。请务必先参看这篇简介循环矩阵性质的博客。

思想

一般化的跟踪问题可以分解成如下几步：

1. 在$I_t$帧中，在当前位置$p_t$附近采样，训练一个回归器。这个回归器能计算一个小窗口采样的响应。
2. 在$I_{t+1}$帧中，在前一帧位置$p_t$附近采样，用前述回归器判断每个采样的响应。
3. 响应最强的采样作为本帧位置$p_{t+1}$。

循环矩阵表示图像块

在图像中，循环位移操作可以用来近似采样窗口的位移。

训练时，围绕着当前位置进行的一系列位移采样可以用二维分块循环矩阵$X$表示，第ij块表示原始图像下移i行右移j列的结果。类似地，测试时，前一帧结果附近的一系列位移采样也可以用$X$表示。

这样的$X$可以利用傅里叶变换快速完成许多线性运算。

线性回归训练提速

此部分频繁用到了循环矩阵的各类性质，请参看这篇博客。
线性回归的最小二乘方法解为：
$$w={\left(X^{H}X+\lambda I\right)}^{-1}{X}^{H}y$$

根据循环矩阵乘法性质，$X^{H}X$的特征值为$\hat{x}\odot {\hat{x}}^{\ast }$。$I$本身就是一个循环矩阵，其生成向量为$[1,0,0...0]$，这个生成向量的傅里叶变换为全1向量，记为$\delta$。
$$w={\left(Fdiag(\hat{x}\odot {\hat{x}}^{\ast })F^H+\lambda Fdiag(\delta )F^H\right)}^{-1}{X}^{H}y$$

$$={\left(Fdiag(\hat{x}\odot {\hat{x}}^{\ast }+\lambda \delta )F^H\right)}^{-1}{X}^{H}y$$

根据循环矩阵求逆性质，可以把矩阵求逆转换为特征值求逆。
$$w=F\cdot diag\left(\frac{1}{\hat{x}\odot {\hat{x}}^{\ast }+\lambda \delta }\right)\cdot F^H{X}^{H}y$$

$$w=F\cdot diag\left(\frac{1}{\hat{x}\odot {\hat{x}}^{\ast }+\lambda \delta }\right)\cdot F^H\cdot Fdiag({\hat{x}}^{\ast })F^H\cdot y$$

利用$F$的酉矩阵性质消元：
$$w=F\cdot diag\left(\frac{{\hat{x}}^{\ast }}{\hat{x}\odot {\hat{x}}^{\ast }+\lambda \delta }\right)\cdot F^H\cdot y$$

分号表示用1进行对位相除。
反用对角化性质：$Fdiag(y)F^H=C({\mathcal{F}}^{-1}(y))$，上式的前三项还是一个循环矩阵。
$$w=C\left({\mathcal{F}}^{-1}\left(\frac{{\hat{x}}^{\ast }}{\hat{x}\odot {\hat{x}}^{\ast }+\lambda \delta }\right)\right)\cdot y$$

利用循环矩阵卷积性质$\mathcal{F}(C(x)\cdot y)={\hat{x}}^{\ast }\odot \hat{y}$：
$$\mathcal{F}(w)={\left(\frac{{\hat{x}}^{\ast }}{\hat{x}\odot {\hat{x}}^{\ast }+\lambda \delta }\right)}^{\ast }\odot \mathcal{F}(y)$$
由于$\hat{x}\odot {\hat{x}}^{\ast }$的每个元素都是实数，所以共轭不变：

$$\mathcal{F}(w)=\frac{\hat{x}}{\hat{x}\odot {\hat{x}}^{\ast }+\lambda \delta }\odot \mathcal{F}(y)=\frac{\hat{x}\odot \hat{y}}{\hat{x}\odot {\hat{x}}^{\ast }+\lambda \delta }$$

论文中，最后这一步推导的分子部分写成${\hat{x}}^{\ast }\odot \hat{y}$，是错误的。但代码中没有涉及。

线性回归系数$\omega$可以通过向量的傅里叶变换和对位乘法计算得到。

核回归训练提速

不熟悉核方法的同学可以参看这篇博客的简单说明。核回归方法的回归式为：
$$f(z)={\alpha }^T\kappa (z)$$
其中$\kappa (z)$表示测试样本$z$和所有训练样本的核函数。参数有闭式解：
$$\alpha ={\left(K+\lambda I\right)}^{-1}y$$

$K$为所有训练样本的核相关矩阵：$K_{ij}=\kappa (x_i,x_j)$。如果核函数选择得当，使得$x$内部元素顺序更换不影响核函数取值，则可以保证$K$也是循环矩阵。以下核都满足这样的条件：

设核相关矩阵的生成向量是$k$。推导和之前线性回归的套路非常类似：
$$\alpha ={\left(Fdiag(\hat{k})F^H+Fdiag(\lambda \delta )F^H\right)}^{-1}y={\left(Fdiag(\hat{k}+\lambda \delta )F^H\right)}^{-1}y$$

$$=Fdiag\left(\frac{1}{\hat{k}+\lambda \delta }\right)F^{H}y=C\left({\mathcal{F}}^{-1}\left(\frac{1}{\hat{k}+\lambda \delta }\right)\right)y$$

利用循环矩阵卷积性质$\mathcal{F}(C(x)\cdot y)={\hat{x}}^{\ast }\odot \hat{y}$：
$$\hat{\alpha }={\left(\frac{1}{\hat{k}+\lambda \delta }\right)}^{\ast }\odot \hat{y}$$
这里$k$是核相关矩阵的第一行，表示原始生成向量$x^0$和移位了$i$的向量$x^i$的核函数。考察其处于对称位置上的两个元素：
$$k_i=\kappa (x^0,x^i),k_{N-i}=\kappa (x^0,x^{N-i})$$

两者都是同一个向量和自身位移结果进行运算。因为所有涉及到的核函数都只和位移的绝对值有关，所以$k_i=k_{N-i}$，即$k$是对称向量。

举例：$x^0=[1,2,3,4]$，$x^1=[4,1,2,3]$，$x^3=[2,3,4,1]$。使用多项式核$\kappa (x,y)=x^{T}y$，容易验证$\kappa (x^0,x^1)=\kappa (x^0,x^3)$。

对称向量的傅里叶变换为实数，有：

$$\hat{\alpha }=\left(\frac{1}{\hat{k}+\lambda \delta }\right)\odot \hat{y}=\frac{\hat{y}}{\hat{k}+\lambda \delta }$$

论文中，利用$k$的对称性消除共轭的步骤没有提及。

线性回归系数$\alpha$可以通过向量的傅里叶变换和对位乘法计算得到。

核回归检测提速

所有待检测样本和所有训练样本的核相关矩阵为$K$，每一列对应一个待测样本。可以一次计算所有样本的响应（$N\times 1$向量）：

$$y^{\prime }=K^T\alpha$$

利用循环矩阵的转置性质性质，$C^{(}k)$的特征值为${\hat{k}}^{\ast }$：
$$y^{\prime }=C(k{)}^T\cdot \alpha =C({\hat{k}}^{\ast })\cdot \alpha =k^{\ast }\ast \alpha$$

利用循环矩阵的卷积性质：
$$y^{\prime }=(k^{\ast })\ast \alpha =k\ast \alpha$$

两边傅里叶变换：
$$\widehat{y^{\prime }}=\hat{k}\odot \hat{\alpha }$$

论文中，利用转置消除共轭的步骤没有提及。

所有侯选块的检测响应可以通过向量的傅里叶变换和对位乘法计算得到。

核相关矩阵计算提速

无论训练还是检测，都需要计算核相关矩阵$K$的生成向量$k$。除了直接计算每一个核函数，在某些特定的核函数下可以进一步加速。

多项式核

$$\kappa (x,y)=f\left(x^{T}y\right)$$

其中$f$为多项式函数。写成矩阵形式：
$$K=f\left(X^{T}Y\right)$$

$f$在矩阵的每个元素上单独进行。根据循环矩阵性质，$X^{T}Y$也是一个循环矩阵，其生成向量为${\mathcal{F}}^{-1}(\hat{y}\odot {\hat{x}}^{\ast })$。所以核相关矩阵的生成向量为：
$$k=f\left({\mathcal{F}}^{-1}(\hat{y}\odot {\hat{x}}^{\ast })\right)$$

RBF核

$$\kappa (x,y)=f(||x-y|{|}^2)$$
其中$f$是线性函数。简单展开：
$$\kappa (x,y)=f(||x-y|{|}^2)=f(||x|{|}^2+||y|{|}^2+2x^{T}y)$$
由于$X$中的所有$x$都通过循环移位获得，故$||x|{|}^2$对于所有$x$是常数，同理$||y|{|}^2$也是。所以核相关矩阵的生成向量为：
$$k=f\left(||x|{|}^2+||y|{|}^2+{\mathcal{F}}^{-1}(\hat{y}\odot {\hat{x}}^{\ast })\right)$$

其他核

有一些核函数，虽然能保证$K$是循环矩阵，但无法直接拆解出其特征值，快速得到生成向量。比如Hellinger核：${\sum }_i\sqrt{x_i{y}_i}$，Intersection核：${\sum }_i\min (x_i,y_i)$。

多通道

在多通道情况下（例如使用了HOG特征），生成向量$x$变成$M\times L$，其中$M$是样本像素数，$L$是特征维度。在上述所有计算中，需要更改的只有向量的内积：
$$x^{T}y=\sum _l(x^l{)}^T{y}^l$$

注：非常感谢GX1415926535和大家的帮助，发现原文一处错误。(21)式中不应有转置，应为：
$$f(z)=K^z\alpha$$