次型の最小化

ドキュメント内 optimize_print.book (Page 79-92)

lsqnonlin を使用する Simulink の例題

ヘッセ行列が密であるが構造化されている 2 次型の最小化

quadprog の大規模法では、ヘッセ行列が密であるが構造化されている大規模問

題を解くことができます。 これらの問題に対して、quadprogは、中規模問題やス パース行列H を使った大規模問題で行うようにヘッセ行列 H から直接H*Y を計 算するという方法は行いません。 H の作成には多大なメモリを要するからです。

その代わり、行列 Y と H についての情報をもとに W = H*Y を計算する関数を quadprogに与えなければなりません。

この例題では、ヘッセ行列 Hは、H = B + A*A' という構造をしています。 ここ

で、B は 512 行 512 列の対称スパース行列で、A は、密な列から構成される 512 行

10 列のスパース行列です。H が密であるため、直接 H を使うことで起こる過度の

メモリ使用を避けるために、ヘッセ行列を乗算する関数 qpbox4mult を与えます。

この関数は、行列 Y を渡されると、スパース行列 A および B を使って、ヘッセ行 列の乗算 W = H*Y = (B + A*AO)*Y を計算します。

この例題では、ヘッセ行列を乗算する関数qpbox4multに引き渡される行列Aと

B が必要です。quadprogに対する第一引数として、ヘッセ行列を乗算する関数と

して渡される1つに行列を渡すことができます。 第2 の行列の値を与えるために 入れ子関数を使用できます。

ステップ 1: 1番目の引数として quadprog に渡す H の一部を決定

A またはB を、quadprog の1番目の引数として渡すことができます。 例題では、

より良い前処理の結果から ( 2-66ページの"前処理"を参照 )、1番目の引数とし

て B を渡すことにします。

quadprog(B,f,[],[],[],[],l,u,xstart,options)

ステップ 2: H に対する ヘッセ行列を計算する関数の記述 ここで、つぎのような関数 rungpbox4t を定義します。

• A および Bを用いて、W = H*Y = (B + A*AO)*Y となるヘッセ 行列の積 W を計

算する関数 qpbox4mult を定義します。 この関数は、以下の形式でなければな りません。

W = qpbox4mult(Hinfo,Y,...)

最初の2つの引数 Hinfo と Y が必要です。

• qpbox4.matより、問題のパラメータをロードします。

• optimset を使って、HessMult オプションに 関数ハンドルで qpbox4mult を指 定します。

• 第1引数を Bして quadprog をコールします。

入れ子関数qpbox4mult に対する第1引数は、quadprogに対して渡される第1引 数と同じである必要があります。 個の場合、行列 Bです。

qpbox4mult に対する第2引数は、行列 Y (of W = H*Y)です。quadprogでは、ヘッ セ行列の積の構成に使われるY を想定しているため、問題の次元数をn としたと きに、Y は常にn行の行列です。Y の列数は可変です。 関数 qpbox4mult は、行列 Aの値が外側の関数から来るように入れ子になります。

function [fval, exitflag, output, x] = runqpbox4

% RUNQPBOX4 は、範囲を伴うQUADPROGに対し、

% 'HessMult' オプションのデモンストレーションを行います。

% Copyright 1984-2004 The MathWorks, Inc.

% $Revision: 1.5.6.1 $ $Date: 2004/02/11 14:43:46 $ problem = load('qpbox4'); % xstart, u, l, B, A, fの取得 xstart = problem.xstart; u = problem.u; l = problem.l;

B = problem.B; A = problem.A; f = problem.f;

mtxmpy = @qpbox4mult;

% 関数 qpbox4mult を 関数ハンドルとして設定

% サブ関数

% HessMult オプションの選択

options = optimset('HessMult',mtxmpy);

% H 引数を通して、qpbox4mult に B を渡します。また、B は、

% PCGのpreconditionerの計算に使用されます。

% A は、'options'の後に、追加の引数として渡されます。

[x, fval, exitflag, output] = ...

quadprog(B,f,[],[],[],[],l,u,xstart,options);

function W = qpbox4mult(B,Y);

大規模例題

%QPBOX4MULT 構造化された ヘッセ行列の積 %ヘッセ行列 % W = qpbox4mult(B,Y) は、W = (B + A*A')*Y を計算します。

% INPUT:

% B - スパース正方行列 (512 × 512)

% Y - B + A'*A と乗算するベクトル (または行列)

% 外側の関数 runqpbox4からの変数 % A - 512行10列のスパース行列 %

% OUTPUT:

% W - (B + A*A')*Y の積 %

% A*A'の形成を避けるため、乗算の順番を変えています。

W = B*Y + A*(A'*Y);

end end

ステップ 3: 初期値を使った

2

次型の最小化ルーチンの呼び出し

runqpbox4内の2次型の最小化ルーチンを呼び出すためには、つぎのように入力 して、

[fval,exitflag,output] = runqpbox4

前述のコードを実行し、fval, exitflag, outputの値を表示します。 結果は、つ ぎのようになります。

Optimization terminated: relative function value changing by less than sqrt(OPTIONS.TolFun), no negative curvature detected in current

trust region model and the rate of progress (change in f(x)) is slow.

fval =

-1.0538e+003

exitflag = 3

output =

iterations: 18

algorithm: 'large-scale: reflective trust-region' firstorderopt: 0.0028

cgiterations: 50

message: [1x206 char]

30 の PCG 繰り返しによる 18 回の繰り返しの後、関数値が減ります。

fval =

-1.0538e+003

そして、1次最適性は、以下のようになります。

output.firstorderopt = 0.0043

前処理

この例では、H が陰的にのみ存在するため、quadprog は前処理演算にH を使用す ることはできません。 その代わり、quadprog は H の代わりに渡す引数として、前 処理演算で B を使用します。B はH と同じサイズであり、また H と同じ程度の近 似となるため良い選択となります。B が H と同じサイズでない場合、quadprog は アルゴリズムで決定された対角スケール行列をベースに前処理演算を行います。

通常、これはうまく実行されません。

H が陽的に使用可能な場合、前処理演算子はより良い近似になるため、TolPcg パ

ラメータの調整はやや小さめの値にする必要があります。 この例題では、前の例 題と同じですが、TolPcg はデフォルトの 0.1 から 0.01 に減らしています。

function [fval, exitflag, output, x] = runqpbox4prec

% RUNQPBOX4PREC は、境界をもつQUADPROG に対する

% 'HessMult'オプションのデモンストレーションを行います。

% Copyright 1984-2004 The MathWorks, Inc.

% $Revision: 1.5.6.1 $ $Date: 2004/02/11 14:43:46 $ problem = load('qpbox4'); % xstart, u, l, B, A, fの取得 xstart = problem.xstart; u = problem.u; l = problem.l;

B = problem.B; A = problem.A; f = problem.f;

大規模例題

mtxmpy = @qpbox4mult;

% 関数 qpbox4mult を 関数ハンドルとして設定

% HessMult オプションの選択

% TolPCG オプションの上書き

options = optimset('HessMult',mtxmpy,'TolPcg',0.01);

% 引数 H によりqpbox4mult に B を渡します。

% また、B はPCGに対するpreconditionerを計算する場合に使用されます。

% A は'options'の後に追加の引数として渡されます。

[x, fval, exitflag, output] =

quadprog(B,f,[],[],[],[],l,u,xstart,options);

function W = qpbox4mult(B,Y);

%QPBOX4MULT 構造化された ヘッセ行列の積 %ヘッセ行列 % W = qpbox4mult(B,Y) は、W = (B + A*A')*Y を計算します。

% 入力:

% B - スパース正方行列 (512 × 512)

% Y - B + A'*A と乗算するベクトル (または行列)

% 外側の関数 runqpbox4からの変数 % A - 512行10列のスパース行列 %

% 出力:

% W - (B + A*A')*Y の積 %

% A*A'の形成を避けるため、 乗算の順番を変えています。

W = B*Y + A*(A'*Y);

end end

ここで、つぎのように入力して、

[fval,exitflag,output] = runqpbox4prec

以前のコードを実行します。18 回の繰り返しと 50回の PCG 繰り返しの後、関数 値は5桁の有効数字となる同じ値になります。

fval = -1.0538e+003

しかし、1次の最適性は、さらに減ります。

output.firstorderopt = 0.0028

注意 TolPcg を小さくしすぎると、PCG 繰り返し数がかなり増えることがありま す。

範囲による制約を課した線形最小二乗法

多くの状況において、変数に範囲をもたせたスパース線形最小二乗問題が存在し

ます。 つぎの問題では、変数が非負である必要があります。 この問題は、区分的な

線形スプラインに対して関数近似を適合させるものです。 特に、粒子が単位長さ の正方形上で散乱しているものを考えます。 近似される関数はこれらの点で評価 され、区分的な線形スプライン近似が、係数は非負である条件のもとで実行され

ます。 400 個の変数にフィットする 2000 の方程式があります。

load particle % C, d の取得 lb = zeros(400,1);

[x,resnorm,residual,exitflag,output] = ...

lsqlin(C,d,[],[],[],[],lb);

デフォルトの対角型の前処理演算子が、非常に良く機能します。

exitflag = 1 resnorm = 22.5794 output =

algorithm: 'large-scale: trust-region reflective Newton' firstorderopt: 2.7870e-005

iterations: 10 cgiterations: 42

範囲による制約問題について、1次の最適性は、v.*g の無限ノルムです。 ここで、

v は 4-7ページの"ボックス制約" にあるように定義され、g は勾配です。

PrecondBandWidth をinf に設定し、各繰り返しでスパースQR 分解を使用する ことにより、1次の最適化を改善する ( 減らす ) ことができます。

options = optimset('PrecondBandWidth',inf);

大規模例題

[x,resnorm,residual,exitflag,output] = ...

lsqlin(C,d,[],[],[],[],lb,[],[],options);

繰り返し数および1次の最適性は共に減少します。

exitflag = 1 resnorm = 22.5794 output =

algorithm: 'large-scale:trust-region reflective Newton' firstorderopt: 5.5907e-015

iterations: 12 cgiterations: 11

等式および不等式を使った線形計画法

問題は、つぎのとおりです。

つ ぎ の よ う に し て、行 列 お よ び ベ ク ト ル A, Aeq, b, beq, f お よ び 下 限 lb を MATLAB ワークスペースにロードすることができます。

load sc50b

sc50b.mat の問題には、48 の変数、30 の不等式、および 20 の等式があります。

linprogを使用して、この問題を解くことができます。

[x,fval,exitflag,output] = ...

linprog(f,A,b,Aeq,beq,lb,[],[],optimset('Display','iter'));

optimset を使用して、繰り返しでの表示を設定しているので、表示される結果

は、つぎのようになります。

Residuals: Primal Dual Duality Total Infeas Infeas Gap Rel A*x-b A'*y+z-f x'*z Error Iter 0: 1.50e+003 2.19e+001 1.91e+004 1.00e+002 Iter 1: 1.15e+002 2.94e-015 3.62e+003 9.90e-001

fTx

min such that

Aeq x⋅ = beq A x⋅ ≤b

x≥0

Iter 2: 1.16e-012 2.21e-015 4.32e+002 9.48e-001 Iter 3: 3.23e-012 5.16e-015 7.78e+001 6.88e-001 Iter 4: 5.78e-011 7.61e-016 2.38e+001 2.69e-001 Iter 5: 9.31e-011 1.84e-015 5.05e+000 6.89e-002 Iter 6: 2.96e-011 1.62e-016 1.64e-001 2.34e-003 Iter 7: 1.51e-011 2.74e-016 1.09e-005 1.55e-007 Iter 8: 1.51e-012 2.37e-016 1.09e-011 1.51e-013 Optimization terminated successfully.

この問題では、大規模線形計画法アルゴリズムが、スケーリングされた残差をデ フォルトの許容範囲 1e-08 以下まで速やかに減少させることがわかります。

exitflag が正であることは、ユーザの linprog が収束したことを意味します。 ま た、fvalの最終関数値および output.iterations の繰り返し数を得ることもで きます。

exitflag = 1 fval = -70.0000 output =

iterations: 8 cgiterations: 0

algorithm: 'lipsol'

等式の中の密の列をもつ線形計画法

問題は、つぎのとおりです。

つぎのようにして、行列およびベクトル Aeq, beq, f, lb, および ub を MATLAB ワークスペースにロードすることができます。

load densecolumns

densecolumns.mat の問題は1677 の変数と627の等式からなり、すべての変数が 下限をもち、399 の変数が上限をもちます。 等式行列 Aeq は最初の 25 列が密な列

であり、spy プロットで簡単に確認することができます。

spy(Aeq) fTx

min such that Aeq x⋅ = beq lb≤ ≤x ub

大規模例題

linprogを使用して、この問題を解くことができます。

[x,fval,exitflag,output] = ...

linprog(f,[],[],Aeq,beq,lb,ub,[],optimset('Display','iter'));

optimsetを使用して繰り返しでの表示を設定しているので、表示される結果は、

つぎのようになります。

Residuals: Primal Dual Upper Duality Total Infeas Infeas Bounds Gap Rel A*x-b A'*y+z-w-f {x}+s-ub x'*z+s'*w Error ---Iter 0: 1.67e+003 8.11e+002 1.35e+003 5.30e+006 2.92e+001 Iter 1: 1.37e+002 1.33e+002 1.11e+002 1.27e+006 2.48e+000 Iter 2: 3.56e+001 2.38e+001 2.89e+001 3.42e+005 1.99e+000 Iter 3: 4.86e+000 8.88e+000 3.94e+000 1.40e+005 1.89e+000 Iter 4: 4.24e-001 5.89e-001 3.44e-001 1.91e+004 8.41e-001 Iter 5: 1.23e-001 2.02e-001 9.97e-002 8.41e+003 5.79e-001 Iter 6: 3.98e-002 7.91e-002 3.23e-002 4.05e+003 3.52e-001 Iter 7: 7.25e-003 3.83e-002 5.88e-003 1.85e+003 1.85e-001 Iter 8: 1.47e-003 1.34e-002 1.19e-003 8.12e+002 8.52e-002 Iter 9: 2.52e-004 3.39e-003 2.04e-004 2.78e+002 2.99e-002 Iter 10: 3.46e-005 1.08e-003 2.81e-005 1.09e+002 1.18e-002 Iter 11: 6.95e-007 1.53e-012 5.64e-007 1.48e+001 1.62e-003 Iter 12: 1.04e-006 2.26e-012 3.18e-008 8.32e-001 9.09e-005 Iter 13: 3.08e-006 1.23e-012 3.86e-009 7.26e-002 7.94e-006 Iter 14: 3.75e-007 1.09e-012 6.53e-012 1.11e-003 1.21e-007 Iter 15: 5.21e-008 1.30e-012 3.27e-013 8.62e-008 9.15e-010 Optimization terminated successfully.

exitflag, fval および output の戻り値を確認することができます。

exitflag = 1 fval =

9.1464e+003 output =

iterations: 15 cgiterations: 225 algorithm: 'lipsol'

この例題において、Aeq の中の列が密な列であることがわかっているので、PCG 繰り返し数 ( output.cgiterations ) は、非零になります。 スパース Cholesky 分 解 を 使 用 す る 代 わ り に、linprog は、Sherman-Morrison の 公 式 を 使 用 し て Aeq*Aeq'を含む線形システムを解こうと試みます。 Sherman-Morrison の公式が十 分な残差を出さない場合、PCG 繰り返しを使用します。 4-13ページの"大規模線

形計画法"の節の "メインアルゴリズム"を参照ください。

ドキュメント内 optimize_print.book (Page 79-92)