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

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

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

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

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

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

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

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

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

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

type brownvv

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

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

ここで、brownvvとVで計算される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);

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

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

函数(fval、exitflag、输出,x) = runfleq1 %FLEQ1 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 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 = 1.8874e-14

前処理

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