数据分析阈值分析 求阈值的数学模型

        本篇主要用来学习使用，阈值函数包括软阈值和硬阈值的介绍及求解，详细内容可以参考文后文章。

1 硬阈值(Hard Thresholding)函数

1.1 硬阈值(Hard Thresholding)函数的符号

        硬阈值(Hard Thresholding)并没有软阈值(Soft Thresholding)那么常见，这可能是因为硬阈值解决的问题是非凸的原因吧。硬阈值与软阈值由同一篇文献提出，硬阈值公式参见文献【1】的式(11)：

        第一次邂逅硬阈值(HardThresholding)是在文献【2】中：

        在查询软阈值(Soft Thresholding)的过程中，搜到了文献【3】，进而看到了提到了文献【4】：

        文献【4】中提到的Fig 1如图所示：

        硬阈值的符号到底表示什么意思呢？以文献【1】符号为例，清晰一点来说就是这样的：

$w$是变量，$\lambda$是阈值。

1.2 硬阈值(HardThresholding)函数的作用

        弄清楚了硬阈值(HardThresholding)的符号表示以后，接下来说一说它的作用。这里主要是参考了软阈值的推导过程，然后作者经过一番琢磨和推导而得。

        硬阈值(HardThresholding)可以求解如下优化问题：

        其中：

$||X||_0$是求向量$X$的零范数，即向量$X$中非零元素的个数。根据范数的定义，可以将上面优化问题的目标函数拆开：

$|x|_0$的意思是

$N$个独立的形如函数

$f(x)$进一步写为：

$x\not=0$部分，我们知道它的最小值在$x=b$处取得，最小值为$λ$。现在的问题是$λ$与$b^2$到底谁更小？最小者将是函数的最小值。求解不等式$b^2>λ$可得

$x=0$处取得；

$b^2<λ$可得

$x=b$处取得；

        因此

        与前面的硬阈值(Hard Thresholding)对比一下，发现了么？若将上式中的b视为变量，sqrt(λ)视为阈值，上式即为硬阈值(Hard Thresholding)的公式。

        至此，我们可以得到优化问题

的解为

        注： 该式为硬阈值(Hard Thresholding)的矩阵形式，这里的$B$是一个向量，应该是逐个元素分别执行硬阈值函数。

1.3 硬阈值(HardThresholding)的变形

        当优化问题变为

        因为对目标函数乘一个常系数不影响极值点的获得，所以可等价为优化问题

$hard(\pmb{B},\sqrt(2*\lambda)$。

1.4 硬阈值(Hard Thresholding)的MATLAB代码

        硬阈值(Hard Thresholding)的函数代码可以写成专门针对优化问题

        MATLAB函数代码如下（参考了文献【5】倒数第2页）：

function [ hard_thresh ] = hardthresholding( b,lambda )
    sel = (abs(b)>sqrt(lambda));
    hard_thresh = b.*sel;
end

一定要注意：这种写法是针对最开始的优化问题：

但我个人感觉更应该写成这种通用形式：

function [ x ] = hard( b,T )
    sel = (abs(b)>T);
    x = b.*sel;
end

        如此之后，若要解决优化问题

        只需调用hard(B, sqrt(λ))即可；若要解决优化问题

        只需调用hard(B, sqrt(2*λ))即可。

1.5 硬阈值(HardThresholding)测试代码

        硬阈值(Hard Thresholding)要解决的优化问题目标函数是非凸的，不太常见，手边目前没有其它函数求解这个问题，因此测试代码只能测一下这个函数编写的正确与否了：

clear all;close all;clc; 
b = [-0.8487   -0.3349    0.5528    1.0391   -1.1176]';
lambda = 0.5;
x1=hardthresholding(b,lambda)
x2=hard(b,sqrt(lambda))
fprintf('\nError between hardthresholding and hard = %f\n',norm(x1-x2))

        这里就不给出输出结果了。可以运行一下，从输出结果来看，函数的功能是正确的。

        另外，可以在matlab里输入以下命令看一个软阈值的图像：

x=-5:0.01:5;T=1;y=hard(x,T);plot(x,y);grid;

1.6 参考文献

【1】Donoho D L, JohnstoneJ M. Ideal spatial adaptation by wavelet shrinkage[J]. Biometrika, 1994, 81(3):425-455.

【2】Wright SJ, Nowak R D, Figueiredo M A T. Sparse reconstruction by separableapproximation[J]. IEEE Transactions on Signal Processing, 2009, 57(7):2479-2493.

【3】http://blog.sina.com.cn/s/blog_6d0e97bb01015vq3.html

【4】Elad M,Figueiredo M A T, Ma Y. On the Role of Sparse and Redundant Representations inImage Processing[J]. Proceedings of the IEEE, 2010, 98(6):972-982.

【5】http://www.docin.com/p-553314466.html

---

2 软阈值函数(Soft Thresholding)

2.1 软阈值(Soft Thresholding)函数的符号

        软阈值(Soft Thresholding)目前非常常见，文献【1】【2】最早提出了这个概念。软阈值公式的表达方式归纳起来常见的有三种，以下是各文献中的软阈值定义符号：

文献【1】式(12)：

文献【2】：

文献【3】：

文献【4】式(8)：

文献【5】式(1.5)：

文献【6】式(12)注释：

文献【7】：

无论是sgn(x)还是sign(x)都是符号函数，即当$x>0$时为1，当$x<0$时为-1。)

$w$是变量，$λ$是阈值(非负值)，符号$(|w|-λ)_+$表示当$(|w|-λ)>0$时则等于$|w|-λ$，当$(|w|-λ)<0$时则等于0。那么分三种情况来讨论：第一种情况是$w>\lambda>0$，则$sgn(w)=1$，$|w|=w$，$(|w|-λ)$一定大于0，$(|w|-λ)_+=|w|-λ$，所以$η_S(w,λ)=w-λ$；第二种情况是$w<-λ<0$，则$sgn(w)=-1$，$|w|=-w$，$(|w|-λ)$也一定大于0，$(|w|-λ)_+=|w|-λ$，所以$η_S(w,λ)=-1*(-w-λ)= w+λ$；第三种情况是$|w|<λ$，此时$(|w|-λ)$一定小于0，则$(|w|-λ)_+=0$，所以$η_S(w,λ)=0$。因此

$max{|u|-a,0}$的作用与第一种表示方式中的符号$(|w|-λ)_+$的作用一样，即当$(|u|-a)>0$时$max{|u|-a,0}=(|u|-a)$，当$(|u|-a)<0$时，$max{|u|-a,0}=0$，知道了这一点剩下的分析与第一种表示方式相同。

           综上，三种表示方式均是一致的。

2.2 软阈值(Soft Thresholding)函数的作用

        弄清楚了软阈值(Soft Thresholding)的符号表示以后，接下来说一说它的作用。以下内容主要参考了文献【7】，这是一个非常棒的PPT！！！

        软阈值(SoftThresholding)可以求解如下优化问题：

        其中：

        根据范数的定义，可以将上面优化问题的目标函数拆开：

$N$个独立的形如函数

$f(x)$导数：

        这里要解释一下变量x绝对值的导数，当x>0时，|x|=x，因此其导数等于1；当x<0时，|x|=-x，因此其导数等于-1；综合起来，x绝对值的导数等于sgn(x)。令函数f(x)导数等于0，得：

        这个结果等号两端都有变量x，需要再化简一下。下面分三种情况讨论：

(1)当b>λ/2时

        假设x<0，则sgn(x)=-1，所以x=b+λ/2>0，与假设x<0矛盾；

        假设x>0，则sgn(x)=1，所以x=b-λ/2>0，成立；

        所以此时在x=b-λ/2>0处取得极小值：

        即此时极小值小于f(0)，而当x<0时

        即当x<0时函数f(x)为单调降函数（对任意△x<0，f(0)<f(△x)）。因此，函数在x=b-λ/2>0处取得最小值。

(2)当b<-λ/2时

        假设x<0，则sgn(x)=-1，所以x=b+λ/2<0，成立；

        假设x>0，则sgn(x)=1，所以x=b-λ/2<0，与假设x<0矛盾；

        所以此时在x=b+λ/2<0处取得极小值：

        即此时极小值小于f(0)，而当x>0时

        即当x>0时函数f(x)为单调升函数（对任意△x>0，f(△x)>f(0)）。因此，函数在x=b+λ/2<0处取得最小值。

(3)当-λ/2<b<λ/2时(即|b|<λ/2时)

        假设x<0，则sgn(x)=-1，所以x=b+λ/2>0，与假设x<0矛盾；

        假设x>0，则sgn(x)=1，所以x=b-λ/2<0，与假设x<0矛盾；

        即无论x为大于0还是小于0均没有极值点，那么x=0是否为函数f(x)的极值点呢？

        对于△x≠0，

        当△x >0时，利用条件b<λ/2可得

        当△x <0时，利用条件b<λ/2可得(注：此时|△x |=-△x)

        因此，函数在x=0处取得极小值，也是最小值。

        综合以上三种情况，f(x)的最小值在以下位置取得：

        与前面的软阈值(Soft Thresholding)对比一下，发现了么？若将上式中的b视为变量，λ/2视为阈值，上式即为软阈值(SoftThresholding)的公式。

        至此，我们可以得到优化问题

的解为

        注：该式为软阈值(Soft Thresholding)的矩阵形式。

2.3 软阈值(Soft Thresholding)的变形

        当优化问题变为

        因为对目标函数乘一个常系数不影响极值点的获得，所以可等价为优化问题

此时的解为soft(B, λ)。

2.4 软阈值(Soft Thresholding)的MATLAB代码

        软阈值(Soft Thresholding)的函数代码可以写成专门针对优化问题

        软阈值(Soft Thresholding)是如此简单以至于可以用一句代码去实现它[8]：

        当然，如果不习惯这种形式，也可以写成常见的函数形式：

function [ soft_thresh ] = softthresholding( b,lambda )
    soft_thresh = sign(b).*max(abs(b) - lambda/2,0);
end

        一定要注意：这种写法是针对最开始的优化问题：

        但我个人感觉更应该写成这种通用形式：

function [ x ] = soft( b,T )
    x = sign(b).*max(abs(b) - T,0);
end

        如此之后，若要解决优化问题

只需调用soft(B, λ/2)即可；若要解决优化问题

只需调用soft(B, λ)即可。

2.5 软阈值(Soft Thresholding)测试代码

        用以下一小段代码测试一下软阈值，用来求解优化问题：

压缩感知重构算法之基追踪降噪(Basis PursuitDe-Noising, BPDN)

clear all;close all;clc; 
b = [-0.8487   -0.3349    0.5528    1.0391   -1.1176]';
lambda = 1;
x1=soft(b,lambda)
x2=BPDN_quadprog(b,eye(length(b)),lambda)
fprintf('\nError between soft and BPDN = %f\n',norm(x1-x2))

        这里就不给出输出结果了。运行后，观察输出结果可知，soft函数与BPDN_quadprog函数的输结果相同。

        另外，可以在matlab里输入以下命令看一个软阈值的图像：

x=-5:0.1:5;T=1;y=soft(x,T);plot(x,y);grid;

2.6 总结

        可以发现，软阈值解决的优化问题和基追踪降噪问题很像，但并不一样，而且需要格外说明的是，软阈值并能不解决基追踪降噪问题，文献【8】在最后明确说明了这一点：

2.7 参考文献

【1】Donoho D L, JohnstoneJ M. Ideal spatial adaptation by wavelet shrinkage[J]. Biometrika, 1994, 81(3):425-455.

【2】Donoho D L.De-noising by soft-thresholding[J]. IEEE transactions on information theory,1995, 41(3): 613-627.

【3】Bredies K, Lorenz D.Iterative soft-thresholding converges linearly[R]. Zentrum fürTechnomathematik, 2007.

【4】Bioucas-Dias J M,Figueiredo M A T. A new TwIST: two-step iterative shrinkage/thresholdingalgorithms for image restoration[J]. IEEE Transactions on Image processing,2007, 16(12): 2992-3004.

【5】Beck A, Teboulle M. Afast iterative shrinkage-thresholding algorithm for linear inverse problems[J].SIAM journal on imaging sciences, 2009, 2(1): 183-202.

【6】Wright S J, Nowak RD, Figueiredo M A T. Sparse reconstruction by separable approximation[J]. IEEETransactions on Signal Processing, 2009, 57(7): 2479-2493.

【7】谷鹄翔.IteratedSoft-Thresholding Algorithm[Report,slides]. http://www.sigvc.org/bbs/thread-41-1-2.html

【8】http://www.simonlucey.com/soft-thresholding/

【9】http://blog.sina.com.cn/s/blog_6d0e97bb01015vq3.html

---