密に構造化されたヘッシアンと線形等式を使用した最小化
R2026bメモリ使用量を低減するヘッセ乗算関数
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 番目の引数にします。この手法については、追加パラメーターの受け渡しで詳しく説明します。
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); end
手順 3: 開始点と線形な等式制約を使った非線形最小化ルーチンを呼び出す
この例を実行する際に使用できる fleq1.mat より、問題のパラメーター V と等式制約行列であるスパース行列 Aeq および beq を読み込みます。optimoptions を使って SpecifyObjectiveGradient オプションを true に設定し、HessianMultiplyFcn オプションに関数ハンドルで hmfleq1 を指定します。目的関数 brownvv と付加的なパラメーター V を使って fmincon を呼び出します。このプロシージャを実行するコードは次のとおりです。
function [fval,exitflag,output,x] = runfleq1 % RUNFLEQ1 demonstrates the 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); end
このコードを実行するには、次のように入力します。
[fval,exitflag,output,x] = runfleq1;
Norm of First-order
Iteration f(x) step optimality CG-iterations Wall-Clock Time
0 2297.63 1.41e+03 0.110918 sec
1 1084.59 6.3903 578 1 0.209022 sec
2 1084.59 100 578 3 0.248810 sec
3 1084.59 25 578 0 0.262249 sec
4 1084.59 6.25 578 0 0.303379 sec
5 1047.61 1.5625 240 0 0.333068 sec
6 761.592 3.125 62.4 2 0.378456 sec
7 761.592 6.25 62.4 4 0.417943 sec
8 746.478 1.5625 163 0 0.425663 sec
9 546.578 3.125 84.1 2 0.480939 sec
10 274.311 6.25 26.9 2 0.507347 sec
11 55.6193 11.6597 40 2 0.529394 sec
12 55.6193 25 40 3 0.552373 sec
13 22.2964 6.25 26.3 0 0.557550 sec
14 -49.516 6.25 78 1 0.581471 sec
15 -93.2772 1.5625 68 1 0.623601 sec
16 -207.204 3.125 86.5 1 0.652120 sec
17 -434.162 6.25 70.7 1 0.675557 sec
18 -681.359 6.25 43.7 2 0.694861 sec
19 -681.359 6.25 43.7 4 0.715093 sec
20 -698.041 1.5625 191 0 0.720264 sec
21 -723.959 3.125 256 7 0.736098 sec
22 -751.33 0.78125 154 3 0.751555 sec
23 -793.974 1.5625 24.4 3 0.765335 sec
24 -820.831 2.51937 6.11 3 0.781010 sec
25 -823.069 0.562132 2.87 3 0.800709 sec
26 -823.237 0.196753 0.486 3 0.814352 sec
27 -823.245 0.0621202 0.386 3 0.829405 sec
28 -823.246 0.0199951 0.11 6 0.873742 sec
29 -823.246 0.00731333 0.0404 7 0.910990 sec
30 -823.246 0.00505883 0.0185 8 0.937816 sec
31 -823.246 0.00126471 0.00268 9 0.960486 sec
32 -823.246 0.00149326 0.00521 9 0.979944 sec
33 -823.246 0.000373314 0.00091 9 1.001735 sec
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 はアルゴリズムで決められた対角行列をベースに前提条件子を計算します。通常、これはうまく実行されません。
