Main Content

密な構造化されたヘッシアンを使った二次最小化

quadprogtrust-region-reflective 法では、ヘッシアンが密であっても構造化されていれば、大規模な問題を解くことができます。これらの問題に対し quadprog は、スパース行列 H を使った信頼領域 Reflective 法の問題で行うように、ヘッシアン H を使用して直接 H*Y を計算するという方法は取りません。H の形成には多大なメモリを要するからです。その代わり、行列 YH の情報をもとに W = H*Y を計算する関数を quadprog に与えなければなりません。

この例ではヘッセ行列 HH = B + A*A' という構造をしています。ここで、B は 512 行 512 列の対称スパース行列で A は密な列から構成される 512 行 10 列のスパース行列です。H が密であるため、直接 H を使うことで起こる過度のメモリ使用を避けるために、この例ではヘッセ乗算関数 qpbox4mult を与えます。この関数は行列 Y が渡されると、スパース行列 AB を使用してヘッセ行列の積 W = H*Y = (B + A*A')*Y を計算します。

この例の最初の部分では、ヘッセ乗算関数 qpbox4mult に、行列 AB を渡す必要があります。quadprog の第 1 引数に、ヘッセ乗算関数に渡す行列を入力できます。第 2 の行列の値を与えるために入れ子関数を使用できます。

この例の 2 番目の部分では、正確な行列 H ではなく近似の前提条件子を使用することによる影響を補正するために TolPCG 許容誤差を厳しくする方法を説明します。

手順 1: 1 番目の引数として quadprog に渡す H の一部を決定する

A または Bquadprog の 1 番目の引数として渡すことができます。例では良い前提条件子 (前処理を参照) の結果から 1 番目の引数として B を渡すことにします。

手順 2: H に対するヘッセ行列の乗算を計算する関数を記述する

ここで、次のような関数 runqpbox4 を定義します。

  • A および B を用いて、W = H*Y = (B + A*A')*Y となるヘッセ行列の積 W を計算する入れ子関数 qpbox4mult を含みます。この入れ子関数は、以下の形式でなければなりません。ここで、最初の 2 つの引数 HinfoY が必要です。

W = qpbox4mult(Hinfo,Y,...)
  • qpbox4.mat から問題のパラメーターを読み込みます。

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

  • B を第 1 引数として quadprog を呼び出します。

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

qpbox4mult に対する第 2 引数は行列 Y (W = H*Y のもの) です。quadprog ではヘッセ行列の積の構成に使われる Y を想定しているため、問題の次元数を n とした場合、Y は常に n 行の行列になります。Y の列数は可変です。関数 qpbox4mult は、行列 A の値が外側の関数から得られる入れ子になります。runqpbox4.m ファイルは、この例の終わりに掲載しています。

手順 3: 開始点を使った二次の最小化ルーチンを呼び出す

runqpbox4 内の二次の最小化ルーチンを呼び出すには次のように入力してください。

[fval,exitflag,output] = runqpbox4;
Local minimum possible.

quadprog stopped because the relative change in function value is less than the sqrt of the function tolerance, the rate of change in the function value is slow, and no negative curvature was detected.

次に、fvalexitflagoutput.iterations および output.cgiterations の値を表示します。

fval,exitflag,output.iterations, output.cgiterations
fval = -1.0538e+03
exitflag = 3
ans = 18
ans = 30

30 の PCG 反復による 18 回の反復の後、関数値が減ります。

fval
fval = -1.0538e+03

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

output.firstorderopt
ans = 0.0043

前処理

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

H が陽的に使用可能な場合、前提条件子はより良い近似になるため、TolPCG パラメーターの調整はやや小さめの値にしなければなりません。この例は前の例と同じですが、TolPCG は既定の 0.1 から 0.01 に減らしています。関数 runqpbox4prec は、この例の終わりにリストしています。

[fval,exitflag,output] = runqpbox4prec;
Local minimum possible.

quadprog stopped because the relative change in function value is less than the sqrt of the function tolerance, the rate of change in the function value is slow, and no negative curvature was detected.

18 回の反復と 50 回の PCG 反復の後、関数値は 5 桁の有効数字となる同じ値になります。

fval
fval = -1.0538e+03

さらに、1 次の最適性は基本的には同じです。

output.firstorderopt
ans = 0.0028

メモ: TolPCG を小さくしすぎると、PCG 反復数がかなり増えることがあります。

補助関数

次のコードは、補助関数 runqpbox4 を作成します。

function [fval, exitflag, output, x] = runqpbox4
%RUNQPBOX4 demonstrates 'HessianMultiplyFcn' option for QUADPROG with bounds.

problem = load('qpbox4'); % Get 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; % function handle to qpbox4mult nested function

% Choose algorithm and the HessianMultiplyFcn option
options = optimoptions(@quadprog,'Algorithm','trust-region-reflective','HessianMultiplyFcn',mtxmpy);

% Pass B to qpbox4mult via the H argument. Also, B will be used in
% computing a preconditioner for PCG.
[x, fval, exitflag, output] = quadprog(B,f,[],[],[],[],l,u,xstart,options);

    function W = qpbox4mult(B,Y)
        %QPBOX4MULT Hessian matrix product with dense structured Hessian.
        %   W = qpbox4mult(B,Y) computes W = (B + A*A')*Y where
        %   INPUT:
        %       B - sparse square matrix (512 by 512)
        %       Y - vector (or matrix) to be multiplied by B + A'*A.
        %   VARIABLES from outer function runqpbox4:
        %       A - sparse matrix with 512 rows and 10 columns.
        %
        %   OUTPUT:
        %       W - The product (B + A*A')*Y.
        %

        % Order multiplies to avoid forming A*A',
        %   which is large and dense
        W = B*Y + A*(A'*Y);
    end

end

次のコードは、補助関数 runqpbox4prec を作成します。

function [fval, exitflag, output, x] = runqpbox4prec
%RUNQPBOX4PREC demonstrates 'HessianMultiplyFcn' option for QUADPROG with bounds.

problem = load('qpbox4'); % Get 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; % function handle to qpbox4mult nested function

% Choose algorithm, the HessianMultiplyFcn option, and override the TolPCG option
options = optimoptions(@quadprog,'Algorithm','trust-region-reflective',...
    'HessianMultiplyFcn',mtxmpy,'TolPCG',0.01);

% Pass B to qpbox4mult via the H argument. Also, B will be used in
% computing a preconditioner for PCG.
% A is passed as an additional argument after 'options'
[x, fval, exitflag, output] = quadprog(B,f,[],[],[],[],l,u,xstart,options);

    function W = qpbox4mult(B,Y)
        %QPBOX4MULT Hessian matrix product with dense structured Hessian.
        %   W = qpbox4mult(B,Y) computes W = (B + A*A')*Y where
        %   INPUT:
        %       B - sparse square matrix (512 by 512)
        %       Y - vector (or matrix) to be multiplied by B + A'*A.
        %   VARIABLES from outer function runqpbox4prec:
        %       A - sparse matrix with 512 rows and 10 columns.
        %
        %   OUTPUT:
        %       W - The product (B + A*A')*Y.
        %

        % Order multiplies to avoid forming A*A',
        %   which is large and dense
        W = B*Y + A*(A'*Y);
    end

end

参考

関連するトピック