同伦法
題目：同倫法（Homotopy Method）
? ? ? ? 學(xué)習(xí)壓縮感知重構(gòu)算法，經(jīng)常能見到同倫法，但這里首先要特別說明的是，今天這里討論的同倫法僅僅是一種思想，而不是一個(gè)具體的算法，類似于Majorization-Minimization優(yōu)化框架。
1、同倫的基本原理
? ? ? ??首先看一下同倫的基本原理【1】：
? ? ? ??這里說的比較專業(yè)，直接簡單一點(diǎn)，例如，令
注意，此時(shí)H(x,0)=f(x)，而H(x,1)=g(x)，滿足上面的式(3-31)，因此H(x,s)就是一條連接f和g的路徑。
2、牛頓迭代與同倫法
? ? ? ??單純知道了同倫法的思想是不夠的，下面參考【2】舉一個(gè)同倫法與牛頓迭代法[3]相結(jié)合的應(yīng)用實(shí)例，以更好的消化同倫法的思想。
2.1牛頓迭代法
? ? ? ??在我的印象中，首次見到牛頓迭代法是在大一時(shí)學(xué)C語言，譚浩強(qiáng)的教材里有一道題是讓編寫一個(gè)牛頓迭代的C語言程序。以下參考【3】簡單介紹一下牛頓迭代法。
? ? ? ??牛頓迭代法（Newton's method）又稱為牛頓-拉夫遜（拉弗森）方法（Newton-Raphsonmethod），它是牛頓在17世紀(jì)提出的一種在實(shí)數(shù)域和復(fù)數(shù)域上近似求解方程的方法。多數(shù)方程不存在求根公式，因此求精確根非常困難，甚至不可能，從而尋找方程的近似根就顯得特別重要。牛頓迭代法使用函數(shù)f(x)的泰勒級數(shù)的前面幾項(xiàng)來尋找方程f(x) = 0的根。如下圖所示，若要求解方程f(x) = 0的根，即求解y=f(x)與x軸的交點(diǎn)x*：
? ? ? ??設(shè)x*是f(x) = 0的根，選取x0作為x*的初始近似值，過點(diǎn)(x0,f(x0))做曲線y=f(x)的切線L，根據(jù)初中學(xué)過的知識很容易得出L的方程為（已知直線斜率和線上的一個(gè)點(diǎn)）
令y=0，求出L與x軸交點(diǎn)的橫坐標(biāo)
稱x1為x*的一次近似值。過點(diǎn)(x1,f(x1))做曲線y=f(x)的切線，并求該切線與x軸交點(diǎn)的橫坐標(biāo)
稱x2為x*的二次近似值。重復(fù)以上過程，得x*的近似值序列，其中
稱為x*的n+1次近似值，上式稱為牛頓迭代公式。
? ? ? ??下面通過一個(gè)例子來說明牛頓迭代法的應(yīng)用。
? ? ? ??已知f(x)=ex-1，求方程f(x) = 0的根。（這個(gè)例子很簡單，不用求也知道方程的根為x=0。）
? ? ? ??牛頓迭代法需要函數(shù)的導(dǎo)數(shù)，因此首先編寫兩個(gè)函數(shù)：
? ? ? ??1)func.m（即f(x)）
?
function y = func(x)
y = exp(x)-1;
end
? ? ? ??2)func_g.m（即f?(x)的導(dǎo)數(shù)）
?
function y = func_g(x)
y = exp(x);
end
? ? ? ??下面給出牛頓迭代的函數(shù)代碼：
? ? ? ??3)NewtonIter.m
?
function [x] = NewtonIter(x0,f,g,tol,verbose,EnPlot,neg_in,pos_in,step_in)
% Version: 1.0 written by jbb0523 @2016-08-23
% 牛頓迭代法函數(shù)
% x0為初值, f為目標(biāo)求解函數(shù), g為f的導(dǎo)數(shù)
% tol為精度要求,verbose和EnPlot控制是否輸出迭代過程
% x為輸出最終解
% 參考文獻(xiàn)：1)百度百科：牛頓迭代法
% 2)http://blog.sina.com.cn/s/blog_86186c970102vll8.html
if nargin < 9
step_in = 0.05;%默認(rèn)繪圖步長
end
if nargin < 8
pos_in = 3;%默認(rèn)繪圖上限
end
if nargin < 7
neg_in = -3;%默認(rèn)繪圖下限
end
if nargin < 6
EnPlot = 0;%默認(rèn)不輸出迭代過程(figure)
end
if nargin < 5
verbose = 0;%默認(rèn)不輸出迭代過程(Command Window)
end
if nargin < 4
tol = 1e-3;
end
if EnPlot
neginf = neg_in;posinf = pos_in;step = step_in;
X_Plot = neginf:step:posinf;
figure;plot(X_Plot,f(X_Plot),'r');
grid on;hold on;pause;
end
iter = 0;
x = x0;
if verbose
fprintf('Iter %d: x=%f,f(x)=%f\n',iter,x,f(x));
end
if EnPlot
Y_Plot = f(x) + g(x)*(X_Plot-x);%切線方程
plot(X_Plot,Y_Plot);
hold on;pause;
end
while abs(f(x)) > tol
x = x - f(x)/g(x);%牛頓迭代公式
iter = iter + 1;
if verbose
fprintf('Iter %d: x=%f,f(x)=%f\n',iter,x,f(x));
end
if EnPlot
Y_Plot = f(x) + g(x)*(X_Plot-x);%切線方程
plot(X_Plot,Y_Plot);
hold on;pause;
end
end
if verbose
hold off;
end
end
?
? ? ? ??4)Test.m
?
close all;clear all;clc;
x0 = 5;
x =NewtonIter(x0,@func,@func_g);
fprintf('x=%f\n',x)
?
? ? ? ??運(yùn)行Test.m后，最終輸出結(jié)果為x=0.00015，這是由于NewtonIter.m中的精度參數(shù)tol默認(rèn)為10-3，若想提到x的精度，請修改tol的值。
? ? ? ??特別注意，這里說的牛頓迭代法并不是最優(yōu)化方法中的牛頓法，若要學(xué)習(xí)最優(yōu)化方法中的牛頓法，參見博客http://blog.csdn.net/itplus/，該博客中作者用了五篇詳細(xì)說明了從牛頓法到擬牛頓法的一系列問題[4]。
2.2牛頓迭代法的不足
? ? ? ??牛頓迭代法需要輸入一個(gè)初值x0，然而這個(gè)初值有時(shí)并不好選，選的不合適還會(huì)造成迭代不收斂，以一個(gè)例子來說明。
? ? ? ??已知f(x)=arctan(x)，求方程f(x) = 0的根。（這個(gè)例子也很簡單，不用求也知道方程的根為x=0。）
? ? ? ??跟前面一樣，首先編寫兩個(gè)函數(shù)：
? ? ? ??1)fx.m（即f(x)）
?
function y = fx(x)
%f(x)=arctan(x)
y = atan(x);
end
?
? ? ? ??2)gx.m（即f?(x)的導(dǎo)數(shù)）
?
function y = gx(x)
%g(x)=1/(1+x^2)
y = 1./(1+x.^2);
end
?
? ? ? ??3)Test.m
?
close all;clear all;clc;
x0 = 5;
x =NewtonIter(x0,@fx,@gx,1e-3,1,1);
fprintf('x=%f\n',x)
?
? ? ? ??我們發(fā)現(xiàn)輸出結(jié)果為x=NaN，并沒有得到方程的根，這是因?yàn)橛捎诔踔颠x的不合適，而這個(gè)函數(shù)又有其特殊之處，因此最終迭代未收斂。文獻(xiàn)【2】指出，若初始值
? ? ? ??則迭代就會(huì)發(fā)散。可以在測試程序中將x0=1.3，則可以得到輸出結(jié)果為x=-0.000026。
2.3 將同倫法應(yīng)用于牛頓迭代
? ? ? ??文獻(xiàn)【2】指出，同倫法可以用來幫助牛頓迭代產(chǎn)生一個(gè)好的初始值：
? ? ? ??構(gòu)建函數(shù)F0(x)和H(x,s)：
則函數(shù)F0(x)與F(x)同倫，同倫路徑為H(x,s)。尋找一個(gè)滿足條件函數(shù)F0(x)很簡單，例如：
其中x*是一個(gè)給定的起始值，把x=x*代入F0(x)，則有F0(x*)=0 。將F0(x)代入H(x,s)得：
? ? ? ??將同倫法應(yīng)用于牛頓迭代工作流程如下：
? ? ? ??因此，為了“已知f(x)=arctan(x)，求方程f(x) = 0的根”，我們編程如下：
? ? ? ??1)Hxs.m（即H(x,s)）
?
function y = Hxs(x,s,x0)
%H(x,s)=arctan(x) +(s-1)*arctan(x0)
y = atan(x) + (s-1)*atan(x0);
end
?
? ? ? ??2)gx.m（與前面一樣，因?yàn)镠(x,s)相比于f(x)只加了一個(gè)常數(shù)，其導(dǎo)數(shù)不變）
?
function y = gx(x)
%g(x)=1/(1+x^2)
y = 1./(1+x.^2);
end
?
? ? ? ??由于應(yīng)用同倫法函數(shù)H(x,s)除了參數(shù)x還要輸入s和x0，這就需要牛頓迭代函數(shù)需要相應(yīng)的改變。
? ? ? ??3)NewtonIterHomo.m（專門用于同倫法的牛頓迭代函數(shù)）
?
function [x] = NewtonIterHomo(x0,Hxs,gx,s,tol,verbose)
% Version: 1.0 written by jbb0523 @2016-08-21
% 專門用于同倫方法的牛頓迭代法函數(shù)
% 目標(biāo)函數(shù)：H(x,s)=arctan(x) + (s-1)*arctan(x0)
% x0為初值, Hxs為目標(biāo)求解函數(shù), gx為Hxs的導(dǎo)數(shù)，s為目標(biāo)函數(shù)的輸入?yún)?shù)
% tol為精度要求,verbose控制是否輸出迭代過程
% x為輸出最終解
% 參考文獻(xiàn)：1)百度百科：牛頓迭代法
% 2)http://blog.sina.com.cn/s/blog_86186c970102vll8.html
% 3)http://www.maths.lth.se/na/courses/FMN081/FMN081-06/lecture8.pdf
if nargin < 6
verbose = 0;%默認(rèn)不輸出迭代過程
end
if nargin < 5
tol = 1e-3;
end
iter = 0;
x = x0;
if verbose
fprintf('Iter %d: x=%f,f(x)=%f\n',iter,x,Hxs(x,s,x0));
end
while abs(Hxs(x,s,x0)) > tol
x = x - Hxs(x,s,x0)/gx(x);%牛頓迭代公式
iter = iter + 1;
if verbose
fprintf('Iter %d: x=%f,f(x)=%f\n',iter,x,Hxs(x,s,x0));
end
end
end
?
? ? ? ??4)TestHomotopy.m（測試代碼）
?
clear all;close all;clc;
x0 = 5;
s = 0:0.1:1;
x_star = zeros(length(s),1);
neginf = -20;posinf = 20;step = 0.05;
X_Plot = neginf:step:posinf;%作圖橫坐標(biāo)
EnPlot = 1;%是否作圖,1:是，0:否
if EnPlot
figure;plot(X_Plot,Hxs(X_Plot,s(1),x0),'r');%畫出H(x,0)，即F0(x)
xlabel('x');ylabel('H(x,s,x_0)');title('同倫函數(shù)曲線')
grid on;hold on;pause;
end
x_star(1) = x0;
for i=2:1:11
x_star(i) = NewtonIterHomo(x_star(i-1),@Hxs,@gx,s(i),1e-3,1);
if EnPlot
plot(X_Plot,Hxs(X_Plot,s(i),x_star(i)),'b');
hold on;pause;
end
end
if EnPlot
hold off;
end
figure;plot(s,x_star,'.-','MarkerEdgeColor','r','MarkerSize',20);grid;
xlabel('s');ylabel('x^*');title('同倫方法零點(diǎn)收斂過程');
?
? ? ? ??運(yùn)行測試代碼，可以發(fā)現(xiàn)最后的x迭代結(jié)果為x=0.000493，得到了正確的收斂的結(jié)果。
3、結(jié)束語
? ? ? ??以前曾在文獻(xiàn)中看到過同倫法，但也沒有深究，近來在看SpaRSA算法時(shí)里面有個(gè)Continuation的概念，然后往上追蹤到了FPC算法，在FPC的文獻(xiàn)里又一次見到了Homotopymethod這個(gè)概念，看來是繞不過去了，那就花點(diǎn)時(shí)間吃掉它吧……
? ? ? ??同倫法（homotopy method）有幾個(gè)別名【2】：
? ? ? ??注意括號里面的continuation method，那么continuation method翻譯為中文應(yīng)該是什么呢？查到的中文文獻(xiàn)不多，在文獻(xiàn)【5】的3.6節(jié)中將continuation翻譯為“連續(xù)”，而在文獻(xiàn)【6】中將continuation翻譯為“持續(xù)”，有道詞典是這樣翻譯的：
? ? ? ??另外，發(fā)現(xiàn)有些老外的PPT做的很好，看起來很清晰，比如這次的文獻(xiàn)【2】，這已經(jīng)不是第一次受益于老外講課的PPT了，同樣一個(gè)概念，不同的講法，但會(huì)產(chǎn)生截然不同的效果，能不能深入淺出的把一個(gè)復(fù)雜的概念講明白是一項(xiàng)很重要的本事……
4、參考文獻(xiàn)
【1】鄧軍. 基于凸優(yōu)化的壓縮感知信號恢復(fù)算法研究[D]. 哈爾濱工業(yè)大學(xué), 2011.
【2】http://www.maths.lth.se/na/courses/FMN081/FMN081-06/lecture8.pdf
【3】百度百科，牛頓迭代法
【4】http://blog.csdn.net/itplus/article/details/21896453
【5】陳飛. 壓縮感知凸優(yōu)化方法分析[D].浙江大學(xué), 2012.
【6】楊真真, 楊震.信號壓縮與重構(gòu)的交替方向外點(diǎn)持續(xù)法[J]. 電子學(xué)報(bào),2014(3):485-490.
總結(jié)
- 上一篇： Golang基础学习总结
- 下一篇： ARM的FP寄存器