ベクトル化を交えた組み合わせ計算の方法

現在、多変数関数について、各変数を1次元の配列として全ての組み合わせを計算しようとしています。

下記にコードを記しましたが、変数が多くなるにつれて計算に膨大な時間がかかってしまうため高速化したいです。

並列処理はfor文を複数使う際には適用できなかったため、ベクトル化をして、for文を減らせないか検討中です。

下記コードでA,B,Cの各配列から3×3×3のD,Eを変数とした3次元を作成し、計算することは可能でしょうか？

syms A B C D E;

A=1:1:5;

B=1:1:5;

C=1:1:5;

D=1:1:5;

E=1:1:5;

Z=A*B*C*D*E;

OUTPUT=zeros(5,5,5,5,5);

for i=1:3

% for j=1:3

% for k=1:3

% for l=1:3

% for m=1:3

% OUTPUT(i,j,k,l,m)=subs(Z,[A B C D E], [A(1,i) B(1,j) C(1,k) D(1,l) E(1,m)]);

% end

% end

% end

% end

% end

Answer 1

[Akira Agata]

ndgrid を使う方法はいかがでしょうか？

たとえば簡単のため A~E をすべて 1:3 とすると、Z は以下のように計算できます。

% A~E すべての組み合わせに対して　A*B*...*E を計算
[A, B, C, D, E] = ndgrid(1:3);
Z = A.*B.*C.*D.*E;

% 結果を確認
[A(:), B(:), C(:), D(:), E(:), Z(:)]

ans = 243×6

1 1 1 1 1 1 2 1 1 1 1 2 3 1 1 1 1 3 1 2 1 1 1 2 2 2 1 1 1 4 3 2 1 1 1 6 1 3 1 1 1 3 2 3 1 1 1 6 3 3 1 1 1 9 1 1 2 1 1 2

[Shota Ino]

御回答ありがとうございます。

ndgridという関数を初めて見ました。

一点質問なのですが、実際のコードですと、Zを行列計算から算出しており複雑な式となっています。（乗除算時は「.」を使用していません）

上記のようなZに対してもndgridは適用できるのでしょうか？

[Akira Agata]

@Shota Ino様

ご質問の件ですが、おそらく可能と考えております。

具体的にはどのような行列計算になりますでしょうか？

[Shota Ino]

電気回路のF行列の行列計算を実施しており、全ては開示できませんがおおよそ下記のようなコードとなっています。

syms w L C

%実際はF1~F10まであります。以下は例

F1=[1 1i*w*L ; 0 1];

F2=[1 0; 1i*w*C 1];

F_all=F1*F2;

subs(F_all,[w L C],[w1 L1 C1])

[Akira Agata]

ご説明ありがとうございます。

F行列の計算ということは 2次元配列の計算になるかと思いますが、@Shota Ino 様が実施したい「全ての組み合わせの計算」とはどのようなイメージでしょうか。

たとえば上記の F_all = F1*F2 を例に取ると、F1, F2 の各要素（合計 8 個）について、それぞれ取り得る値がN通りあり、すべての組み合わせ（N^8 パターン）の F_all を計算したい、ということでしょうか。

[Shota Ino]

F行列の出力結果は2次元です。

全ての組み合わせとは、F行列内の変数（上記であれば w, L, C)について取りうる値の組み合わせを全て求めたいと思っています。例えば、w, L, Cがそれぞれ5通りの値を取りうるなら、2次元のF行列の結果が5×5×5=125通り出力されるようなイメージです。

[Akira Agata]

ご説明ありがとうございます。

簡単のため、w, L, C がそれぞれ 3 通り (1~3) 、合計 3^3 = 27通り出力されるようなケースを想定すると、以下のようになるかと思います (2通りの手法で出力してみました)。

% w, L, C の全パターンを作成
[w, L, C] = ndgrid(1:3);

% 方法1: pagemtimes を使って一気に計算
n = numel(w);
F1 = repmat([1 nan; 0 1], 1, 1, n);
F2 = repmat([1 0; nan 1], 1, 1, n);

F1(1, 2, :) = 1i*w(:).*L(:);
F2(2, 1, :) = 1i*w(:).*C(:);

Fall = pagemtimes(F1, F2); % -> Fall(:,:,k) が k番目の出力F行列

% 方法2: 後で確認しやすいテーブル型変数に整理しながら作成
tData = table(w(:), L(:), C(:), ...
  'VariableNames', {'w', 'L', 'C'});

tData.F1 = arrayfun(@(x) [1 x; 0 1], 1i*tData.w.*tData.L, ...
  'UniformOutput', false);
tData.F2 = arrayfun(@(x) [1 0; x 1], 1i*tData.w.*tData.C, ...
  'UniformOutput', false);
tData.Fall = cellfun(@(x,y) x*y, tData.F1, tData.F2, ...
  'UniformOutput', false);

disp(tData) % -> tData.Fall{k} が k番目の出力F行列

    w    L    C         F1              F2             Fall    
    _    _    _    ____________    ____________    ____________

    1    1    1    {2×2 double}    {2×2 double}    {2×2 double}
    2    1    1    {2×2 double}    {2×2 double}    {2×2 double}
    3    1    1    {2×2 double}    {2×2 double}    {2×2 double}
    1    2    1    {2×2 double}    {2×2 double}    {2×2 double}
    2    2    1    {2×2 double}    {2×2 double}    {2×2 double}
    3    2    1    {2×2 double}    {2×2 double}    {2×2 double}
    1    3    1    {2×2 double}    {2×2 double}    {2×2 double}
    2    3    1    {2×2 double}    {2×2 double}    {2×2 double}
    3    3    1    {2×2 double}    {2×2 double}    {2×2 double}
    1    1    2    {2×2 double}    {2×2 double}    {2×2 double}
    2    1    2    {2×2 double}    {2×2 double}    {2×2 double}
    3    1    2    {2×2 double}    {2×2 double}    {2×2 double}
    1    2    2    {2×2 double}    {2×2 double}    {2×2 double}
    2    2    2    {2×2 double}    {2×2 double}    {2×2 double}
    3    2    2    {2×2 double}    {2×2 double}    {2×2 double}
    1    3    2    {2×2 double}    {2×2 double}    {2×2 double}
    2    3    2    {2×2 double}    {2×2 double}    {2×2 double}
    3    3    2    {2×2 double}    {2×2 double}    {2×2 double}
    1    1    3    {2×2 double}    {2×2 double}    {2×2 double}
    2    1    3    {2×2 double}    {2×2 double}    {2×2 double}
    3    1    3    {2×2 double}    {2×2 double}    {2×2 double}
    1    2    3    {2×2 double}    {2×2 double}    {2×2 double}
    2    2    3    {2×2 double}    {2×2 double}    {2×2 double}
    3    2    3    {2×2 double}    {2×2 double}    {2×2 double}
    1    3    3    {2×2 double}    {2×2 double}    {2×2 double}
    2    3    3    {2×2 double}    {2×2 double}    {2×2 double}
    3    3    3    {2×2 double}    {2×2 double}    {2×2 double}

Answer 2

[sp6038sy]

変数の組み合わせを事前に計算し行列演算するのはどうでしょうか。

例) 多変数関数 の場合

syms A B Z;

% 多変数関数
Z = A.*B;

% 変数の組み合わせを計算
A = 1:3;
B = 1:3;
[A, B] = ndgrid(A, B)

A = 3×3

1 1 1 2 2 2 3 3 3

B = 3×3

1 2 3 1 2 3 1 2 3

% 多変数関数に代入
out = subs(Z)

out = 

[Shota Ino]

回答ありがとうございます。

ZをA, Bの組み合わせからなる2次元配列にすることができました。ありがとうございました。しかし、C,D,Eの変数についてfor文で組み合わせ計算を実施する際にサイズ違いによるエラーが出てしまいました。

現在解決手段を模索中です