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 を想定しているため、問題の次元数をｎ としたと きに、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ページの"大規模線
形計画法"の節の "メインアルゴリズム"を参照ください。