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

ライブ スクリプトを開く

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

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

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

f(x)=fˆ(x)-12xTVVTx,

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

H=Hˆ-VVT.

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

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

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

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

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

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

type brownvv

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

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

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

W = hmfleq1(Hinfo,Y)

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

なぜなら fmincon では 2 番目の引数としてヘッセ行列の積の構成に使われる Y を想定しているため、問題の次元数を n とした場合、Y は常に n 行の行列になります。Y の列数は可変です。最後に、V を扱うために無名関数の関数ハンドルを使用するため、V を hmfleq1 の 3 番目の引数にします。この手法については、追加パラメーターの受け渡しで詳しく説明します。以下は、ファイル hmfleq1.m のリストです。

type hmfleq1.m
function W = hmfleq1(Hinfo,Y,V)
%HMFLEQ1 Hessian-matrix product function for BROWNVV objective.
% Documentation example.
% W = hmfbx4(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.

%   Copyright 1984-2008 The MathWorks, Inc.
  
W = Hinfo*Y - V*(V'*Y);

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

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

type runfleq1.m
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;
                                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 value of the function tolerance.

<stopping criteria details>

最適化が進むと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 を使用します。Hinfo は H とサイズが同じであり、また H とある程度近似しているため適切な選択となります。Hinfo が H と同じサイズでない場合、fmincon はアルゴリズムで決められた対角行列をベースに前提条件子を計算します。通常、これはうまく実行されません。

参考

|

トピック