ドキュメンテーション

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

密に構造化されたヘッセ行列と線形等式を使用した最小化

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

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 アルゴリズムも同様です。内点法アルゴリズムについては ヘッセ行列'HessMult' を参照してください。目的関数は以下の構造になります。

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.mbrownvv.m は長すぎるため、ここでは記述しません。次のコマンドによりコードを確認することができます。

type brownvv

brownvv は勾配およびヘッセ行列の一部を目的関数と同様に計算するため、例 (手順 3) では optimoptions を使用して GradObj および Hessian オプションを 'on' に設定します。

手順 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 を使って GradObj および Hessian オプションを 'on' に設定し、HessMult オプションに関数ハンドルで 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','GradObj','on', ...
    'Hessian','user-supplied','HessMult',@(Hinfo,Y)hmfleq1(Hinfo,Y,V),'Display','iter', ...
    'TolFun',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            1084.59         6.3903            578           1
     2            1084.59            100            578           3
     3            1084.59             25            578           0
     4            1084.59           6.25            578           0
     5            1047.61         1.5625            240           0
     6            761.592          3.125           62.4           2
     7            761.592           6.25           62.4           4
     8            746.478         1.5625            163           0
     9            546.578          3.125           84.1           2
    10            274.311           6.25           26.9           2
    11            55.6193        11.6597             40           2
    12            55.6193             25             40           3
    13            22.2964           6.25           26.3           0
    14            -49.516           6.25             78           1
    15           -93.2772         1.5625             68           1
    16           -207.204          3.125           86.5           1
    17           -434.162           6.25           70.7           1
    18           -681.359           6.25           43.7           2
    19           -681.359           6.25           43.7           4
    20           -698.041         1.5625            191           0
    21           -723.959          3.125            256           7
    22            -751.33        0.78125            154           3
    23           -793.974         1.5625           24.4           3
    24           -820.831        2.51937           6.11           3
    25           -823.069       0.562132           2.87           3
    26           -823.237       0.196753          0.486           3
    27           -823.245      0.0621202          0.386           3
    28           -823.246      0.0199951           0.11           6
    29           -823.246     0.00731333         0.0404           7
    30           -823.246     0.00505883         0.0185           8
    31           -823.246     0.00126471        0.00268           9
    32           -823.246     0.00149326        0.00521           9
    33           -823.246    0.000373314        0.00091           9

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.3093e-14

前処理

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

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