Largest Small n-Polygons 生成プログラムの紹介

この記事では，Largest Small $n$-Polygons を生成する，以下のプログラムを紹介します．

https://github.com/tomtkg/Largest-Small-Polygon

Largest Small $n$-Polygons の詳細については，以下の記事をご確認ください．

https://qiita.com/tomtkg/items/55ccaa24c04c22d5baa9

紹介プログラムは，入力$n$に対して「いずれの2点間の距離も1以下で面積が最大となる$n$角形 ⇒ LSP($n$)を出力します．言語はMATLABです．

# LSP(n): nが奇数の場合
$n=$奇数の場合，LSP($n$)は正多角形（＝正$n$角形）になります．
プログラムでは，2種類の方法が書かれています．

1つ目は，非線形最適化ソルバー[fmincon](https://jp.mathworks.com/help/optim/ug/fmincon.html)を利用する方法です．
後述する$n=$偶数の場合と同じ方法で，最適なLSP($n$)の近似解を求めます．
動作方法は，`X = LSP(n)`のように，`LSP.m`ファイルを呼び出すだけです．

2つ目は，直径1の円周上に正$n$角形を作成し，それをスケーリング・移動する方法です．
`main.m`ファイルには，以下のように書かれています．

```MATLAB:正多角形変換法
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](https://jp.mathworks.com/help/stats/pdist.html) で計算し，その最大距離が1になるように図形をスケーリングします．
最後に，他の図形との比較を容易にするため，$y$座標の最大値が$0.5$になるように全体を移動します．

$n=5,7$で得られる図形は，以下の通りです．

|![](https://github.com/tomtkg/Largest-Small-Polygon/blob/main/image/special-5.png?raw=true)|![](https://github.com/tomtkg/Largest-Small-Polygon/blob/main/image/special-6.png?raw=true)|
|:-:|:-:|
|スケーリング前の正五角形 ($r=0.5$)<br>面積：$0.594410$|スケーリング前の正七角形 ($r=0.5$)<br>面積：$0.684103$|
|![](https://github.com/tomtkg/Largest-Small-Polygon/blob/main/image/Regular_polygon(5)-2.png?raw=true)|![](https://github.com/tomtkg/Largest-Small-Polygon/blob/main/image/Regular_polygon(7)-2.png?raw=true)|
|LSP(5)：正五角形 ($r=0.525731$)<br>面積：$0.657164$|LSP(7)：正七角形 ($r=0.512858$)<br>面積：$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](https://jp.mathworks.com/help/optim/ug/fmincon.html)を利用します．
そのため，[Optimization Toolbox](https://jp.mathworks.com/help/optim)が必要です．以下は，[LSP.m](https://github.com/tomtkg/Largest-Small-Polygon/blob/main/LSP.m)の全文です．

```MATLAB: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として，以下のように書けます．
```MATLAB:LSP(6) Graham's Biggest Little Hexagon with one variable
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$に設定した時の六角形は，以下になります．
|![](https://github.com/tomtkg/Largest-Small-Polygon/blob/main/image/GBLH-(0)-2.png?raw=true)|![](https://github.com/tomtkg/Largest-Small-Polygon/blob/main/image/GBLH-(1)-2.png?raw=true)|
|:-:|:-:|
|$x_1=0.0$，面積：$0.500000$|$x_1=0.1$，面積：$0.577995$|
|![](https://github.com/tomtkg/Largest-Small-Polygon/blob/main/image/GBLH-(2)-2.png?raw=true)|![](https://github.com/tomtkg/Largest-Small-Polygon/blob/main/image/GBLH-(3)-2.png?raw=true)|
|$x_1=0.2$，面積：$0.636767$|$x_1=0.3$，面積：$0.670788$|
|![](https://github.com/tomtkg/Largest-Small-Polygon/blob/main/image/GBLH-(4)-2.png?raw=true)|![](https://github.com/tomtkg/Largest-Small-Polygon/blob/main/image/GBLH-(5)-2.png?raw=true)|
|$x_1=0.4$，面積：$0.666007$|$x_1=0.5$，面積：$0.500000$|

$x_1$と六角形の面積$A$の関係は，以下になります．

<img src="https://github.com/tomtkg/Largest-Small-Polygon/blob/main/image/GBLH.png?raw=true" width="600">

図の横軸が$x_1$の値，縦軸が面積$A$の値を表します．
この図からも，$x_1=0.343771453...$の時に最大の面積$A=0.674981...$になることが直感的に理解できるかと思います．

上記のGraham's Biggest Little Hexagonに関する検証は，[関数 GBLH]() としてまとめられています．

# 終わりに
Largest Small $n$-Polygons の生成プログラムを紹介しました．
プログラムの行数に対して，説明を丁寧に書き過ぎたかもしれません．
LSPの説明と記事を分割すると決めてから，分量がかなり増大してしまいました．

紹介したプログラムですが，MATLABが優秀なのか，予想よりも簡潔に書けました．
[LSP.m](https://github.com/tomtkg/Largest-Small-Polygon/blob/main/LSP.m)は，幾つかの2点間の距離が$1$になる等式条件を追加したり，1点は必ず$(0.5,y)$の値を取る条件を追加したら，解の精度が向上すると考えられます．
是非，チューニングや新しい工夫を試して，[2021年の論文](https://arxiv.org/abs/2101.01263)と結果を比較してみてください．
どこかで工夫をしないと，LSP(1000)は解けないと思います．
（最後に紹介したLSP(6)の1変数への変換，[ドロネー三角形分割](https://jp.mathworks.com/help/matlab/ref/delaunaytriangulation.html)の利用が鍵になると思います）

関連記事：

https://github.com/tomtkg/Largest-Small-Polygon

https://qiita.com/tomtkg/items/55ccaa24c04c22d5baa9