密に構造化されたヘッシアンと線形等式を使用した最小化
メモリ使用量を低減するヘッセ乗算関数
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
オプションを参照してください。目的関数は以下の構造になります。
ここで は 1000 行 2 列の行列です。 のヘッシアンは密ですが、 のヘッシアンはスパースです。 のヘッシアンが の場合、 のヘッシアン は以下のようになります。
直接 を使うことで生じる過度のメモリ使用を避けるために、この例ではヘッセ乗算関数 hmfleq1
を提供します。この関数は行列 Y
が渡されると、次のヘッセ行列の積を計算するために、 に対応するスパース行列 Hinfo
と V
を使用します。
W = H*Y = (Hinfo - V*V')*Y.
この例でヘッセ乗算関数は、ヘッセ行列の積を計算するために と V
を必要とします。V
は定数であるので V
を無名関数の関数ハンドルで扱うことができます。
しかし、 は定数ではありません。そのため、現在の x
で計算しなければなりません。これは目的関数の中で を計算することによって実行され、第 3 出力引数 Hinfo
として を返します。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.
最適化が進むと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
はアルゴリズムで決められた対角行列をベースに前提条件子を計算します。通常、これはうまく実行されません。