中身が2行だけ書かれたMATLABプログラムを全力解説

※こちらは，会社の技術ブログとのクロスポスト記事です．元の記事 ➡ 前編，中編，後編

中身が2行だけ書かれた，以下のMATLABプログラムについて解説します．　

function W = SLD(H,M)
  W = nchoosek(1:H+M-1,M-1) - repmat(0:M-2,nchoosek(H+M-1,M-1),1) - 1;
  W = ([W,zeros(size(W,1),1)+H]-[zeros(size(W,1),1),W])/H;
end

こちらは，とあるプログラムから2行だけ抜き出して作成した関数です．
プログラムを動かすには，$H$: 1以上の整数，$M$: 2以上の整数 を指定する必要があります．

以上の情報だけで，どういうプログラムか読み解けた方は，かなり凄いです．
それでは，全力解説していきます．

MATLABの前提と入出力

MATLABの関数はfunctionで始まり，endで終わるのが作法です．
関数の先頭は左からfunction，[出力]，=，関数名，(入力)です．

適当に値を入力すると，以下の出力が得られます．

>> SLD(1,2)
ans =
     0     1
     1     0

>> SLD(1,3)
ans =
     0     0     1
     0     1     0
     1     0     0

>> SLD(1,4)
ans =
     0     0     0     1
     0     0     1     0
     0     1     0     0
     1     0     0     0

>> SLD(2,2)
ans =
         0    1.0000
    0.5000    0.5000
    1.0000         0

>> SLD(3,2)
ans =
         0    1.0000
    0.3333    0.6667
    0.6667    0.3333
    1.0000         0

$M=2$，$H=\{9,14,20\}$のとき，出力$W$の散布図は以下になります．

$M=3$，$H=\{3,4,5\}$のとき，出力$W$の散布図は以下になります．

点の数は，左から順に10, 15, 21です．処理内容が見えてきましたね．
このプログラムは，$M$次元空間の単位超平面上に均一分布する点群の生成プログラムです．
$M=\{2,3\}$次元の場合，

• $M=2$: 2つの点$(0,1)$と$(1,0)$を結ぶ直線上に均一分布する点群が得られます．
• $M=3$: 3つの点$(0,0,1),(0,1,0),(1,0,0)$の3点を頂点とする平面上に均一分布する点群が得られます．

Simplex-lattice design

このプログラムが生成する点群をSimplex-lattice design (SLD)といいます．
点群SLDは，「進化計算による多目的最適化」という研究分野でよく利用されています．
今回解説する2行は，PlatEMOという，同研究分野でよく利用されるライブラリから抜き出したものです．抜き出したプログラムはこちらです．

点の各要素は$\{0/H, 1/H,\dots, H/H\}$のいずれかの値を取る，点の要素の総和は1になる，点群サイズは$N=C^{M-1}_{H+M-1}$で計算できるという特徴があります．
$M$，$H$と点群サイズ$N$の関係性を表にすると以下になります．

$M$ \ $H$	1	2	3	4	5	6	7	8	9
2	2	3	4	5	6	7	8	9	10
3	3	6	10	15	21	28	36	45	55
4	4	10	20	35	56	84	120	165	220
5	5	15	35	70	126	210	330	495	715
6	6	21	56	126	252	462	792	1287	2002
7	7	28	84	210	462	924	1716	3003	5005
8	8	36	120	330	792	1716	3432	6435	11440

以上で，プログラムの入出力とその性質，利用される分野の解説は終了です．
ここから，プログラムを読み解くために必要なMATLABの知識を解説します．

必要なMATLAB知識

前提

MATLABは，多次元配列を扱えるように設計されています．
そのため，説明に「配列」と書かれることが多いです．
この記事では，分かり易さのため，任意の次元の配列を指す場合のみ「配列」と書きます．
各用語の対応は，以下の通りです．

• 1次元配列，1階のテンソル　→　ベクトル
• 2次元配列，2階のテンソル　→　行列
• 特に，行数が1の2次元配列　→　行ベクトル
• 特に，列数が1の2次元配列　→　列ベクトル
• 任意の次元の配列　　　　　→　配列

MATLABは，1つの関数が複数の機能を持っていることがあります．
nchoosekは，2つの機能を持っており，どちらも解説します．
それ以外の関数は，プログラムを読み解くために必要な機能のみ解説します．

解説の理解に詰まった場合，まずは演算子と特殊文字や配列と行列の演算を見てみることをお勧めします．算術演算子(+-*/)の説明は，省略します．

明示されていない組み込み関数

MATLABは，変数を同じように操作する方法が複数提供されている場合があります．
単純操作の多くに代替法が設定されており，関数として意識することなく利用できます．
ここでは，明示されていない基本的な組み込み関数を4つ解説します．

• colon, :：2変数を始点と終点とする，公差1の等差数列を返します．行ベクトルとして返します．
  • 右の3つは，全て同じ変数を表します．[1 2 3 4]　1:4　colon(1,4)
• vertcat, ;：垂直に連結した配列を返します．
  • 右の2つは，同じ変数を表します．[A; B]　vertcat(A,B)
• horzcat：水平に連結した配列を返します．
  • 右の3つは，全て同じ変数を表します．[A B]　[A,B]　horzcat(A,B)
• ctranspose, '： 要素に複素数が無ければ，ベクトルまたは行列の転置を返します．（transpose, .'と同じ）
  • 右の5つは，全て同じ変数を表します．[1; 2; 3]　[1 2 3]'　[1 2 3].'　transpose([1 2 3])　ctranspose([1 2 3])

: はベクトルの作成の他，インデックス，for ループの反復でも使われます．
; は配列の連結の他，行末の指定，コード行の出力の非表示でも使われます．

明示されている組み込み関数

明示されているMATLABの組み込み関数は，プログラム内で色付き表示されます．(Qiitaでは未対応？)
ここでは，プログラムに登場する4つの組み込み関数を解説します．

• nchoosek：二項係数 or すべての組み合わせを返します．
  • 右の2つは，同じ変数を表します．6　nchoosek(4,2)
  • 右の2つは，同じ変数を表します．[[1 2];[1 3];[1 4];[2 3];[2 4];[3 4]]　nchoosek(1:4,2)
• repmat：配列のコピーを返します．第三引数が1の場合，第一引数をコピー元として，行数が第二引数倍されたコピーを返します．
  • 右の2つは，同じ変数を表します．[[9 8];[9 8];[9 8]]　repmat([9 8],3,1)
• size：配列サイズを返します．第二引数が1の場合，配列の行数を返します．
  • 右の2つは，同じ変数を表します．3　size([9; 8; 7],1)
• zeros：すべての要素が0の配列を返します．左の引数で行数，右の引数で列数を指定します．
  • 右の2つは，同じ変数を表します．[0; 0; 0]　zeros(3,1)

演算子の優先順位

演算子の優先順位は，ここを参照ください．
今回のプログラムでは，以下の3つの基準だけ覚えれば大丈夫です．

• 小かっこ ( ) の中の演算が最優先．
• コロン演算子 (:)より加算 (+)，減算 (-)を優先する．
• 演算子は，左から右へ評価する．

以上の知識があれば，プログラムを読み解くことができるはずです．

さて，準備が整いました．いよいよプログラムの解説に入ります．

プログラム解説

$H=4$，$M=3$を例にプログラムを解説します．
処理を1個ずつ進めながら，プログラムの変化を見ていきます．
1行目，2行目で分けて考えます．

1行目の処理

処理開始前．初期状態．

  W = nchoosek(1:H+M-1,M-1) - repmat(0:M-2,nchoosek(H+M-1,M-1),1) - 1;

最初に，プログラムを以下のように変形します．

  A = nchoosek(1:H+M-1,M-1);
  B = repmat(0:M-2,nchoosek(H+M-1,M-1),1);
  W = A - B - 1;

次に，括弧の中でカンマで区切られた要素の四則演算をします．

  A = nchoosek(1:6,2);
  B = repmat(0:1,nchoosek(6,2),1);
  W = A - B - 1;

次に，括弧の中でカンマで区切られた要素について，コロン演算子(:)を適用し行ベクトルを作成します．

  A = nchoosek([1 2 3 4 5 6],2);
  B = repmat([0 1],nchoosek(6,2),1);
  W = A - B - 1;

次に，nchoosekを利用して，行ベクトル[1 2 3 4 5 6]から2要素取り出す組み合わせを計算します．

  A = [[1 2];[1 3];[1 4];[1 5];[1 6];[2 3];[2 4];[2 5];[2 6];[3 4];[3 5];[3 6];[4 5];[4 6];[5 6]];
  B = repmat([0 1],nchoosek(6,2),1);
  W = A - B - 1;

結果，$A$は15行2列の行列になります．散布図としてプロットすると以下が得られます．

次に，nchoosekを利用して，6個の対象から2個選ぶ選び方の総数，二項係数を計算します．

  A = [[1 2];[1 3];[1 4];[1 5];[1 6];[2 3];[2 4];[2 5];[2 6];[3 4];[3 5];[3 6];[4 5];[4 6];[5 6]];
  B = repmat([0 1],15,1);
  W = A - B - 1;

次に，repmatを利用して，行ベクトル[0 1]のコピーを作成します．$B$も15行2列の行列になります．

  A = [[1 2];[1 3];[1 4];[1 5];[1 6];[2 3];[2 4];[2 5];[2 6];[3 4];[3 5];[3 6];[4 5];[4 6];[5 6]];
  B = [[0 1];[0 1];[0 1];[0 1];[0 1];[0 1];[0 1];[0 1];[0 1];[0 1];[0 1];[0 1];[0 1];[0 1];[0 1]];
  W = A - B - 1;

次に，同じ要素数の行列$A$，$B$の各要素について，引き算をします．

  W = [[1 1];[1 2];[1 3];[1 4];[1 5];[2 2];[2 3];[2 4];[2 5];[3 3];[3 4];[3 5];[4 4];[4 5];[5 5]] - 1;

そして，行列の各要素を-1します．

  W = [[0 0];[0 1];[0 2];[0 3];[0 4];[1 1];[1 2];[1 3];[1 4];[2 2];[2 3];[2 4];[3 3];[3 4];[4 4]];

以上で，1行目の処理が全て終了し，15行2列の行列$W$が作成されます．
散布図としてプロットすると以下が得られます．

2行目の処理

処理開始前．初期状態．

  W = [[0 0];[0 1];[0 2];[0 3];[0 4];[1 1];[1 2];[1 3];[1 4];[2 2];[2 3];[2 4];[3 3];[3 4];[4 4]];
  W = ([W,zeros(size(W,1),1)+H]-[zeros(size(W,1),1),W])/H;

最初に，プログラムを以下のように変形します．

  W = [[0 0];[0 1];[0 2];[0 3];[0 4];[1 1];[1 2];[1 3];[1 4];[2 2];[2 3];[2 4];[3 3];[3 4];[4 4]];
  C = zeros(size(W,1),1);
  D = C+H;
  E = [W,D];
  F = [C,W];
  W = (E-F)/H;

次に，sizeを利用して，行列$W$の行数を計算します．

  W = [[0 0];[0 1];[0 2];[0 3];[0 4];[1 1];[1 2];[1 3];[1 4];[2 2];[2 3];[2 4];[3 3];[3 4];[4 4]];
  C = zeros(15,1);
  D = C+H;
  E = [W,D];
  F = [C,W];
  W = (E-F)/H;

次に，$C$を計算します．（zerosを利用して，すべての要素が0の列ベクトルを作成します．）

  W = [[0 0];[0 1];[0 2];[0 3];[0 4];[1 1];[1 2];[1 3];[1 4];[2 2];[2 3];[2 4];[3 3];[3 4];[4 4]];
  C = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]';
  D = C+H;
  E = [W,D];
  F = [C,W];
  W = (E-F)/H;

次に，$D$を計算します．（列ベクトル$C$の各要素を$+H(=4)$します．）

  W = [[0 0];[0 1];[0 2];[0 3];[0 4];[1 1];[1 2];[1 3];[1 4];[2 2];[2 3];[2 4];[3 3];[3 4];[4 4]];
  C = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]';
  D = [4 4 4 4 4 4 4 4 4 4 4 4 4 4 4]';
  E = [W,D];
  F = [C,W];
  W = (E-F)/H;

次に，$E$と$F$を計算します．（同じ行数である$W$と$D$，$C$と$W$を水平連結します．）

  E = [[0 0 4];[0 1 4];[0 2 4];[0 3 4];[0 4 4];[1 1 4];[1 2 4];[1 3 4];[1 4 4];[2 2 4];[2 3 4];[2 4 4];[3 3 4];[3 4 4];[4 4 4]];
  F = [[0 0 0];[0 0 1];[0 0 2];[0 0 3];[0 0 4];[0 1 1];[0 1 2];[0 1 3];[0 1 4];[0 2 2];[0 2 3];[0 2 4];[0 3 3];[0 3 4];[0 4 4]];
  W = (F-G)/H;

最後に，行列を引き算し，各要素を$H$で割って正規化します．
これで，単位超平面に均一分布する点群SLDが得られます．

$$
W\quad=\quad\left(
\begin{matrix}
[0\ 0\ 4]\cr
[0\ 1\ 4]\cr
[0\ 2\ 4]\cr
[0\ 3\ 4]\cr
[0\ 4\ 4]\cr
[1\ 1\ 4]\cr
[1\ 2\ 4]\cr
[1\ 3\ 4]\cr
[1\ 4\ 4]\cr
[2\ 2\ 4]\cr
[2\ 3\ 4]\cr
[2\ 4\ 4]\cr
[3\ 3\ 4]\cr
[3\ 4\ 4]\cr
[4\ 4\ 4]\cr
\end{matrix}\;-\;\;
\begin{matrix}
[0\ 0\ 0]\cr
[0\ 0\ 1]\cr
[0\ 0\ 2]\cr
[0\ 0\ 3]\cr
[0\ 0\ 4]\cr
[0\ 1\ 1]\cr
[0\ 1\ 2]\cr
[0\ 1\ 3]\cr
[0\ 1\ 4]\cr
[0\ 2\ 2]\cr
[0\ 2\ 3]\cr
[0\ 2\ 4]\cr
[0\ 3\ 3]\cr
[0\ 3\ 4]\cr
[0\ 4\ 4]\cr
\end{matrix}
\right)\;/\;4\quad=\quad
\left(
\begin{matrix}
[0\ 0\ 4]\cr
[0\ 1\ 3]\cr
[0\ 2\ 2]\cr
[0\ 3\ 1]\cr
[0\ 4\ 0]\cr
[1\ 0\ 3]\cr
[1\ 1\ 2]\cr
[1\ 2\ 1]\cr
[1\ 3\ 0]\cr
[2\ 0\ 2]\cr
[2\ 1\ 1]\cr
[2\ 2\ 0]\cr
[3\ 0\ 1]\cr
[3\ 1\ 0]\cr
[4\ 0\ 0]\cr
\end{matrix}
\right)\;/\;4\quad=\quad
\begin{matrix}
[\quad\ 0\quad\ \ \ 0\quad\ \ \ 1]\cr
[\quad\ 0\ \ 0.25\ \ 0.75]\cr
[\quad\ 0\quad0.5\quad0.5]\cr
[\quad\ 0\ \ 0.75\ \ 0.25]\cr
[\quad\ 0\quad\ \ \ 1\quad\ \ \ 0]\cr
[0.25\quad\ \ \ 0\ \ 0.75]\cr
[0.25\ \ 0.25\quad0.5]\cr
[0.25\quad0.5\ \ 0.25]\cr
[0.25\ \ 0.75\quad\ \ \ 0]\cr
[\ \ 0.5\quad\ \ \ 0\quad0.5]\cr
[\ \ 0.5\ \ 0.25\ \ 0.25]\cr
[\ \ 0.5\quad0.5\quad\ \ \ 0]\cr
[0.75\quad\ \ \ 0\ \ 0.25]\cr
[0.75\ \ 0.25\quad\ \ \ 0]\cr
[\quad\ 1\quad\ \ \ 0\quad\ \ \ 0]\cr
\end{matrix}
$$

以上で，プログラムは終了です．

終わりに

解説するのはたった2行なのに，内容がすごく膨らんで，自分でも驚きました．
基礎から解説したため，MATLABに触れたことがない人でもわかる記事が書けたと思います．

MATLABは，非常に短く簡潔にプログラムが書ける利点があります．
それが，今回解説した2行のプログラムに良く表れていると思います．
また，MATLABは，日本語で書かれた公式のドキュメント，教材が充実しており，学習も非常にしやすいです．
更に，公式が提供するファイル共有サイト，GitHubに多数のコードが公開されています．
ぜひ一度，MATLABを触ってみてください．

最後に，この内容を更に詳しく知りたいという方は，元のプログラムを眺めてみることをお勧めします．プログラムでは以下の論文がReferenceとして紹介されています．

• Y. Tian, X. Xiang, X. Zhang, R. Cheng, and Y. Jin, “Sampling reference points on the Pareto fronts of benchmark multi-objective optimization problems,” Proceedings of the IEEE Congress on Evolutionary Computation (CEC2018), 2018.
• T. Takagi, K. Takadama, and H. Sato, “Incremental lattice design of weight vector set,” Proceedings of the Genetic and Evolutionary Computation Conference Companion (GECCO2020), pp. 1486-1494, 2020.

利用したコードリポジトリのReferenceは以下の通りです．
Y. Tian, R. Cheng, X. Zhang, and Y. Jin, “PlatEMO: A MATLAB Platform for Evolutionary Multi-Objective Optimization [Educational Forum],” IEEE Computational Intelligence Magazine, Vol. 12, No. 4, pp. 73-87, 2017.