Main Content

非負の線形最小二乗法、問題ベース

ライブ スクリプトを開く

この例では、いくつかのアルゴリズムを使用して、解が非負となる範囲制約のある、線形最小二乗問題を解く方法を説明します。線形最小二乗問題は次の形式をとります。

minx‖Cx-d‖2.

この場合、解が非負 x≥0 となるように制約します。

はじめに、配列 C および d をワークスペースに読み込みます。

load particle

各配列のサイズを表示します。

sizec = size(C)
sizec = 1×2

        2000         400


sized = size(d)
sized = 1×2

        2000           1

C による乗算用に適切なサイズの最適化変数 x を作成します。x の要素に対して下限 0 を設定します。

x = optimvar('x',sizec(2),'LowerBound',0);

目的関数の式を作成します。

residual = C*x - d;
obj = sum(residual.^2);

nonneglsq という名前の最適化問題を作成し、目的関数を問題に含めます。

nonneglsq = optimproblem('Objective',obj);

問題に対する既定のソルバーを求めます。

opts = optimoptions(nonneglsq)
opts = 
  lsqlin options:

   Options used by current Algorithm ('interior-point'):
   (Other available algorithms: 'active-set', 'trust-region-reflective')

   Set properties:
     No options set.

   Default properties:
              Algorithm: 'interior-point'
    ConstraintTolerance: 1.0000e-08
                Display: 'final'
           LinearSolver: 'auto'
          MaxIterations: 200
    OptimalityTolerance: 1.0000e-08
          StepTolerance: 1.0000e-12

   Options not used by current Algorithm ('interior-point')
   Default properties:
      FunctionTolerance: 2.2204e-14
    JacobianMultiplyFcn: []
         ObjectiveLimit: -1.0000e+20
    SubproblemAlgorithm: 'cg'
               TypicalX: 'ones(numberOfVariables,1)'

既定のソルバーを使用して、問題を解きます。

[sol,fval,exitflag,output] = solve(nonneglsq);
Solving problem using lsqlin.

Minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

最適化プロセスの詳細を確認するには、出力構造体を検証します。

disp(output)
            message: 'Minimum found that satisfies the constraints....'
          algorithm: 'interior-point'
      firstorderopt: 9.9673e-07
    constrviolation: 0
         iterations: 9
       linearsolver: 'sparse'
       cgiterations: []
             solver: 'lsqlin'

この出力構造体は、lsqlin ソルバーが内点法アルゴリズムの内部的な線形ソルバーとしてスパースを使用し、反復を 9 回実行して、1 次の最適性の尺度が約 1e-6 に到達したことを示しています。

アルゴリズムの変更

信頼領域 Reflective 法アルゴリズムは、範囲制約付き問題を扱います。次の問題に対して、うまく機能することを確認します。

opts.Algorithm = 'trust-region-reflective';
[sol2,fval2,exitflag2,output2] = solve(nonneglsq,'Options',opts);
Solving problem using lsqlin.

Local minimum possible.

lsqlin stopped because the relative change in function value is less than the function tolerance.

disp(output2)
          algorithm: 'trust-region-reflective'
         iterations: 14
      firstorderopt: 5.2620e-08
       cgiterations: 54
    constrviolation: []
       linearsolver: []
            message: 'Local minimum possible....'
             solver: 'lsqlin'

このとき、ソルバーは反復回数を増やし、1 次の最適性の尺度がより低い (良い) 解に到達します。

さらに良好な 1 次の最適性の尺度を得るには、SubproblemAlgorithm オプションを 'factorization' に設定します。

opts.SubproblemAlgorithm = 'factorization';
[sol3,fval3,exitflag3,output3] = solve(nonneglsq,'Options',opts);
Solving problem using lsqlin.

Optimal solution found.

disp(output3)
          algorithm: 'trust-region-reflective'
         iterations: 11
      firstorderopt: 1.3973e-14
       cgiterations: 0
    constrviolation: []
       linearsolver: []
            message: 'Optimal solution found.'
             solver: 'lsqlin'

このオプションを使用することにより、1 次の最適性の尺度は、最良であるほぼゼロになります。

ソルバーの変更

lsqnonneg ソルバーは、非負の線形最小二乗法を処理するために特別に設計されています。このソルバーを試してみましょう。

[sol4,fval4,exitflag4,output4] = solve(nonneglsq,'Solver','lsqnonneg');
Solving problem using lsqnonneg.

disp(output4)
    iterations: 184
     algorithm: 'active-set'
       message: 'Optimization terminated.'
        solver: 'lsqnonneg'

lsqnonneg は 1 次の最適性の尺度を報告しません。代わりに、fval 出力で返される残差ノルムを調査します。下位の有効数字を確認するには、各残差ノルムから 22.5794 を減算します。

t = table(fval - 22.5794, fval2 - 22.5794, fval3 - 22.5794, fval4 - 22.5794,...
    'VariableNames',{'default','trust-region-reflective','factorization','lsqnonneg'})
t=1×4 table
     default      trust-region-reflective    factorization    lsqnonneg 
    __________    _______________________    _____________    __________

    5.0804e-05          4.9179e-05            4.9179e-05      4.9179e-05

既定のソルバーは、残差ノルムが他の 3 つよりも高く (悪く) なっています。他の 3 つの残差ノルムは、このレベルの表示精度では区別できないほどの違いしかありません。最も低い残差ノルムを確認するには、2 つの結果から lsqnonneg の結果を減算します。

disp([fval2 - fval4,fval3 - fval4])
   1.0e-12 *

    0.7176    0.7034

量の違いはほとんど無視できるレベルであるため、lsqnonneg 残差ノルムは最小値となります。しかし、lsqnonneg は、収束までの反復回数が最大です。

参考

関連するトピック