この記事では，Largest Small $n$-Polygons を生成する，以下のプログラムを紹介します．
Largest Small $n$-Polygons の詳細については，以下の記事をご確認ください．
紹介プログラムは，入力$n$に対して「いずれの2点間の距離も1以下で面積が最大となる$n$角形 ⇒ LSP($n$)を出力します．言語はMATLABです．
LSP(n): nが奇数の場合
$n=$奇数の場合，LSP($n$)は正多角形（＝正$n$角形）になります．
プログラムでは，2種類の方法が書かれています．
1つ目は，非線形最適化ソルバーfminconを利用する方法です．
後述する$n=$偶数の場合と同じ方法で，最適なLSP($n$)の近似解を求めます．
動作方法は，X = LSP(n)のように，LSP.mファイルを呼び出すだけです．
2つ目は，直径1の円周上に正$n$角形を作成し，それをスケーリング・移動する方法です．
main.mファイルには，以下のように書かれています．
theta = (1 : n) * 2*pi/n;
X = 0.5 * [sin(theta)' cos(theta)'];
X = X / max(pdist(X)); % Scaling
X(:,2) = X(:,2) - max(X(:,2)) + 0.5; % Moving
まず，角度$\boldsymbol{\theta}=(\theta_1,\theta_2,\dots,\theta_n)$を作成します．角度$\boldsymbol{\theta}$は，以下の通りです．
$$\boldsymbol{\theta}=\left(\frac{1}{n}2\pi,\frac{2}{n}2\pi,\dots,\frac{n}{n}2\pi\right)$$
極座標を直交座標に変更し，半径$r=0.5$の円周上に分布する正$n$角形の座標を獲得します．
このままでは面積が最大とならないので，2点間の距離を関数 pdist で計算し，その最大距離が1になるように図形をスケーリングします．
最後に，他の図形との比較を容易にするため，$y$座標の最大値が$0.5$になるように全体を移動します．
$n=5,7$で得られる図形は，以下の通りです．
 |
 |
|---|---|
| スケーリング前の正五角形 ($r=0.5$) 面積：$0.594410$ |
スケーリング前の正七角形 ($r=0.5$) 面積：$0.684103$ |
 |
 |
| LSP(5)：正五角形 ($r=0.525731$) 面積：$0.657164$ |
LSP(7)：正七角形 ($r=0.512858$) 面積：$0.719741$ |
スケーリングによって正五角形の外接円の半径$r$が大きくなり，正五角形の面積も大きくなっています．
また，正七角形も同様に$r$と面積がも大きくなっていますが，その上がり幅は減少しています．
$n$が大きくなると正多角形は円に近づくので，最終的に$r=0.5$，面積$\frac{1}{4}\pi=0.78539816$に収束します．
ソルバーを利用する方法で得られた結果と正$n$角形を利用する方法で得られた結果は，少数第6位くらいで値に違いが見られます．
厳密な解が欲しいときは，LSP.mファイルを利用せず，main.mファイルに記述のある正多角形変換法を利用した方が良いです．
LSP(n): nが偶数の場合
$n=$偶数の場合，LSP($n$)は正多角形（＝正$n$角形）ではありません．（$n=4$を除く）
$n=4$の場合のみ，LSP($n$)は無限に存在し，その1つが正方形＝正多角形になります．
そのため，奇数の場合のように正多角形の変換，一意の操作でLSP($n$)を作成できません．
LSP($n$)の生成には，一般的に非線形最適化ソルバーを利用します．
この記事で紹介するLSP生成プログラムLSP.mファイルは，ソルバー fminconを利用します．
そのため，Optimization Toolboxが必要です．以下は，LSP.mの全文です．
function X = LSP(n)
% Largest Small Polygon (LSP)
if mod(n, 2)
update_x = @(x) [0.0 0.5; x; flipud([-x(:,1) x(:,2)])];
else
update_x = @(x) [0.0 0.5; x; 0 -0.5; flipud([-x(:,1) x(:,2)])];
end
rng(42);
d = ceil(n/2 - 1);
theta = (1 : d) * 2*pi/n + 0.1*rand(1, d);
x0 = 0.5 * [sin(theta)' cos(theta)'];
lb = -0.5 * ones(d, 2);
ub = 0.5 * ones(d, 2);
options = optimoptions('fmincon', 'Display', 'none', ...
'Algorithm', 'sqp', 'TolFun', 1e-9, 'TolX', 1e-9, ...
'MaxFunEvals', 1e6, 'MaxIter', 1e6);
x = fmincon(@(x) -polygon_area(x, update_x), x0, [], [], [], [], ...
lb, ub, @(x) distance_constraints(x, update_x), options);
X = update_x(x);
end
function A = polygon_area(x, update_x)
X = update_x(x);
A = polyarea(X(:,1), X(:,2));
end
function [c, ceq] = distance_constraints(x, update_x)
X = update_x(x);
c = pdist(X)' - 1;
ceq = [];
end
LSP.mは，3つの関数から構成されており，空行を含めて全部で36行です．
関数LSPは，値$n$を受け取り，$n\times2$の点群$X$を返します．
点群$X$にはLSP($n$)の頂点が格納されており，その1行目は$(0.0, 0.5)$で固定されています．
点群$X$の中身は，時計周りにソートされています．
LSP.mでは，LSPの性質を利用した変数$\boldsymbol{x}$の次元削減トリックを使っています．
まず，1つの頂点を座標$(0.0, 0.5)$に固定することで，変数から除外します．
$n$が偶数の場合，更にもう1つの頂点を$(0.0, -0.5)$に固定することで，変数から除外します．
そして，LSPの線対称性を利用して，探索する変数を更に半分に削減します．
例えばLSP(6)の場合，変数$\boldsymbol{x}$に該当するのは，2点のxy座標だけです．
2点の座標$(x_1, y_1)$，$(x_2, y_2)$が決まれば，残りの点の座標は$(0.0, 0.5)$，$(0.0, -0.5)$，$(-x_1, y_1)$，$(-x_2, y_2)$として一意に定まります．
最終的に，変数$\boldsymbol{x}$の次元$d=\lfloor n/2-1\rfloor$になります．
変数$\boldsymbol{x}$から$n$角形のすべての頂点座標（＝点群$X$）を導出する処理がupdate_xになります．
ソルバーが最適化する目的関数値は，得られた点群$X$から作成した$n$角形の面積Aになります．
最小化問題として解くため，目的関数は-polygon_areaにします．
制約関数のdistance_constraintsですが，等式制約ceqは設定しません．
不等式制約cは，各点の間の距離を設定します．
あらゆる2点間の距離が1以下なら，制約を満足します．
最適化の初期値x0には，正多角形にランダムノイズを加えた点群を利用します．
ソルバー内でも乱数は使われますが，初期値にも乱数を利用します．
乱数シードを固定しないと，毎回微妙に（$10^{-5}$程度）異なる結果が得られます．
変数の下限lbと上限ubには，それぞれ$\boldsymbol{-0.5}$と$\boldsymbol{0.5}$を設定します．
optionsは，最適化アルゴリズムやその終了条件，ログ設定などを記載します．
最後に，関数fminconに上記の変数を入れて，最適化を実行します．
導出された面積最大となる最適な点群$X$をLSP($N$)として出力します．
以上のプログラムでLSPが導出できます．
$n\leq100$までなら，少数第5位まで正確な近似解を導出できるはずです．
$n>100$の場合，実行に時間がかかったりエラーが出たりします．
$n=6,8,10,20,50,100$のとき，その面積は以下になりました．
| 図形 | $n=6$ | $n=8$ | $n=10$ | $n=20$ | $n=50$ | $n=100$ |
|---|---|---|---|---|---|---|
| 正$n$角形 | 0.649519 | 0.707107 | 0.734732 | 0.772542 | 0.783333 | 0.784881 |
| LSP($n$) | 0.674981 | 0.726868 | 0.749137 | 0.776859 | 0.784017 | 0.785048 |
LSP(6): Graham's Biggest Little Hexagon
LSP(6)は，Graham's Biggest Little Hexagonという名前がついています．
この多角形のみ，1変数の最適化によって作図可能です．
LSP($n$)に対する汎用性はありませんが，1変数によって六角形の面積がどのように変化するか，検証します．
まず，$n=6$は偶数なので，2つの頂点は，$(0.0,0.5)$と$(0.0,-0.5)$になります．
次に，点$(0.0, 0.5)$からの距離が1で，点$(0.0, -0.5)$と隣接する頂点の座標を決定します．
ここで，$x$座標を変数$x_1$すると，$y$座標$y_1=0.5-\sqrt{1-x_1^2}$と書けます．
点$(0.0, -0.5)$と隣接するもう一つの頂点は，$(-x_1,y_1)$に決まります．
残りの2頂点は，$x$座標＝-0.5と0.5で，$y$座標$y_2=y_1+\sqrt{1-(0.5+x_1)^2}$になります．
以上をまとめると，最終的に六角形の頂点座標は，以下になります．（時計周り）
| 番号 | 座標$x$ | 座標$y$ |
|---|---|---|
| 1 | $0.0$ | $0.5$ |
| 2 | $0.5$ | $y_2=0.5-\sqrt{1-x_1^2}+\sqrt{1-(0.5+x_1)^2}$ |
| 3 | $x_1$ | $y_1=0.5-\sqrt{1-x_1^2}$ |
| 4 | $0.0$ | $-0.5$ |
| 5 | $-x_1$ | $y_1=0.5-\sqrt{1-x_1^2}$ |
| 6 | $-0.5$ | $y_2=0.5-\sqrt{1-x_1^2}+\sqrt{1-(0.5+x_1)^2}$ |
プログラムなら，$x_1\Rightarrow$xとして，以下のように書けます．
y1 = 0.5 - sqrt(1 - x^2);
y2 = y1 + sqrt(1 - (0.5 + x)^2);
X = [0.0 0.5 x 0.0 -x -0.5; 0.5 y2 y1 -0.5 y1 y2]';
$x_1=0.0,0.1,\dots,0.5$に設定した時の六角形は，以下になります．
 |
 |
|---|---|
| $x_1=0.0$，面積：$0.500000$ | $x_1=0.1$，面積：$0.577995$ |
 |
 |
| $x_1=0.2$，面積：$0.636767$ | $x_1=0.3$，面積：$0.670788$ |
 |
 |
| $x_1=0.4$，面積：$0.666007$ | $x_1=0.5$，面積：$0.500000$ |
$x_1$と六角形の面積$A$の関係は，以下になります．
図の横軸が$x_1$の値，縦軸が面積$A$の値を表します．
この図からも，$x_1=0.343771453...$の時に最大の面積$A=0.674981...$になることが直感的に理解できるかと思います．
上記のGraham's Biggest Little Hexagonに関する検証は，関数 GBLH としてまとめられています．
終わりに
Largest Small $n$-Polygons の生成プログラムを紹介しました．
プログラムの行数に対して，説明を丁寧に書き過ぎたかもしれません．
LSPの説明と記事を分割すると決めてから，分量がかなり増大してしまいました．
紹介したプログラムですが，MATLABが優秀なのか，予想よりも簡潔に書けました．
LSP.mは，幾つかの2点間の距離が$1$になる等式条件を追加したり，1点は必ず$(0.5,y)$の値を取る条件を追加したら，解の精度が向上すると考えられます．
是非，チューニングや新しい工夫を試して，2021年の論文と結果を比較してみてください．
どこかで工夫をしないと，LSP(1000)は解けないと思います．
（最後に紹介したLSP(6)の1変数への変換，ドロネー三角形分割の利用が鍵になると思います）
関連記事：