ドキュメンテーション

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

密に構造化されたヘッシアンと線形等式を使用した最小化

メモリ使用量を低減するヘッセ乗算関数

fmincon interior-pointtrust-region-reflective のアルゴリズム、および fminunc trust-region アルゴリズムは、ヘッシアンが密で構造化されている問題を解くことができます。これらの問題に対して、ヘッシアン H を作成するとメモリに負荷がかかるため fminconfminunc はヘッシアン H を直接使って H*Y を計算しません。その代わり、W = H*Y の計算では関数を使って H の情報と行列 Y を fmincon または fminunc に与えなければなりません。

この例は目的関数に非線形等式および線形等式があるため、fmincon を使います。この説明は信頼領域 Reflective 法にも該当し、fminunc trust-region アルゴリズムも同様です。内点法アルゴリズムについては ヘッセ乗算関数HessianMultiplyFcn を参照してください。目的関数は以下の構造になります。

f(x)=f^(x)12xTVVTx,

ここで V は 1000 行 2 列の行列です。f のヘッシアンは密ですが、f^ のヘッシアンはスパースです。 のヘッシアンが H^ の場合、f のヘッシアン H は以下のようになります。

H=H^VVT.

直接 H を使うことで生じる過度のメモリ使用を避けるために、この例ではヘッセ乗算関数 hmfleq1 を提供します。この関数は行列 Y が渡されると、次のヘッセ行列の積を計算するために、 に対応するスパース行列 HinfoV を使用します。

W = H*Y = (Hinfo - V*V')*Y

この例でヘッセ乗算関数は、ヘッセ行列積を計算するために V を必要とします。V は定数であるので V を無名関数の関数ハンドルで扱うことができます。

しかし、 は定数ではありません。そのため、現在の x で計算しなければなりません。これは目的関数の中で を計算することによって実行され、第 3 出力引数 Hinfo として を返します。optimoptions を使って 'Hessian' オプションを 'on' に設定することによって、fmincon は目的関数から Hinfo を得ることができます。そしてそれをヘッセ乗算関数の hmfleq1 に渡します。

手順 1: 目的関数、勾配およびヘッシアンのスパース部分を計算するファイル brownvv.m を記述する

例では目的関数として brownvvfmincon に渡します。ファイル brownvv.m は長すぎるため、ここでは記述しません。次のコマンドによりコードを確認することができます。

type brownvv

brownvv で勾配と目的関数を計算するため、例 (手順 3) では optimoptions を使って、SpecifyObjectiveGradient オプションを true に設定します。

手順 2: H と与えられる行列 Y の積に対するヘッセ行列を計算する関数を記述する

ここで、brownvvV で計算される Hinfo を使用する関数 hmfleq1 を定義します。これは W = H*Y = (Hinfo - V*V')*Y となるヘッセ行列の積である W を計算するための無名関数の関数ハンドルで取得することができます。この関数は、以下の形式でなければなりません。

W = hmfleq1(Hinfo,Y)

1 番目の引数は目的関数 brownvv で返される 3 番目の引数と同じでなければなりません。ヘッセ乗算関数に対する 2 番目の引数は、行列 (W = H*Y) の Y です。

なぜなら fmincon では 2 番目の引数としてヘッセ行列の積の構成に使われる Y を想定しているため、問題の次元数を n とした場合、Y は常に n 行の行列になります。Y の列数は可変です。最後に、V を扱うために無名関数の関数ハンドルを使用するため、V を 'hmfleqq' の 3 番目の引数にします。

function W = hmfleq1(Hinfo,Y,V);
%HMFLEQ1 Hessian-matrix product function for BROWNVV objective.
%   W = hmfleq1(Hinfo,Y,V) computes W = (Hinfo-V*V')*Y
%   where Hinfo is a sparse matrix computed by BROWNVV 
%   and V is a 2 column matrix.
W = Hinfo*Y - V*(V'*Y);

    メモ:   関数 hmfleq1hmfleq1.m として optimdemos フォルダーで使用できます。

手順 3: 開始値と線形な等式制約を使った非線形最小化ルーチンを呼び出す

optimdemos フォルダーにある fleq1.mat より、問題のパラメーター V と等式制約行列であるスパース行列 Aeq および beq を読み込みます。optimoptions を使って SpecifyObjectiveGradient オプションを true に設定し、HessianMultiplyFcn オプションに関数ハンドルで hmfleq1 を指定します。目的関数 brownvv と付加的なパラメーター V を使って fmincon を呼び出します。

function [fval,exitflag,output,x] = runfleq1
% RUNFLEQ1 demonstrates 'HessMult' option for FMINCON with linear
% equalities.

problem = load('fleq1');   % Get V, Aeq, beq
V = problem.V; Aeq = problem.Aeq; beq = problem.beq;
n = 1000;             % problem dimension
xstart = -ones(n,1); xstart(2:2:n,1) = ones(length(2:2:n),1); % starting point
options = optimoptions(@fmincon,...
    'Algorithm','trust-region-reflective',...
    'SpecifyObjectiveGradient',true, ...
    'HessianMultiplyFcn',@(Hinfo,Y)hmfleq1(Hinfo,Y,V),...
    'Display','iter',...
    'OptimalityTolerance',1e-9,...
    'FunctionTolerance',1e-9);
[x,fval,exitflag,output] = fmincon(@(x)brownvv(x,V),xstart,[],[],Aeq,beq,[],[], ...
                                    [],options);

上記のコードを実行するには、次のように入力します。

[fval,exitflag,output,x] = runfleq1;

optimoptions を使用して反復表示が設定されているので、このコマンドでは次の反復表示が生成されます。

                                Norm of      First-order 
 Iteration        f(x)          step          optimality   CG-iterations
     0            2297.63                      1.41e+03                
     1            1081.92        7.03821            555           2
     2            474.935        8.17769            223           2
     3            116.236        11.3752           72.7           2
     4            44.5572        7.44815           23.5           2
     5            44.5572            100           23.5           3
     6            44.5572             25           23.5           0
     7            12.5815           6.25             40           0
     8           -94.8062           6.25             41           1
     9           -510.124           12.5            318           1
    10           -510.124           12.5            318           3
    11           -571.121          3.125           66.6           0
    12           -581.933        0.78125             99           4
    13           -621.295         1.5625           75.6           5
    14            -746.16          3.125            102           3
    15           -802.344        5.15367           33.5           3
    16           -822.362        1.39916           5.37           2
    17           -823.207       0.353415           1.06           2
    18           -823.245       0.124695          0.287           3
    19           -823.246      0.0278802          0.109           5
    20           -823.246     0.00337474        0.00122           7
    21           -823.246    0.000131293       7.67e-05           6

Local minimum possible.

fmincon stopped because the final change in function value relative to 
its initial value is less than the selected value of the function tolerance.

最適化が進むとPCG の反復コストが適度に増加するこのサイズの問題に対しては、収束が速くなります。等式制約の可能領域は、下記の解において維持されます。

problem = load('fleq1');   % Get V, Aeq, beq
V = problem.V; Aeq = problem.Aeq; beq = problem.beq;
norm(Aeq*x-beq,inf)

ans =
   2.6645e-14

前処理

この例では H が陰的に存在するため fmincon は前提条件子の計算に H を使用することができません。H の代わりに fmincon は前提条件子の計算に brownvv から返される 3 番目の引数の Hinfo を使用します。HinfoH とサイズが同じであり、また H とある程度近似しているため適切な選択となります。HinfoH と同じサイズでなかった場合、fmincon はアルゴリズムで決められた対角行列をベースに前提条件子を計算します。通常、これはうまく実行されません。

この情報は役に立ちましたか?