目标跟踪之ADMM求解简介

近两年，基于相关滤波的目标跟踪开始大量利用ADMM进行求解，从BACF、STRCF、ARCF、ASRCF到AutoTrack，基本上是一脉相承。SRDCF虽然采用高斯塞尔德进行迭代求解，但效率低下，但其实SRDCF也是可以利用ADMM求解的，速度还快很多。
通过观察近两年的论文公式推导，比如BACF、STRCF，AutoTrack我发现论文中给的公式与代码的公式略有差异，主要是表现在系数上多了个1/T。本文主要是以BACF、SRDCF求解为例，简要分析一下ADMM的求解思路，有很多不足之处，仅供参考，共同学习，欢迎指正。

一、ADMM基本原理

ADMM的基本形式如下：
argminx,yf(x)+g(y)arg \underset{x,y}{min} f(x)+g(y)argx,yminf(x)+g(y)
s.tAx+By=Cs.t Ax+By=Cs.tAx+By=C
要是换成以前的思路，我肯定是利用x,y的关系，求出x等于多少y，然后消去y变量，对f(x)+g(y)f(x)+g(y)f(x)+g(y)进行求导，转为对变量x的求解，但ADMM的思路主要是将两个变量的约束条件（Ax+By=C）引入问题求解，将大问题的最优求解变为两个子问题，即对x，y分开求解，迭代数次，得到最优解。
首先构造增广拉格朗日公式：
Lρ(x,z,y)=f(x)+g(z)+ζT(Ax+Bz−C)+ρ2∥Ax+By−C∥22L_{\rho }(x,z,y)=f(x)+g(z)+\zeta^{T}(Ax+Bz-C)+\frac{\rho }{2}\left \| Ax+By-C \right \|_{2}^{2}Lρ(x,z,y)=f(x)+g(z)+ζT(Ax+Bz−C)+2ρ∥Ax+By−C∥22
ρ2∥Ax+By−C∥22\frac{\rho }{2}\left \| Ax+By-C \right \|_{2}^{2}2ρ∥Ax+By−C∥22的引入主要是为了使上诉凸函数变得更凸，可以这么理解，越凸越好，也因为该项的存在，所以叫增广拉格朗日。接下来，将问题分为子问题进行迭代求解：
xk+1=argminx,yLρ(xk,yk,ςk)=f(x)+ζT(Ax+Bz−C)+ρ2∥Ax+By−C∥22x^{k+1}=arg \underset{x,y}{min} L_{\rho }(x^{k},y^{k},\varsigma^{k})=f(x)+\zeta^{T}(Ax+Bz-C)+\frac{\rho }{2}\left \| Ax+By-C \right \|_{2}^{2}xk+1=argx,yminLρ(xk,yk,ςk)=f(x)+ζT(Ax+Bz−C)+2ρ∥Ax+By−C∥22
yk+1=argminx,yLρ(xk+1,yk,ςk)=g(y)+ζT(Axk+1+By−C)+ρ2∥Axk+1+By−C∥22y^{k+1}=arg \underset{x,y}{min} L_{\rho }(x^{k+1},y^{k},\varsigma^{k})=g(y)+\zeta^{T}(Ax^{k+1}+By-C)+\frac{\rho }{2}\left \| Ax^{k+1}+By-C \right \|_{2}^{2}yk+1=argx,yminLρ(xk+1,yk,ςk)=g(y)+ζT(Axk+1+By−C)+2ρ∥∥Axk+1+By−C∥∥22
ςk+1=ςk+ρ(Axk+1+Byk+1−C)\varsigma^{k+1}=\varsigma^{k}+\rho (Ax^{k+1}+By^{k+1}-C)ςk+1=ςk+ρ(Axk+1+Byk+1−C)
注：
xk+1x^{k+1}xk+1第k+1次迭代得到的x
xkx^{k}xk第k次迭代得到的x
通过事先设定好迭代次数，既可以求出最优解

二、ADMM与相关滤波的求解

相关滤波是一种利用图像连续两帧相关程度来进行目标跟踪的方法，近年来，相关滤波主要面临的问题是边界效应，解决方法主要有两种主流：一是以SRDCF为代表的正则项约束，另一个是以BACF为代表的区域约束，其他方法大都是以这两个算法为baseline进行修改。
这两种方法都可以使用ADMM进行求解，详情可以参阅下面两篇博客，介绍的很详细。
BACF跟踪算法与代码一致的推导
SRDCF相关滤波跟踪算法中ADMM的使用

经典的相关滤波误差函数如下：
E(f)=12∥y−∑k=1Kfk∗xk∥22+λ2∑k=1K∥fk∥22E(f)=\frac{1}{2}\left \|y- \sum_{k=1}^{K}f_{k}*x_{k}\right \|_{2}^{2}+\frac{\lambda }{2}\sum_{k=1}^{K}\left \| f_{k} \right \|_{2}^{2}E(f)=21∥∥∥y−∑k=1Kfk∗xk∥∥∥22+2λ∑k=1K∥fk∥22
f：滤波器模板，fkf_{k}fk第k维特征对应的滤波器模板
K：总共有K维特征
但是该方法为解决训练样本不足的缺点，引入循环矩阵产生大量样本，但也带来边界效应，主要是因为循环矩阵的拼接引入的，解决方法主要有：

2.1 BACF

BACF缩小训练样本尺寸，保证尽可能多的真实样本
主要体现在剪切矩阵P，使用矩阵P来对有边界效应的大图像进行剪切，产生小样本但是真实的训练图像
误差函数：
E(f)=12∑j=1T∥y(j)−∑k=1KfkT∗Pxk[Δτj]∥22+λ2∑k=1K∥fk∥22E(f)=\frac{1}{2} \sum_{j=1}^{T} \left \|y(j)- \sum_{k=1}^{K}f_{k}^{T}*Px_{k}[\Delta \tau_{j }] \right \|_{2}^{2}+\frac{\lambda }{2}\sum_{k=1}^{K}\left \| f_{k} \right \|_{2}^{2}E(f)=21∑j=1T∥∥∥y(j)−∑k=1KfkT∗Pxk[Δτj]∥∥∥22+2λ∑k=1K∥fk∥22
T：样本x的尺寸大小，可以理解为有多少像素
P：D×T的矩阵，使样本x的尺寸从T变为D（T>>D）
Δτj\Delta \tau_{j }Δτj：循环矩阵的第j个样本
由于误差函数只有一个变量f，实际可以直接求解，在KCF、CSK当中，就是直接对f求导，计算最佳的f使得误差函数最小，但从SRDCF开始，不再这么简单粗暴，而是通过构造辅助变量，引入g，将问题求解变成两个子问题
辅助变量的构造主要有两种：
1、g−f=0g-f=0g−f=0
对应的增广拉格朗日：
E(f,g)=12∑j=1T∥y(j)−∑k=1KgkT∗Pxk[Δτj]∥22+λ2∑k=1K∥fk∥22+ςT∑k=1K(fk−gk)+ρ2∑k=1K∥fk−gk∥22E(f,g)=\frac{1}{2} \sum_{j=1}^{T} \left \|y(j)- \sum_{k=1}^{K}g_{k}^{T}*Px_{k}[\Delta \tau_{j }] \right \|_{2}^{2}+\frac{\lambda }{2}\sum_{k=1}^{K}\left \| f_{k} \right \|_{2}^{2}+\varsigma ^{T}\sum_{k=1}^{K}(f_{k}-g_{k})+\frac{\rho }{2}\sum_{k=1}^{K}\left \| f_{k}-g_{k}\right \|_{2}^{2}E(f,g)=21∑j=1T∥∥∥y(j)−∑k=1KgkT∗Pxk[Δτj]∥∥∥22+2λ∑k=1K∥fk∥22+ςT∑k=1K(fk−gk)+2ρ∑k=1K∥fk−gk∥22
2、g^−TFf\hat{g}-\sqrt{T}Ffg^−TFf
注意F与f不一样，F是正交矩阵，FFT=IFF^{T}=IFFT=I,而f是滤波器该辅助变量我个人认为更加适合相关滤波的求解，直接将涉及到卷积的子问题抛到频域，而没有卷积操作的子问题留在空域求解
对应的增广拉格朗日：
L(f,g^)=12T∑j=1T∥y(j)−∑k=1Kg^kTPxk[Δτj]∥22+λ2∑k=1K∥fk∥22+ς^T∑k=1K(g^−T(FPT⊗IK)f)+ρ2∑k=1K∥g^−T(FPT⊗IK)f∥22L(f,\hat{g})=\frac{1}{2T} \sum_{j=1}^{T} \left \|y(j)- \sum_{k=1}^{K}\hat{g}_{k}^{T}Px_{k}[\Delta \tau_{j }] \right \|_{2}^{2}+\frac{\lambda }{2}\sum_{k=1}^{K}\left \| f_{k} \right \|_{2}^{2}+\hat{\varsigma} ^{T}\sum_{k=1}^{K}(\hat{g}-\sqrt{T}(FP^{T}\otimes I_{K})f)+\frac{\rho }{2}\sum_{k=1}^{K}\left \| \hat{g}-\sqrt{T}(FP^{T}\otimes I_{K})f\right \|_{2}^{2}L(f,g^)=2T1∑j=1T∥∥∥y(j)−∑k=1Kg^kTPxk[Δτj]∥∥∥22+2λ∑k=1K∥fk∥22+ς^T∑k=1K(g^−T(FPT⊗IK)f)+2ρ∑k=1K∥∥∥g^−T(FPT⊗IK)f∥∥∥22
注意，第一项系数为12T\frac{1}{2T}2T1 很多论文推导公式的时候都没有1/T，但在代码中又无缘无故多出了1/T，主要是在这里少了1/T
另外ADMM还有一种重要形式（scale form）（ASRCF、AutoTrack、STRCF均采用这种形式）：
引入变量vT=ςTρv^{T}=\frac{\varsigma ^{T}}{\rho}vT=ρςT
增广拉格朗日变为：
L(f,g^)=12T∑j=1T∥y(j)−∑k=1Kg^kTPxk[Δτj]∥22+λ2∑k=1K∥fk∥22+ρ2∑k=1K∥g^−T(FPT⊗IK)f+v^T∥22−ρ2∥v^T∥22L(f,\hat{g})=\frac{1}{2T} \sum_{j=1}^{T} \left \|y(j)- \sum_{k=1}^{K}\hat{g}_{k}^{T}Px_{k}[\Delta \tau_{j }] \right \|_{2}^{2}+\frac{\lambda }{2}\sum_{k=1}^{K}\left \| f_{k} \right \|_{2}^{2}+\frac{\rho }{2}\sum_{k=1}^{K}\left \| \hat{g}-\sqrt{T}(FP^{T}\otimes I_{K})f +\hat{v}^{T}\right \|_{2}^{2}-\frac{\rho }{2}\left \| \hat{v}^{T} \right \|_{2}^{2}L(f,g^)=2T1∑j=1T∥∥∥y(j)−∑k=1Kg^kTPxk[Δτj]∥∥∥22+2λ∑k=1K∥fk∥22+2ρ∑k=1K∥∥∥g^−T(FPT⊗IK)f+v^T∥∥∥22−2ρ∥∥v^T∥∥22
但因为最后一项ρ2∥v^T∥22\frac{\rho }{2}\left \| \hat{v}^{T} \right \|_{2}^{2}2ρ∥∥v^T∥∥22即不包含g^\hat{g}g^，也不包含hhh，所以论文中通常将其舍去。
这两种形式都可以得到最优解，实际使用没有谁优谁劣之说
接下来采用形式1进行求解，先将求和求掉：
E(f,g^)=12T∥y^−g^X^∥22+λ2∥f∥22+ς^T(g^−T(FPT⊗IK)f)+ρ2∥g^−T(FPT⊗IK)f∥22E(f,\hat{g})=\frac{1}{2T}\left \|\hat{y}- \hat{g}\hat{X} \right \|_{2}^{2}+\frac{\lambda }{2}\left \| f \right \|_{2}^{2}+\hat{\varsigma} ^{T}(\hat{g}-\sqrt{T}(FP^{T}\otimes I_{K})f)+\frac{\rho }{2}\left \| \hat{g}-\sqrt{T}(FP^{T}\otimes I_{K})f\right \|_{2}^{2}E(f,g^)=2T1∥∥∥y^−g^X^∥∥∥22+2λ∥f∥22+ς^T(g^−T(FPT⊗IK)f)+2ρ∥∥∥g^−T(FPT⊗IK)f∥∥∥22
其实和原式子是一样的，只不过这里的g^X^\hat{g}\hat{X}g^X^相当于是K维的集合，换个方式写-----（在matlab代码里面都不需要考虑这么多）
子问题g^\hat{g}g^的求解
g^=argming^12T∥y^−g^X^∥22+ς^T(g^−T(FPT⊗IK)f)+ρ2∥g^−T(FPT⊗IK)f∥22\hat{g}=arg \underset{\hat{g}}{min} \frac{1}{2T}\left \|\hat{y}- \hat{g}\hat{X} \right \|_{2}^{2}+\hat{\varsigma} ^{T}(\hat{g}-\sqrt{T}(FP^{T}\otimes I_{K})f)+\frac{\rho }{2}\left \| \hat{g}-\sqrt{T}(FP^{T}\otimes I_{K})f\right \|_{2}^{2}g^=argg^min2T1∥∥∥y^−g^X^∥∥∥22+ς^T(g^−T(FPT⊗IK)f)+2ρ∥∥∥g^−T(FPT⊗IK)f∥∥∥22
对g^\hat{g}g^求导并使其为0，解得：
g^=y^X^T−Tς^T+Tρf^X^X^T+TρIK\hat{g}=\frac{\hat{y}\hat{X}^{T}-T\hat{\varsigma} ^{T}+T\rho\hat{f}}{\hat{X}\hat{X}^T+T\rho I_{K}}g^=X^X^T+TρIKy^X^T−Tς^T+Tρf^
注：f^=T(FPT⊗IK)f\hat{f}=\sqrt{T}(FP^{T}\otimes I_{K})ff^=T(FPT⊗IK)f
但该公式在matlab中的计算量仍然庞大，可以使用Sherman-Morrison公式进一步化简：
(A+uvT)−1=A−1−A−1uvTA−11+vTA−1u(A+uv^{T})^{-1}=A^{-1}-\frac{A^{-1}uv^{T}A^{-1}}{1+v^{T}A^{-1}u}(A+uvT)−1=A−1−1+vTA−1uA−1uvTA−1
在此处A=Tρ,u=X^,vT=X^TA=T\rho,u=\hat{X},v^{T}=\hat{X}^{T}A=Tρ,u=X^,vT=X^T解得：
g^=1ρT(I−X^X^TρT+X^X^T)(y^X^T−Tς^T+Tρf^)\hat{g}=\frac{1}{\rho T}(I-\frac{\hat{X}\hat{X}^{T}}{\rho T+\hat{X}\hat{X}^{T}})(\hat{y}\hat{X}^{T}-T\hat{\varsigma} ^{T}+T\rho\hat{f})g^=ρT1(I−ρT+X^X^TX^X^T)(y^X^T−Tς^T+Tρf^)
子问题f的求解：
f=argminfλ2∥f∥22+ς^T(g^−T(FPT⊗IK)f)+ρ2∥g^−T(FPT⊗IK)f∥22f=arg \underset{f}{min} \frac{\lambda }{2}\left \| f \right \|_{2}^{2}+\hat{\varsigma} ^{T}(\hat{g}-\sqrt{T}(FP^{T}\otimes I_{K})f)+\frac{\rho }{2}\left \| \hat{g}-\sqrt{T}(FP^{T}\otimes I_{K})f\right \|_{2}^{2}f=argfmin2λ∥f∥22+ς^T(g^−T(FPT⊗IK)f)+2ρ∥∥∥g^−T(FPT⊗IK)f∥∥∥22
同样对f求导并使其为零，解得：
f=ρTg+ς^TTλ+ρTf=\frac{\rho Tg+\hat{\varsigma} ^{T}T}{\lambda+\rho T}f=λ+ρTρTg+ς^TT
之后进行多次迭代完成求解，代码如下：

 %迭代次数为2while (i <= params.admm_iterations)%   solve for G- please refer to the paper for more detailsB = S_xx + (T * mu);S_lx = sum(conj(model_xf) .* l_f, 3);S_hx = sum(conj(model_xf) .* h_f, 3);g_f = (((1/(T*mu)) * bsxfun(@times, yf, model_xf)) - ((1/mu) * l_f)  + h_f) - ...bsxfun(@rdivide,(((1/(T*mu)) * bsxfun(@times, model_xf, (S_xx .* yf))) ...- ((1/mu) * bsxfun(@times, model_xf, S_lx)) + (bsxfun(@times, model_xf, S_hx))), B);%   solve for Hh = (T/((mu*T)+ params.admm_lambda))  *   ifft2((mu*g_f) + l_f);     [sx,sy,h] = get_subwindow_no_window(h, floor(use_sz/2) , small_filter_sz);t = single(zeros(use_sz(1), use_sz(2), size(h,3)));t(sx,sy,:) = h;h_f = fft2(t);%   update Ll_f = l_f + (mu * (g_f - h_f));%   update mu- betha = 10.mu = min(betha * mu, mumax);i = i+1;end

注：
代码中的mu对应公式中的ρ\rhoρ
代码中的T取值为2500（50×50）
l_f对应公式中的ς^T\hat{\varsigma} ^{T}ς^T
h、h_f对应公式中的fff和f^\hat{f}f^

2.2 SRDCF的ADMM求解如下：

SRDCF解决边界效应主要是引入正则项惩罚w，w是一个倒立的高斯函数（近两年都是使用矩形），边界值比较大，惩罚滤波器边界的系数，即使训练出来的滤波器更加关注图像中心的细节，忽视边缘不真实的细节。
误差函数：
E(f)=12∥y−∑k=1Kfk∗xk∥22+12∑k=1K∥wfk∥22E(f)=\frac{1}{2}\left \|y- \sum_{k=1}^{K}f_{k}*x_{k}\right \|_{2}^{2}+\frac{1}{2}\sum_{k=1}^{K}\left \| wf_{k} \right \|_{2}^{2}E(f)=21∥∥∥y−∑k=1Kfk∗xk∥∥∥22+21∑k=1K∥wfk∥22
增广拉格朗日:
L(f,g^)=12T∥y^−g^X^∥22+12∥wf∥22+ς^T(g^−TFf)+ρ2∥g^−TFf∥22L(f,\hat{g})=\frac{1}{2T}\left \|\hat{y}- \hat{g}\hat{X} \right \|_{2}^{2}+\frac{1 }{2}\left \| wf \right \|_{2}^{2}+\hat{\varsigma} ^{T}(\hat{g}-\sqrt{T}Ff)+\frac{\rho }{2}\left \| \hat{g}-\sqrt{T}Ff\right \|_{2}^{2}L(f,g^)=2T1∥∥∥y^−g^X^∥∥∥22+21∥wf∥22+ς^T(g^−TFf)+2ρ∥∥∥g^−TFf∥∥∥22
子问题g^\hat{g}g^：
g^=argming^12T∥y^−g^X^∥22+ς^T(g^−TFf)+ρ2∥g^−TFf∥22\hat{g}=arg \underset{\hat{g}}{min}\frac{1}{2T}\left \|\hat{y}- \hat{g}\hat{X} \right \|_{2}^{2}+\hat{\varsigma} ^{T}(\hat{g}-\sqrt{T}Ff)+\frac{\rho }{2}\left \| \hat{g}-\sqrt{T}Ff\right \|_{2}^{2}g^=argg^min2T1∥∥∥y^−g^X^∥∥∥22+ς^T(g^−TFf)+2ρ∥∥∥g^−TFf∥∥∥22
解得：
g^=1ρT(I−X^X^TρT+X^X^T)(X^Ty^−Tς^T+ρTf^)\hat{g}=\frac{1}{\rho T}(I-\frac{\hat{X}\hat{X}^{T}}{\rho T+\hat{X}\hat{X}^{T}})(\hat{X}^{T}\hat{y}-T\hat{\varsigma} ^{T}+\rho T \hat{f})g^=ρT1(I−ρT+X^X^TX^X^T)(X^Ty^−Tς^T+ρTf^)
子问题h

h=argminh12∥wf∥22+ς^T(g^−TFf)+ρ2∥g^−TFf∥22h=arg \underset{h}{min} \frac{1 }{2}\left \| wf \right \|_{2}^{2} +\hat{\varsigma} ^{T}(\hat{g}-\sqrt{T}Ff)+\frac{\rho }{2}\left \| \hat{g}-\sqrt{T}Ff\right \|_{2}^{2}h=arghmin21∥wf∥22+ς^T(g^−TFf)+2ρ∥∥∥g^−TFf∥∥∥22
解得：
h=TFς^T+ρTFg^wTw+ρT=TςT+ρTgwTw+ρTh=\frac{\sqrt{T}F\hat{\varsigma} ^{T}+\rho \sqrt{T}F\hat{g}}{w^{T}w+\rho T} =\frac{T{\varsigma} ^{T}+\rho T{g}}{w^{T}w+\rho T}h=wTw+ρTTFς^T+ρTFg^=wTw+ρTTςT+ρTg