多くの変数を使用した非線形問題の求解

この例では、非線形問題で多数の変数を処理する方法を説明します。一般的に、このタイプの問題は次のようになります。

低メモリの BFGS (LBFGS) ヘッセ近似を使用します。このオプションは、既定の fminunc アルゴリズムおよび fmincon アルゴリズムで使用できます。
明示的な勾配がある場合は、有限差分ヘッシアンと 'cg' 部分問題アルゴリズムも使用できます。
明示的なヘッシアンがある場合は、そのヘッシアンをスパースとして定式化します。
この例には含まれていませんが、構造化された問題があり、ヘッシアンの積を任意のベクトルで評価できる場合は、ヘッセ乗算関数の使用を試みます。詳細は、密に構造化されたヘッシアンと線形等式を使用した最小化を参照してください。

この例では、一般的な非線形ソルバー fminunc および fmincon に補助関数 hfminunc0obj (この例の終わりに掲載) を使用します。この関数は、Rosenbrock 関数の N 次元一般化であり、数値的に最小化するのが難しい関数です。最小値 0 は一意の点 x = ones(N,1) において得られます。

この関数は、明示的な二乗和です。したがって、この例は最小二乗法ソルバー使用の効率も示します。最小二乗法ソルバー lsqnonlin の場合、この例では関数 hfminunc0obj と等価のベクトル目的関数として補助関数 hlsqnonlin0obj (この例の終わりに掲載) を使用します。

問題の設定

1,000 個の変数に対して目的関数 hfminunc0obj を使用する問題を設定します。各変数の初期点 x0 を -2 に設定します。

fun = @hfminunc0obj;
N = 1e3;
x0 = -2*ones(N,1);

初期オプションについて、関数評価または反復の回数の表示も制限も指定しません。

options = optimoptions("fminunc",Display="none",...
    MaxFunctionEvaluations=Inf,MaxIterations=Inf);

データを保持するテーブルを設定します。ソルバー実行 8 回分のラベルを指定し、実行時間を収集する場所、返される関数値、終了フラグ、反復の回数、および反復ごとの時間を定義します。

ExperimentLabels = ["BFGS_NoGrad", "LBFGS_NoGrad",...
    "BFGS_Grad", "LBFGS_Grad", "Analytic", "fin-diff-grads",...
    "LSQ_NoJacob", "LSQ_Jacob"];
timetable = table('Size',[8 5],'VariableTypes',["double" "double" "double" "double" "double"],...
    'VariableNames',["Time" "Fval" "Eflag" "Iters" "TimePerIter"],...
    'RowNames',ExperimentLabels);

次のコードセクションは、ソルバー実行ごとの適切なオプションを作成し、出力を可能な限りテーブルに直接収集します。

BFGS ヘッセ近似、勾配なし

opts = options;
opts.HessianApproximation = 'bfgs';
opts.SpecifyObjectiveGradient = false;
overallTime = tic;
[~,timetable.Fval("BFGS_NoGrad"),timetable.Eflag("BFGS_NoGrad"),output] =...
    fminunc(fun, x0, opts);
timetable.Time("BFGS_NoGrad") = toc(overallTime);
timetable.Iters("BFGS_NoGrad") = output.iterations;
timetable.TimePerIter("BFGS_NoGrad") =...
    timetable.Time("BFGS_NoGrad")/timetable.Iters("BFGS_NoGrad");

LBFGS ヘッセ近似、勾配なし

opts = options;
opts.HessianApproximation = 'lbfgs';
opts.SpecifyObjectiveGradient = false;
overallTime = tic;
[~,timetable.Fval("LBFGS_NoGrad"),timetable.Eflag("LBFGS_NoGrad"),output] =...
    fminunc(fun, x0, opts);
timetable.Time("LBFGS_NoGrad") = toc(overallTime);
timetable.Iters("LBFGS_NoGrad") = output.iterations;
timetable.TimePerIter("LBFGS_NoGrad") =...
    timetable.Time("LBFGS_NoGrad")/timetable.Iters("LBFGS_NoGrad");

勾配のある BFGS

opts = options;
opts.HessianApproximation = 'bfgs';
opts.SpecifyObjectiveGradient = true;
overallTime = tic;
[~,timetable.Fval("BFGS_Grad"),timetable.Eflag("BFGS_Grad"),output] =...
    fminunc(fun, x0, opts);
timetable.Time("BFGS_Grad") = toc(overallTime);
timetable.Iters("BFGS_Grad") = output.iterations;
timetable.TimePerIter("BFGS_Grad") =...
    timetable.Time("BFGS_Grad")/timetable.Iters("BFGS_Grad");

勾配のある LBFGS

opts = options;
opts.HessianApproximation = 'lbfgs';
opts.SpecifyObjectiveGradient = true;
overallTime = tic;
[~,timetable.Fval("LBFGS_Grad"),timetable.Eflag("LBFGS_Grad"),output] =...
    fminunc(fun, x0, opts);
timetable.Time("LBFGS_Grad") = toc(overallTime);
timetable.Iters("LBFGS_Grad") = output.iterations;
timetable.TimePerIter("LBFGS_Grad") =...
    timetable.Time("LBFGS_Grad")/timetable.Iters("LBFGS_Grad");

解析的ヘッシアン、`'trust-region'` アルゴリズム

opts = options;
opts.Algorithm = 'trust-region';
opts.SpecifyObjectiveGradient = true;
opts.HessianFcn = "objective";
overallTime = tic;
[~,timetable.Fval("Analytic"),timetable.Eflag("Analytic"),output] =...
    fminunc(fun, x0, opts);
timetable.Time("Analytic") = toc(overallTime);
timetable.Iters("Analytic") = output.iterations;
timetable.TimePerIter("Analytic") =...
    timetable.Time("Analytic")/timetable.Iters("Analytic");

勾配のある有限差分ヘッシアン、`fmincon` ソルバー

opts = optimoptions("fmincon","SpecifyObjectiveGradient",true,...
    "Display","none","HessianApproximation","finite-difference",...
    SubproblemAlgorithm="cg",MaxFunctionEvaluations=Inf,MaxIterations=Inf);
overallTime = tic;
[~,timetable.Fval("fin-diff-grads"),timetable.Eflag("fin-diff-grads"),output] =...
    fmincon(fun, x0, [],[],[],[],[],[],[],opts);
timetable.Time("fin-diff-grads") = toc(overallTime);
timetable.Iters("fin-diff-grads") = output.iterations;
timetable.TimePerIter("fin-diff-grads") =...
    timetable.Time("fin-diff-grads")/timetable.Iters("fin-diff-grads");

最小二乗法、ヤコビアンなし

lsqopts = optimoptions("lsqnonlin","Display","none",...
    "MaxFunctionEvaluations",Inf,"MaxIterations",Inf);
fun = @hlsqnonlin0obj;
overallTime = tic;
[~,timetable.Fval("LSQ_NoJacob"),~,timetable.Eflag("LSQ_NoJacob"),output] =...
    lsqnonlin(fun, x0, [],[],lsqopts);
timetable.Time("LSQ_NoJacob") = toc(overallTime);
timetable.Iters("LSQ_NoJacob") = output.iterations;
timetable.TimePerIter("LSQ_NoJacob") =...
    timetable.Time("LSQ_NoJacob")/timetable.Iters("LSQ_NoJacob");

ヤコビアンを使用した最小二乗法

lsqopts.SpecifyObjectiveGradient = true;
overallTime = tic;
[~,timetable.Fval("LSQ_Jacob"),~,timetable.Eflag("LSQ_Jacob"),output] =...
    lsqnonlin(fun, x0, [],[],lsqopts);
timetable.Time("LSQ_Jacob") = toc(overallTime);
timetable.Iters("LSQ_Jacob") = output.iterations;
timetable.TimePerIter("LSQ_Jacob") =...
    timetable.Time("LSQ_Jacob")/timetable.Iters("LSQ_Jacob");

検証結果

disp(timetable)

                       Time        Fval       Eflag    Iters    TimePerIter
                      ______    __________    _____    _____    ___________

    BFGS_NoGrad       110.44    5.0083e-08      1      7137       0.015475 
    LBFGS_NoGrad      53.143     2.476e-07      1      4902       0.010841 
    BFGS_Grad         35.491    2.9865e-08      1      7105      0.0049952 
    LBFGS_Grad        1.2056    9.7505e-08      1      4907      0.0002457 
    Analytic          7.0991     1.671e-10      3      2301      0.0030852 
    fin-diff-grads     5.217    1.1422e-15      1      1382       0.003775 
    LSQ_NoJacob       94.708    3.7969e-25      1      1664       0.056916 
    LSQ_Jacob         6.5225    3.0056e-25      1      1664      0.0039197

タイミング結果から、次のことが示されます。

この問題では、勾配のある LBFGS ヘッセ近似が飛び抜けて最速です。
次に高速なソルバー実行は、勾配のある有限差分ヘッシアンを使用する fmincon、解析的かつ勾配のあるヘッシアンを使用する信頼領域 fminunc、および解析的ヤコビアンを使用する lsqnonlin です。
勾配のない fminunc BFGS アルゴリズムは、ヤコビアンなしの lsqnonlin ソルバーと同じくらいの速度です。この問題の場合、lsqnonlin に必要な反復回数が fminunc より大幅に少なくなっていますが、各反復にかかる時間はかなり長くなっています。
どのソルバーも、微分によって速度に大きな違いが出ています。

補助関数

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

function [f,G,H] = hfminunc0obj(x)
% Rosenbrock function in N dimensions
N = numel(x);
xx = x(1:N-1);
xx_plus = x(2:N);
f_vec = 100*(xx.^2 - xx_plus).^2 + (xx - 1).^2;
f = sum(f_vec);
if (nargout >= 2) % Gradient
    G = zeros(N,1);
    for k = 1:N
        if (k == 1)
            G(k) = 2*(x(k)-1) + 400*x(k)*(x(k)^2-x(k+1));
        elseif (k == N)
            G(k) = 200*x(k) - 200*x(k-1)^2;
        else
            G(k) = 202*x(k) - 200*x(k-1)^2 - 400*x(k)*(x(k+1) - x(k)^2) - 2;
        end
    end
    if nargout >= 3 % Hessian
        H = spalloc(N,N,3*N);
        for i = 2:(N-1)
            H(i,i) = 202 + 1200*x(i)^2 - 400*x(i+1);
            H(i,i-1) = -400*x(i-1);
            H(i,i+1) = -400*x(i);
        end
        H(1,1) = 2 + 1200*x(1)^2 - 400*x(2);
        H(1,2) = -400*x(1);
        H(N,N) = 200;
        H(N,N-1) = -400*x(N-1);
    end
end
end

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

function [f,G] = hlsqnonlin0obj(x)
% Rosenbrock function in N dimensions
N = numel(x);
xx = x(1:N-1);
xx_plus = x(2:N);
f_vec = [10*(xx.^2 - xx_plus), (xx - 1)];
f = reshape(f_vec',[],1); % Vector of length 2*(N-1)
% Jacobian
if (nargout >= 2)
    G = spalloc(2*(N-1),N,3*N); % Jacobian size 2*(N-1)-by-N with 3*N nonzeros
    for k = 1:(N-1)
        G(2*k-1,k) = 10*2*x(k);
        G(2*k-1,k+1) = -10;
        G(2*k,k) = 1;
    end
end
end