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