minimize.m:共軛梯度法更新BP算法權值

時間 2019-12-28

標籤 minimize.m minimize 共軛梯度更新算法權值简体版

原文原文鏈接

minimize.m:共軛梯度法更新BP算法權值

做者：凱魯嘎吉 - 博客園 http://www.cnblogs.com/kailugaji/html

Carl Edward Rasmussen在高斯機器學習的MATLAB代碼中寫到一個優化類的函數：minimize.m，同時，Geoff Hinton在用BP算法精調深度自編碼網絡時，也借鑑了這個函數minimize.m，下面來簡單聊一聊這個函數的大體機理。java

matlab函數minimum.m用來查找（非線性）多元函數的（局部）最小值。用戶必須提供一個函數，該函數返回全部變量的值和偏導數。該函數基於具備Wolfe-Powel條件的多項式插值，使用Polak-Ribiere共軛梯度和近似線性搜索。算法

做用：Minimize a differentiable multivariate function.網絡

1. 線性搜索技術——肯定迭代步長

2. 非線性共軛梯度法——肯定搜索方向

3. MATLAB代碼詳解

function [X, fX, i] = minimize(X, f, length, varargin)
%X是權值偏置 f輸出的是代價函數和偏導 3次線性搜索 每層網絡對應的節點數Dim和訓練數據data
%f是一個函數的名稱，它主要是用來計算網絡中的代價函數以及代價函數對各個參數X的偏導函數，f的參數值分別爲X，以及minimize函數後面的P1,P2,P3,…使用共軛梯度法進行優化的最大線性搜索長度爲length。
%返回值X爲3次線性搜索最優化後獲得的權值參數，是一個列向量，fX爲在此最優參數X下的代價函數，i爲線性搜索的長度（即迭代的次數）。

% Minimize a differentiable multivariate function. 
%
% Usage: [X, fX, i] = minimize(X, f, length, P1, P2, P3, ... )

%更新參數W和b，ΔW=步長*方向,Δb=步長*方向。步長用Wolfe不肯定線性搜索進行計算，而降低的方向用Polack-Ribiere共軛梯度進行計算。最終輸出更新完以後的參數W,b

%最小化連續微分多元函數。

%起點由「 X」（D乘1）給定，而且在字符串「 f」中命名的函數必須返回函數值和偏導數向量。
%共軛梯度的Polack-Ribiere風格用於計算搜索方向，而且使用二次多項式和三次多項式逼近以及Wolfe-Powell中止準則的線搜索以及斜率比方法來猜想初始步長。
%此外，還要進行大量檢查，以確保正在進行探索，而且推斷不會一望無際。 
%「length」給出了運行的長度：若是爲正，則給出最大的線性搜索次數；若是爲負，則其絕對值給出最大的函數求值次數。
%當函數的長度變長或沒法進一步進行處理時（即，咱們處於最小狀態，或因爲數值問題而接近時，咱們沒法進一步接近），該函數將返回。
%若是函數在幾回迭代中終止，則可能代表函數值和導數不一致（即，「 f」函數的實現中可能存在錯誤）。
%函數返回找到的解「 X」，函數值「 fX」的向量表示進展，「 i」使用的迭代次數（線性搜索或函數評估，取決於「length」的符號）。

%當函數的長度增長或沒法進一步處理時（即，咱們處於（局部）最小值，或因爲數值問題而接近），函數將返回。 
%注意：若是函數在幾回迭代中終止，則可能代表函數值和導數不一致（即，「 f」函數的實現中可能存在錯誤）。 
%函數返回找到的解「 X」，函數值「 fX」的向量表示進展，「 i」使用的迭代次數（行搜索或函數評估，取決於「長度」的符號）。

INT = 0.1;    % don't reevaluate within 0.1 of the limit of the current bracket不要在當前括號限制的0.1之內從新評估
EXT = 3.0;                  % extrapolate maximum 3 times the current step-size外推最大值爲當前步長的3倍
MAX = 20;                         % max 20 function evaluations per line search每次線性搜索最多20個函數求值
RATIO = 10;                                       % maximum allowed slope ratio最大容許斜率
SIG = 0.1; RHO = SIG/2; 
% SIG和RHO是控制Wolfe-Powell條件的常數。 
% SIG是先前斜率和新斜率（搜索方向上的導數）之間容許的最大絕對比率，所以將SIG設置爲低（正）值將強制線搜索中的更高精度。 
% RHO是指望值的最小容許分數（從線性搜索中起始點的斜率開始）。 
% 常數必須知足0 <RHO <SIG <1。調整SIG（取決於要優化的函數的性質）可能會加快最小化；使用rho可能不值得。

%在開始沿着最陡降低的方向進行初始行搜索以後，代碼天然分爲3部分。
%1）咱們首先進入一個while循環，它使用點1（p1）和（p2）來計算外推（p3），直到咱們外推足夠遠（wolfe-powell條件）。
%2）若有必要，咱們進入第二個循環，其中p二、p3和p4選擇包含（局部）最小值的子區間，並對其進行插值，找到一個可接受的點（wolfe-powell條件）。請注意，點始終保持順序p0<=p1<=p2<p3<p4。
%3）使用共軛梯度（polack-ribiere-flavor）計算新的搜索方向，或者在前一線性搜索中出現問題時恢復到最陡。
%若是兩個連續的線性搜索失敗，或者當函數計算或線性搜索用完時，返回迄今爲止的最佳值。
%在外推過程當中，「f」函數可能會因錯誤或返回nan或inf而失敗，minimize應該能很好地處理這個問題。
if max(size(length)) == 2
    red=length(2); 
    length=length(1); 
else %length=3
    red=1;
end
if length>0 
    S='Linesearch';  %線性搜索
else
    S='Function evaluation';  %函數求值
end

i = 0;                                            % zero the run length counter 運行長度計數器清零
ls_failed = 0;                             % no previous line search has failed先前的線性搜索沒有失敗
[f0 df0] = feval(f, X, varargin{:});          % get function value and gradient
fX = f0;
i = i + (length<0);                                            % count epochs?!
s = -df0; d0 = -s'*s;           % initial search direction (steepest) and slope初始搜索方向(最陡，負梯度方向)和斜率
x3 = red/(1-d0);                                  % initial step is red/(|s|+1) 初始步長

while i < abs(length)                                      % while not finished
  i = i + (length>0);                                      % count iterations?!

  X0 = X; F0 = f0; dF0 = df0;                   % make a copy of current values
  if length>0, M = MAX; else M = min(MAX, -length-i); end
  %用p一、p2外推p3
  while 1                             % keep extrapolating as long as necessary
    x2 = 0; f2 = f0; d2 = d0; f3 = f0; df3 = df0;
    success = 0;
    while ~success && M > 0
      try
        M = M - 1; i = i + (length<0);                         % count epochs?!
        [f3 df3] = feval(f, X+x3*s, varargin{:});    %權值(t+1)=權值(t)+初始步長*初始搜索方向
        if isnan(f3) || isinf(f3) || any(isnan(df3)+isinf(df3)), error(''), end
        success = 1;
      catch                                % catch any error which occured in f
        x3 = (x2+x3)/2;                                  % bisect and try again  %步長等分，選取新搜索點
      end
    end
    if f3 < F0, X0 = X+x3*s; F0 = f3; dF0 = df3; end         % keep best values
    d3 = df3'*s;                                                    % new slope
    if d3 > SIG*d0 || f3 > f0+x3*RHO*d0 || M == 0  % are we done extrapolating?
      break
    end
    x1 = x2; f1 = f2; d1 = d2;                        % move point 2 to point 1
    x2 = x3; f2 = f3; d2 = d3;                        % move point 3 to point 2 
    A = 6*(f1-f2)+3*(d2+d1)*(x2-x1);                 % make cubic extrapolation
    B = 3*(f2-f1)-(2*d1+d2)*(x2-x1);
    x3 = x1-d1*(x2-x1)^2/(B+sqrt(B*B-A*d1*(x2-x1))); % num. error possible, ok!
    if ~isreal(x3) || isnan(x3) || isinf(x3) || x3 < 0 % num prob | wrong sign?
      x3 = x2*EXT;                                 % extrapolate maximum amount
    elseif x3 > x2*EXT                  % new point beyond extrapolation limit?
      x3 = x2*EXT;                                 % extrapolate maximum amount
    elseif x3 < x2+INT*(x2-x1)         % new point too close to previous point?
      x3 = x2+INT*(x2-x1);
    end
  end                                                       % end extrapolation
  %插值p二、p3和p4
  while (abs(d3) > -SIG*d0 || f3 > f0+x3*RHO*d0) && M > 0  % keep interpolating
    if d3 > 0 || f3 > f0+x3*RHO*d0                         % choose subinterval
      x4 = x3; f4 = f3; d4 = d3;                      % move point 3 to point 4
    else
      x2 = x3; f2 = f3; d2 = d3;                      % move point 3 to point 2
    end
    if f4 > f0           
      x3 = x2-(0.5*d2*(x4-x2)^2)/(f4-f2-d2*(x4-x2));  % quadratic interpolation 二次插值
    else
      A = 6*(f2-f4)/(x4-x2)+3*(d4+d2);                    % cubic interpolation 三次插值
      B = 3*(f4-f2)-(2*d2+d4)*(x4-x2);
      x3 = x2+(sqrt(B*B-A*d2*(x4-x2)^2)-B)/A;        % num. error possible, ok!
    end
    if isnan(x3) || isinf(x3)
      x3 = (x2+x4)/2;               % if we had a numerical problem then bisect
    end
    x3 = max(min(x3, x4-INT*(x4-x2)),x2+INT*(x4-x2));  % don't accept too close
    [f3 df3] = feval(f, X+x3*s, varargin{:});
    if f3 < F0, X0 = X+x3*s; F0 = f3; dF0 = df3; end         % keep best values
    M = M - 1; i = i + (length<0);                             % count epochs?!
    d3 = df3'*s;                                                    % new slope
  end                                                       % end interpolation
  %用Polack-Ribiere共軛梯度法更新搜索方向
  if abs(d3) < -SIG*d0 && f3 < f0+x3*RHO*d0          % if line search succeeded
    X = X+x3*s; f0 = f3; fX = [fX' f0]';                     % update variables
    fprintf('%s %6i;  Value %4.6e\r', S, i, f0);
    s = (df3'*df3-df0'*df3)/(df0'*df0)*s - df3;   % Polack-Ribiere 共軛梯度方向  搜索方向的更新公式
    df0 = df3;                                               % swap derivatives
    d3 = d0; d0 = df0'*s;
    if d0 > 0                                      % new slope must be negative
      s = -df0; d0 = -s'*s;                  % otherwise use steepest direction 負梯度方向
    end
    x3 = x3 * min(RATIO, d3/(d0-realmin));          % slope ratio but max RATIO
    ls_failed = 0;                              % this line search did not fail
  else
    X = X0; f0 = F0; df0 = dF0;                     % restore best point so far
    if ls_failed || i > abs(length)         % line search failed twice in a row
      break;                             % or we ran out of time, so we give up
    end
    s = -df0; d0 = -s'*s;                                        % try steepest
    x3 = 1/(1-d0);                     
    ls_failed = 1;                                    % this line search failed
  end
end
fprintf('\n');