効率的な最適化問題の作成

線形問題に整数制約がある場合、solve は intlinprog を呼び出して、解を得ます。より迅速に解を得たり、より多くの整数実行可能点を得たりするためのヒントについては、整数線形計画法の調整を参照してください。

問題を解き始める前に、問題の制約または目的の定式を改善できる場合があります。一般的に、ソフトウェアでは、目的関数または制約用の式をループで作成するより、ベクトル化された形式で作成するほうが高速です。この速度の差は、特に、最適化式が自動微分の影響を受ける場合に大きくなります。Optimization Toolbox の自動微分を参照してください。

通常は、静的解析のための for ループの作成で説明されているように、別個の関数として書き込むことで効率的なループを作成する方が簡単です。この場合は、最適化式の静的解析で説明されているように、fcn2optimexpr を使用して関数を含む最適化式を作成します。

次の目的関数について考えます。

$\sum_{i = 1}^{30} \sum_{j = 1}^{30} \sum_{k = 1}^{10} x_{i, j, k} b_{k} c_{i, j},$

ここで、x、b、c は最適化変数です。この目的関数を定式化する一般的な 2 つの方法は、次のとおりです。

for ループを使用する。

x = optimvar('x',30,30,10);
b = optimvar('b',10);
c = optimvar('c',30,30);
tic
expr = optimexpr;
for i = 1:30
    for j = 1:30
        for k = 1:10
            expr = expr + x(i,j,k)*b(k)*c(i,j);
        end
    end
end
toc

Elapsed time is 307.459465 seconds.

ここで、expr には、目的関数の式が含まれます。ソフトウェアにこの式を効率的に評価させるには、別個の関数としてこれを書き込みます。

function expr = loopme(x,b,c)
expr = 0;
ni = size(x,1);
nj = size(x,2);
nk = size(x,3);
for i = 1:ni
    for j = 1:nj
        for k = 1:nk
            expr = expr + x(i,j,k)*b(k)*c(i,j);
        end
    end
end
end

fcn2optimexpr を使用して式を含めます。

tic
expr = fcn2optimexpr(@loopme,x,b,c,OutputSize=[1,1]);
toc

Elapsed time is 23.888763 seconds.

目的関数を問題に含めます。

problem = optimproblem("Objective",expr);

ベクトル化されたステートメントを使用する。ベクトル化されたステートメントは、一般に静的解析のために変更されていない for ループより迅速に実行されます。ベクトル化されたステートメントは、いくつかの方法で作成できます。
- b と c を展開する。項単位の乗算を可能にするため、x と同じサイズの定数を作成します。
  tic bigb = reshape(b,1,1,10); bigb = repmat(bigb,30,30,1); bigc = repmat(c,1,1,10); expr = sum(sum(sum(x.*bigb.*bigc))); toc
  Elapsed time is 0.013631 seconds.
- b を 1 回ループする。
  tic expr = optimexpr; for k = 1:10 expr = expr + sum(sum(x(:,:,k).*c))*b(k); end toc
  Elapsed time is 0.044985 seconds.
- b をループしてから、ループ後に項の和を求めることによって式を作成する。
  tic expr = optimexpr(30,30,10); for k = 1:10 expr(:,:,k) = x(:,:,k).*c*b(k); end expr = sum(expr(:)); toc
  Elapsed time is 0.039518 seconds.

例最適化変数を使用した静電学における制約付き非線形最適化のベクトル化された実装とベクトル化されていない実装の速度の違いを観察する。この例は、R2020b の自動微分を使用して計測されました。

N = 30;
x = optimvar('x',N,'LowerBound',-1,'UpperBound',1);
y = optimvar('y',N,'LowerBound',-1,'UpperBound',1);
z = optimvar('z',N,'LowerBound',-2,'UpperBound',0);
elecprob = optimproblem;
elecprob.Constraints.spherec = (x.^2 + y.^2 + (z+1).^2) <= 1;
elecprob.Constraints.plane1 = z <= -x-y;
elecprob.Constraints.plane2 = z <= -x+y;
elecprob.Constraints.plane3 = z <= x-y;
elecprob.Constraints.plane4 = z <= x+y;

rng default % For reproducibility
x0 = randn(N,3);
for ii=1:N
    x0(ii,:) = x0(ii,:)/norm(x0(ii,:))/2;
    x0(ii,3) = x0(ii,3) - 1;
end
init.x = x0(:,1);
init.y = x0(:,2);
init.z = x0(:,3);
opts = optimoptions('fmincon','Display','off');

tic
energy = optimexpr(1);
for ii = 1:(N-1)
    jj = (ii+1):N; % Vectorized
    tempe = (x(ii) - x(jj)).^2 + (y(ii) - y(jj)).^2 + (z(ii) - z(jj)).^2;
    energy = energy + sum(tempe.^(-1/2));
end
elecprob.Objective = energy;
disp('Vectorized computation time:')
[sol,fval,exitflag,output] = solve(elecprob,init,'Options',opts);
toc

Vectorized computation time:
Elapsed time is 1.838136 seconds.

tic
energy2 = optimexpr(1); % For nonvectorized comparison
for ii = 1:(N-1)
    for jjj = (ii+1):N; % Not vectorized
        energy2 = energy2 + ((x(ii) - x(jjj))^2 + (y(ii) - y(jjj))^2 + (z(ii) - z(jjj))^2)^(-1/2);
    end
end
elecprob.Objective = energy2;
disp('Non-vectorized computation time:')
[sol,fval,exitflag,output] = solve(elecprob,init,'Options',opts);
toc

Non-vectorized computation time:
Elapsed time is 204.615210 seconds.

ベクトル化されたバージョンは、ベクトル化されていないバージョンより約 100 倍高速です。

静的解析バージョンと比較します。関数 myenergy のコードは、この例の終わりに掲載しています。

elecprob.Objective = fcn2optimexpr(@myenergy,N,x,y,z);

solve を呼び出して問題を解きます。求解時間を測定します。

disp('Static analysis computation time:')
tic
[sol,fval,exitflag,output] = solve(elecprob,init,'Options',opts);
toc

Static analysis computation time:
Elapsed time is 1.293513 seconds.

静的解析は、計算時間を最も抑えることができます。

このコードは関数 myenergy を作成します。

function energy = myenergy(N, x, y, z)
energy = 0;
for ii = 1:(N-1)
    for jj = (ii+1):N
        energy = energy + ((x(ii) - x(jj))^2 + (y(ii) - y(jj))^2 + (z(ii) - z(jj))^2)^(-1/2);
    end
end
end

効率的な最適化問題の作成

関連するトピック