数理最適化でアナグラム生成AIを作る 〜クラウドはXIII機関になれるか〜

[エイチームフィナジー](https://finergy.a-tm.co.jp/)のアドベントカレンダー1日目は、[s2terminal](https://qiita.com/suzuki_sh)が担当します！

※本稿はスクウェア・エニックス社が2005年に発売したプレイステーション2用ゲームソフト[「キングダムハーツ2」](https://ja.wikipedia.org/wiki/%E3%82%AD%E3%83%B3%E3%82%B0%E3%83%80%E3%83%A0_%E3%83%8F%E3%83%BC%E3%83%84II)について一部ネタバレを含みます。

----

現在私のチーム(というかこの問題に取り組んでいるのはほぼ私一人ですが)は、ウェブサイト内のコンテンツ掲載順序の最適化について研究開発しています。この記事では、そこで使った技術や使わなかった技術たちを応用して、**アナグラムの自動生成**に取り組んでみたいと思います。アナグラムの自動生成プログラムの開発を通して、最適化の手法や用語を少しずつ紹介していきます。最適化の手法がなにかの役に立つきっかけとなれば幸いです。

本稿では題材として、スクウェア・エニックス社から発売されているゲームシリーズ[「キングダムハーツ」シリーズ](https://www.jp.square-enix.com/kingdom/)に登場するXIII機関(じゅうさんきかん、Organization Thirteen)のような、**「元の名前に`X`(エックス)を足して並び替えた結果を、新しい名前とする」アナグラム**を題材とします。

たとえば、`sora`(ソラ)という文字列を受け取って、`roxas`(ロクサス)のような文字列を返す関数の作成を目標とします。

## 実行環境
- Python 3.8.6

またPythonの実行環境構築には下記を利用しています。

- Docker 19.03.13
- docker compose 1.27.4
- Poetry 1.1.4

Pythonプログラムの実行に必要なコマンドが本文記載の内容とは一部異なりますが、本文中では省略しています。例えば`$ python main.py`と表記している箇所は、実際には`$ docker-compose exec app poetry run python main.py`というコマンドを実行しています。
これらを含めた実際のプログラムはGitHubに公開しています。

https://github.com/s2terminal/xname

## 総当たりでのアナグラムの生成
さっそく「元の名前に`X`(エックス)を足して並び替えた結果を、新しい名前とする」アナグラムの自動生成をする関数を開発していきます。

まずはシンプルに、与えられた文字列のアナグラムを全パターン列挙するプログラムを書いていきます。これは文字列を配列化して、[itertools.permutations](https://docs.python.org/ja/3/library/itertools.html#itertools.permutations)を使えば簡単に実現できます。

```python
from itertools import permutations

original_name = "abe"
for x in permutations([ a for a in original_name ]):
    print("".join(x))
```

```
abe
aeb
bae
bea
eab
eba
```

取りうるアナグラムの列挙はできました。この中からひとつ選んで返す必要があります。そのために、これらのアナグラムのうちどの文字列が良いか、あるいは悪いのか、**評価**を定義していきます。

### 「名前らしさ」をあらわすアナグラムの評価関数の定義
たとえば`ROXAS`のアナグラムとして、ただ適当に並び替えると`SRXAO`といった読めない名前も登場してしまいます。そこで「文字列が"名前らしい"かどうか」を何らかの方法で判定する必要があります。

ここでは簡単のため、条件はただひとつだけ

- 子音が2つ以上連続していないこと

とします。また、母音も`a,e,i,o,u`の5文字のみとします。

より厳密には、長さ$n$の文字列$\\{x_i\\}$から隣接する2文字$x_i,x_{i+1}$を取って

- $x_i$と$x_{i+1}$のいずれも子音だった場合は、0
- $x_i$と$x_{i+1}$のいずれかが母音だった場合は、1

...とした数列(長さ$n-1$)の、平均値を評価関数とします。例えば`SRXAO`は0.5、`ROXAS`は1.0となります。

プログラムにすると以下のようになります。

```python:src/bruteforce.py
vowels = [ a for a in "aiueo" ]
def name_score(name: str, original: str) -> int:
    score = 0
    is_vowel = None
    for s in name:
        if (is_vowel is not None and (is_vowel or (s in vowels))): score += 1
        is_vowel = (s in vowels)
    score /= (len(name) - 1)
    return score
```

「この値が最大になる文字列が、アナグラムとして適切である」として、探索していきます。

### 探索の実行
先述の評価関数部分を除いた、プログラムの全体像は下記のようになります。

```python:src/bruteforce.py
def search_xname(original_name: str):
    best_score = 0
    best_name = original_name
    for v in permutations([ n for n in original_name ]):
        name = "".join(v)
        score = name_score(name, original_name)
        if (best_score < score):
            best_score = score
            best_name = name
    return (best_name, best_score)

if __name__ == '__main__':
    original_name = input('input name: ')
    max_names, score = search_xname(original_name)
    print("{}は{}になりました。スコア: {}".format(original_name, max_names, score))
```

実行してみましょう。

```bash
$ python src/bruteforce.py
input name: sora
soraはsoraxになりました。スコア: 1.0
```

できました！`sora`に`x`を足したアナグラムは`sorax`(ソーラックス)になりました。

このように、ある条件を満たす解候補の中から最も最適な解を見つけることを**最適化**と言います。特に、最適化の対象や解候補の条件が数式で表されてる最適化を**数理最適化**と言います。[^optimise] このプログラムでいう`name_score()`のように、数理最適化において最適化対象となる関数を**目的関数**と呼びます。

[^optimise]: https://www.slideshare.net/amp-kyoto/it-59618687

### レーベンシュタイン距離による目的関数へのペナルティの導入で、元の名前と"似ていない"名前を探索する
しかし。`sora`に`x`を足して探索したアナグラムが`sorax`というのは、「末尾に文字を付けただけでは...」となって、なんだか**"面白み"に欠ける結果**です。

なぜ面白くないのでしょうか。ここでは「元の名前と似ているから面白くない」と仮定します。もう少し面白いアナグラムになるよう、プログラムを改良してみましょう。

そこで「元の名前と似ていないアナグラム」を探すために、「元の名前とどれくらい似ているか」を定量化します。ここでは**レーベンシュタイン距離**を導入します。

レーベンシュタイン距離とは、元の文字列から挿入・削除・置換・交換といった操作を何回行えば目的の文字列を得られるか、という指標です。「距離」とあるだけあって、多ければ多いほど遠いです。例えば「sora」と「sorax」のレーベンシュタイン距離は、「x」を挿入する1回の操作だけで目的の文字列が得られるため、1です。

ここでは、文字列の長さが変わっても同じ指標で比べられるように、レーベンシュタイン距離を目的の文字列の長さで割ることで0〜1の値を取るように標準化した値を用います。

各文字列のレーベンシュタイン距離の例は以下です。

元の文字列|目的の文字列|レーベンシュタイン距離|標準化レーベンシュタイン距離|元の文字と似てるように見えるか
:--|:--|:--|:--|:--
sora|sorax|1|0.2|似てる
sora|roxas|3|0.6|似てない
sora|srxao|3|0.6|似てない
ansem|xansem|1|0.16666666666666666|似てる
ansem|xemnas|6|1.0|似てない

標準化したレーベンシュタイン距離の小ささをペナルティとして与えた改善評価関数は以下のようになります。事前に[python\-Levenshtein](https://pypi.org/project/python-Levenshtein/)ライブラリをインストール・インポートしておく必要があります。

```diff:src/bruteforce.py
 vowels = [ a for a in "aiueo" ]
 def name_score(name: str, original: str) -> int:
     score = 0
     is_vowel = None
     for s in name:
         if (is_vowel is not None and (is_vowel or (s in vowels))): score += 1
         is_vowel = (s in vowels)
     score /= (len(name) - 1)
+    score -= (1 - (Levenshtein.distance(name, original)) / (len(name)))
     return score
```

```bash
$ python src/bruteforce.py
input name: Sora
SoraはSaxorになりました。スコア: 0.8
```

`Saxor`(サクソー)という名前が得られました。元のプログラムで得られた`Sorax`よりは、結構"面白い"結果だと言えるのではないでしょうか。

このように目的関数にペナルティを与える等で調整していくと、最適化によって得られる結果がより求められる結果に近づいていきます。

### この方法での問題点
さあ、これでどんな名前でもアナグラムが生成できるようになりました。

試しに、スクウェア・エニックスから発売されている[ゲーム「ファイナルファンタジーVII」](https://www.jp.square-enix.com/ffvii_remake/)に登場する架空のキャラクター「クラウド(Cloud)」の名前を入力してみましょう。

```bash
$ python src/bruteforce.py
input name: cloud
cloudはluxcodになりました。スコア: 0.8
```

おっと、いけません。[クラウドの本名はクラウドストライフ(Cloud Strife)](https://www.jp.square-enix.com/ffvii_remake/character/index.html)と言うのでした。こちらで入力してみましょう。

```bash
$ python src/bruteforce.py
input name: cloudstrife



```

...なかなか結果が返ってこないと思います。

それでは、バンダイナムコゲームスから発売されている[「エースコンバット7」](https://ace7.acecombat.jp/)に登場するライバルキャラクター「ミハイ」の本名**「Mihaly Dumitru Margareta Corneliu Leopold Blanca Karol Aeon Ignatius Raphael Maria Niketas A Shilage（ミハイドゥミトルマルガレータコルネリウレオポルドブランカカロルイオンイグナチウスラファエルマリアニケタスアシラージ）」**の名前を入力してみましょう。

```bash
$ python src/bruteforce.py
input name: mihalydumitrumargaretacorneliuleopoldblancakarolaeonignatiusraphaelmarianiketasashilage


```

......。

............。

........................結果は全然返って来ないと思います。

このプログラムは**順列**を探索しています。順列というのは、要素の数が増えるほどパターン数は爆発的に多くなります。以下、パターン数と私のPCでのプログラム実行時間を記します。

文字列|パターン数|実行時間
--:|--:|--:
sora + x|5!=120|0.0018秒
ansem + x|6!=720|0.0133秒
lauriam + x|8!=40,320|0.9546秒
elizabeth + x|10!=3,628,800|145.5917秒
cloudstrife + x|12!=479,001,600|--
mihaly(中略)shilage + x|88!≒$10^{134}$|--

「mihalydumitrumargaretacorneliuleopoldblancakarolaeonignatiusraphaelmarianiketasashilage」という名前にxを足したアナグラムの探索は、全部で**88!=185482642257398439114796845645546284380220968949399346684421580986889562184028199319100141244804501828416633516851200000000000000000000通り**の探索となります。これはとても一般的なコンピュータで全部の数値を計算することはできない量です。このプログラムで解を求められるのは、高々10文字くらいの名前が限界だと思います。

このように長い文字列のアナグラムにおいて、今まで述べたような総当たりで厳密な最適解を求めるのは現実的ではありません。このような時でも"だいたいそれらしい"解が得られるように、**厳密解**ではなく**近似解**、最大値ではないかもしれない(最大値かもしれない)解を得る方法を見ていきます。

## ランダム法による近似解法
解になりうるパターンから適当にいくつか選んで目的関数を計算し、最大値となるパターンを選んでみます。ここで「順列」の中から適当にパターンを選ぶ必要があるのですが、まず「文字数と同じ数の乱数を生成し、昇順に並び替える」方法で、ランダムな順列を得ていきます。各文字に重み付けをしていくようなイメージです。乱数の列をいくつか**探索**し、その中で最適な物を選んでいきます。

たとえば元の文字列が`sorax`であれば、乱数を5個生成し、「5」「2」「1」「4」「3」となれば`r`が1番目、`o`が2番目...となって最終的に`roxas`という文字列が解候補として得られます。

なお乱数で同じ数字が出た時の挙動を一意にするために、長さ$n$の文字列の$i$番目の文字に対して、0から$n-1$までの乱数$r(n)$を得たときに、$r(n) + i/n$で重み付けをしていきます。

```python:random_search.py
import random

def name_generator(name: str):
    size = len(name)
    x = [random.randint(0, size - 1) + i/size for i in range(0, size)]
    name_list = list(zip([a for a in name], x))
    return "".join([x[0] for x in sorted(name_list, key=lambda x: x[1])])
```

これを適当な回数繰り返して探索し、最もいいスコアになった物を最適解として返します。

```python:random_search.py
from typing import Callable, List

def random_search_base(name: str, name_generator: Callable[[str], str], times: int):
    original_name = name
    best_score = -float('inf')
    best_name = name
    for _ in range(0, times):
        tmp_name = name_generator(name)
        tmp_score = name_score(tmp_name, original_name)
        if (best_score < tmp_score):
            best_score = tmp_score
            best_name = tmp_name
    return best_name, best_score

def search_ary(name: str, times=10):
    return random_search_base(name, name_generator, times)
```

```bash
$ python src/random_search.py 
input name: cloudstrife
('lscetfdorxui', 0.46212121212121204)
```

`cloudstrife`に`x`を足したアナグラムとして`lscetfdorxui`(ルスセトフドークシー)という結果が得られました。探索回数にもよりますが、1,000回程度であればすぐに結果が返ってきます。(記載したプログラムは10回)

### 実行可能領域への写像を調整する
ここで言う「与えられた文字列のアナグラム」のような、解の候補となる条件を満たした要素のことを数理最適化の用語で**実行可能解**と呼びます。実行可能解の集合を**実行可能領域**と呼びます。

このプログラムでの探索空間である「n個の乱数の集合」というのは、「実行可能領域」と1対1で対応しておらず、重複があります。例えば乱数のうち「1,2,3,4,5」も「1,1,1,2,3」も、同じ実行可能解になってしまいます。厳密に言うと、探索空間から実行可能領域への写像が、全射ではありますが単射になっていません。単射になるような写像を定義することで、探索の効率を変えてみましょう。

いままでは$n$通りの乱数を$n$個受け取っていました。
代わりに$n!$通りの乱数をひとつだけ受け取って、$n$で割った余り(0〜$n-1$の$n$通り)をキーとする配列の要素を抽出し、$n$で割った商について同じ手続きを繰り返す、という関数を定義してみます。文字で説明するとなんだかややこしいですが、プログラムにすると下記です。

```python
def func(ary: typing.List, i: int) -> typing.List:
    n = len(ary)
    if n <= 1:
        return ary
    val = ary.pop(i % n)
    ret = func(ary, i // n)
    ret.extend([val])
    return ret
```

この関数はこんな挙動をします。

```python
l = [1,2,3,4]
print(func(l, 0))
print(func(l, 1))
print(func(l, 2))
print(func(l, 3))
print(func(l, 4))
```

```
python src/opt.py 
[4, 3, 2, 1]
[1, 2, 4, 3]
[3, 2, 1, 4]
[1, 2, 3, 4]
[4, 2, 3, 1]
```

写像$f:N^1\rightarrow N^n$をひとつ挟むことで、探索空間をn次元から1次元に削減しました。プログラムで言うと、生成する乱数がn個から1個に減り、単射になっているので効率が良くなります。

```python:random_search.py
def search_map(name: str, times=DEFAULT_SEARCH_TIMES):
    name_generator = lambda name: "".join(
        func([a for a in name], random.randint(0, math.factorial(len(name)) - 1))
    )
    return random_search_base(name, name_generator, times)
```

実行してみます。

```bash
$ python src/random_search.py 
input name: cloudstrife
('tslrxoiceduf', 0.6363636363636364)
```

`tslrxoiceduf`(ツルルクソイセダフ)という結果が得られました。


長い文字列でもちゃんと結果が返ってきます。

```bash
$ python src/opt.py 
input name: mihalydumitrumargaretacorneliuleopoldblancakarolaeonignatiusraphaelmarianiketasashilage
('culrlruetailiolkgleapyixanaaphoittnaurimrbsheanlasohamlagsoraritmacdakuenraenelgodmeaaii', 0.5993991640543365)
```

「カルルルルータイリアルクグリーピクサナーフォイトノーリムルブシーンラソハムラグソーラーリトマクダクーンラエネルゴドミーエーイー」という結果が得られました。

このように実行可能領域とは異なる探索空間を定義し、実行可能領域への写像を定義することで探索ができるようになります。この方法を取った場合は、写像や探索空間の定義を工夫することで、最適化の性能を調整できます。実際には**Bottom-Left法[^bottom-left]**のように、実行可能領域が広大で探索が大変な場合等に、より"狭い"探索空間を探索する方法として使われています。

[^bottom-left]: http://www.msi.co.jp/nuopt/glossary/term_4497146a690870abd1038155739dc60e78d8f20c.html

また本稿では乱数を使いましたが、たとえば乱数が1と2とでそれぞれ評価関数にかけた場合に"近い"スコアが結果として得られる場合など、探索空間上の評価関数値が何らかの法則(関数)に沿って分布している場合は、乱数ではなく**ベイズ最適化(Bayesian Optimization)[^bayesian-optimization]**など別のブラックボックス最適化の手法を使うことで、より効率的に探索できると思います。

[^bayesian-optimization]: https://tech.preferred.jp/ja/blog/high-dimensional-black-box-optimization-using-bayesian-optimization/

## 交換近傍による局所探索法
いままでは「実行可能領域は未知のものである」として、乱数列を探索空間として実行可能領域への写像を定義することで、最適化できるようにしてきました。ですが本稿で題材としているのはアナグラムであり、列挙可能な有限個の順列からなる集合です。この実行可能領域の特徴を踏まえて直接探索する方法を紹介します。

アナグラムの探索において、元の文字列に対して「元の文字列をちょっとだけ変えた文字列」を探索していって、目的関数に改善が見られたらその文字列に移動する、という方法が取れます。例えば最初の出発点が`sorax`であれば、`sorax`の目的関数のスコアは0.0ですが、最初と最後の文字だけ入れ替えた`xoras`という文字列のスコアは0.4です。元の文字列をちょっとだけ変えた文字列の探索を続けていくことで、実行可能領域が広大だったとしてもより良い解を探索できます。

このように、ある実行可能解を少しだけ変更した解のことを**近傍(Neighborhood)**と呼び、解候補の近傍の中からより良い解があれば移動する、という手順を繰り返していく最適化手法のことを**局所探索法(Local Search)**と呼びます。

まずは「元の文字列から、任意の2文字を選択して入れ替えたもの」を近傍として定義しましょう。このような近傍を**交換近傍(Swap-based Neighborhood)**と呼びます。Pythonのプログラムで書くと、以下のようなイテレータで表現できます。

```python:local_search.py
def swap_two_any_neighbor(name: str) -> Iterator[int]:
    size = len(name)
    for i in range(0, size - 1):
        for j in range(i + 1, size):
            yield (i, j)
```

`swap_two_any_neighbor()`による近傍を一通り探索するメソッド`swap_two_base()`を定義します。

```python:local_search.py
def swap_two_base(local_best_name: str, original_name: str):
    best_score = -float('inf')
    best_name = local_best_name
    for (i, j) in swap_two_any_neighbor(local_best_name):
        tmp_name = swap_name(local_best_name, [i,j])
        tmp_score = name_score(tmp_name, original_name)
        if (best_score < tmp_score):
            best_score = tmp_score
            best_name = tmp_name
    return (best_name, best_score)

def swap_name(name: str, swaps: Tuple[int, int]) -> str:
    name_ary = [a for a in name]
    ret = ""
    for i in range(len(name)):
        key = i
        if   i == swaps[0]: key = swaps[1]
        elif i == swaps[1]: key = swaps[0]
        ret += name_ary[key]
    return ret
```

そして、`swap_two_base()`による手続きを適当な回数(一般的には、目的関数の改善が見られなくなるまで)繰り返します。

```python:local_search.py
def local_search_base(original_name: str, times=None):
    if times is None: times = len(original_name)
    best_score = -float('inf')
    best_name = original_name
    for _ in range(times):
        (tmp_name, tmp_score) = swap_two_base(best_name, original_name)
        if (best_score < tmp_score):
            best_score = tmp_score
            best_name = tmp_name
        else:
            break
    return (best_name, best_score)
```

実行してみましょう。

```
$ python src/local_search.py 
input name: cloudstrife
('dilecufxotsr', 0.7272727272727273)
```

`dilecufxotsr`(ディレカフクソツ)という結果が得られました。

### 隣接交換近傍による局所探索法
先のプログラムでは近傍の定義を「任意の2文字を入れ替える」としましたが、ここを変更して「隣接した2文字を入れ替える」としてみます。このような近傍を**隣接交換近傍**と呼びます。プログラムで表現すると下記です。

```python:local_search.py
def swap_two_adjacent_neighbor(name: str) -> Iterator[int]:
    size = len(name)
    for i in range(0, size - 1):
        yield (i, i + 1)
```

実行してみましょう。

```
$ python src/local_search.py 
input name: cloudstrife
('xolusdtirefx', 0.31060606060606066)
```

`colusdtirefx`(コラスドティレフクス)という結果が得られました。なんとなく元の文字列に似ているように見えますし、スコアも悪くなってしまっています。

このように局所探索法において近傍の定義は複数あり、差し替えることで最適化の結果も変わってきます。本稿で紹介した交換近傍の他に**挿入近傍**や**2-opt近傍**といった物があり、問題によって最適な近傍も変わってきます。

### 局所探索法のスタート点
実際に局所探索法を行う時は、出発点としてそれなりに"良い"解候補から探索を始める事が大切です。しかしこのプログラムでは出発点が元の文字列、つまりレーベンシュタイン距離が最も小さく、スコアの悪い出発点になっています。例えば以下のように、元の文字列を逆順にした文字列を出発点とするような工夫をすると、より良い結果が得られます。

```diff:local_search.py
 def local_search_base(original_name: str, neighbor_generator: Callable[[int], Iterator[int]], times=None):
     if times is None: times = len(original_name)
     best_score = -float('inf')
-    best_name = original_name
+    best_name = ''.join(list(reversed(original_name)))
````

この方法の応用として、複数の出発点から局所探索法で解候補を探索していく**多スタート局所探索法(multi-start local search)**という手法もあります。

## 線形計画問題としての方法
最後に[PuLP](https://pypi.org/project/PuLP/)というライブラリを用いて、**[巡回セールスマン問題(Travelling Salesman Problem、TSP)](https://ja.wikipedia.org/wiki/%E5%B7%A1%E5%9B%9E%E3%82%BB%E3%83%BC%E3%83%AB%E3%82%B9%E3%83%9E%E3%83%B3%E5%95%8F%E9%A1%8C)**の解き方を応用してみます。本稿で取り上げるアナグラム探索問題はTSPではないためTSPそのものについては説明しませんが、少し知っていたほうが本稿は読みやすいと思います。

数理最適化問題のうち、目的関数と制約条件とが線形(1次式。$c_1x_1+c_2x_2+...$のような形式)で表現できる問題を**線形計画問題**と言います。PuLPはざっくりいうと、**線形計画問題(Linear Programming、LP)**を解くライブラリ[^pulp]です。軽くPuLP自体の紹介をします。

[^pulp]: 厳密にはPuLPは線形計画問題ソルバーのラッパーライブラリであり、デフォルトでCBCというソルバーが設定されています。

以下は、PuLPのREADMEに載っている線形計画問題の例です。

```math
\begin{eqnarray}
\text{Minimize}& &-4x+y \\
\text{Subjet to}& &x+y \leq 2\\
& & 0 \leq x \leq 3 \\
& & 0 \leq y \leq 1 \\
\end{eqnarray}
```

これは、変数$x$は0から3、変数$y$は0から1の値を取るとして、目的関数$x+y$が2以下になるような$(x,y)$の中で、$-4x + y$を最小化するような$(x,y)$を求める、という問題の定式化です。ちなみに答えは$x=2,y=0$です。

この問題をPuLPを使って解くと、以下のようになります。

```python
import pulp

# 変数x,yの定義
x = pulp.LpVariable("x", 0, 3)
y = pulp.LpVariable("y", 0, 1)
# 最小化(Minimize)問題オブジェクトを生成
prob = pulp.LpProblem("myProblem", pulp.LpMinimize)
# 制約条件の追加 (問題オブジェクトに式を+=すると条件を追加できる)
prob += x + y <= 2
# 目的関数の追加 (等号・不等号の無い式を+=すると制約条件ではなく目的関数の追加となる)
prob += -4*x + y
# 解の探索
status = prob.solve()

# 結果の出力
pulp.LpStatus[status], x.value(), y.value()
```

`('Optimal', 2.0, 0.0)`という出力が得られます。Optimal(status=1)というのは、厳密解が得られたという事です。厳密解ではなく近似解が得られた場合は別の出力が得られます。

PuLP自体はそこまで複雑ではないものの、通常のプログラミングとはちょっと変わった使い方のライブラリなので戸惑う人もいるかもしれません(筆者がそうでした)。PuLPの詳しい使い方については、[PuLPによるモデル作成方法 — Pythonオンライン学習サービス PyQ（パイキュー）](https://docs.pyq.jp/python/math_opt/pulp.html)等を参考にしてください。

さて、アナグラムの問題に話を戻します。結論から言うと、本稿で取り上げているアナグラムの探索問題を少し簡略化したものが、以下のように一次式で**定式化**できます。

```math
\begin{eqnarray}
&\text{Maximize}\ & CX&&\\
&\text{Subject to}\ &
   \sum_{i\in V} \sum_{j\in V} x_{ij} = n - 1\\
&& \sum_{j\in V} x_{ij}\leq 1 ,\:\sum_{j\in V} x_{ji}\leq 1\:& &(\forall i \in V)\\
&& \sum_{j\in V} x_{ij} + \sum_{j\in V} x_{ji} \geq 1\:& &(\forall i \in V) \\
&& x_{ii} = 0\:& &(\forall i \in V)\\
\end{eqnarray}\\
```

順に説明していきます。

### バイナリ変数から文字列を得る
数理最適化の定式化において、変数をゼロかイチかの**バイナリ変数**で定義すると、表現力が高まります。そこで、アナグラムにおいて「ある文字の次がある文字である」時に1、そうでない時に0、を取るような変数を使って、この問題を定式化してみます。

$n$文字のアナグラムにおいて、ゼロかイチかを取る変数$x$を$n^2$個用意します。元の文字列の$i$文字目の次が、元の文字列の$j$文字目であるときに$x_{ij}=1$、そうでないときは$x_{ij}=0$となるように定義します。するとアナグラムを探索する問題は、行列$X = \\{ x_{ij} \mid x_{ij} \in \\{0,1\\} \\}$を探索する問題であると捉えることができます。

例えば元の文字列として`sorax`という文字列があって、アナグラムの最初は`ro`から始まるとしたら、`r`は3文字目、`o`は2文字目なので、$x_{32}=1$というように定義していきます。

```math
x_{31}=0\\
x_{32}=1\\
x_{33}=0\\
x_{34}=0\\
x_{35}=0\\
```

つまり

```math
x_{3j} = \begin{pmatrix}
0 & 1 & 0 & 0 & 0
\end{pmatrix}
```

同様にして、soraxからroxasというアナグラムを得ることを表す$x$は、以下のように行列で定義できます。

```math
X = \{x_{ij}\} = \begin{pmatrix}
0 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & 1\\
0 & 1 & 0 & 0 & 0\\
1 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 1 & 0\\
\end{pmatrix}
```

このような名前のバイナリ行列$X$を探索する問題であると定義することができます。

行列の扱いはプログラムでは簡単です。今回は多次元Listを用いて定義します。名前のバイナリ行列から実際の名前を得る関数は下記のように定義できます。

```python:pulp_search.py
def target_name(xs, name: str):
    size = len(name)
    start = 0
    # 最初の1文字を探す
    for j in range(size):
        if (all([xs[i][j] == 0 for i in range(size)])):
            start = j

    ret = name[start]
    current = start
    for _ in range(size - 1):
        for j in range(size):
            if xs[current][j] == 1:
                ret += name[j]
                current = j
                break
    return ret
```

### 目的関数の定式化
名前のバイナリ行列と同様の手順で、元の文字列に対して母音が連続している箇所を示す下記のような定数行列を導入します。

```python
cs = []
for i in range(size):
    cs.append([])
    for j in range(size):
        if name[i] in vowels or name[j] in vowels:
            cs[i].append(1)
        else:
            cs[i].append(0)
```

この$X$と$C$の内積を最大化する問題である、と言えます。

```math
\text{Maximize}\ CX\\
```


これをPulPプログラムで表すと、下記のようになります。

```python
prob = pulp.LpProblem('myprob', sense=pulp.LpMaximize)
prob += pulp.lpSum([pulp.lpDot(c, x) for c, x in zip(cs, xs)])
```

### 制約条件の定義
定義したバイナリ文字列は制約が無さすぎます。今回のアナグラム問題に適用できるよう、制約を追加していきます。

まず「探索する対象はn文字の文字列である」という条件は「文字の連続がn-1回だけ起きる」と言えるため、つまりイチになっている箇所がn-1個だけあれば良いわけです。数式だと以下のように表現できます。

なお、以下$V$を1から$n$までの整数の集合$V=\\{ i \mid i \in \mathbb{N}, 1 \leq i \leq n \\}$とします。

```math
\sum_{i\in V} \sum_{j\in V} x_{ij} = n - 1\\
```

```python
prob += pulp.lpSum([pulp.lpSum([x2 for x2 in x1]) for x1 in xs]) == size - 1
```

「各文字への行きと帰りは、それぞれ一度しか選ばれない」という条件は、「バイナリ行列$X$のどの行または列も、合計値が1以下」と言えます。(「1になる」ではなく「1以下」としているのは、「最初の文字へ帰る文字」と「最後の文字から行く文字」が無いからです)

```math
\sum_{j\in V} x_{ij}\leq 1 ,\:\sum_{j\in V} x_{ji}\leq 1\: (\forall i \in V)\\
```

プログラムにすると以下のように書けます。

```python
for i in range(size):
    prob += pulp.lpSum([xs[j][i] for j in range(size)]) <= 1
    prob += pulp.lpSum([xs[i][j] for j in range(size)]) <= 1
```

また、「すべての文字は、行きか帰りに最低一度は選ばれる」という制約が必要です。

```math
\sum_{j\in V} x_{ij} + \sum_{j\in V} x_{ji} \geq 1\: (\forall i \in V) \\
```

```python
for i in range(size):
    prob += pulp.lpSum([xs[i][j] for j in range(size)]) + pulp.lpSum([xs[j][i] for j in range(size)]) >= 1
```

「同じ文字には連続しない」という制約は下記です。

```math
x_{ii} = 0\: (\forall i \in V)
```

```python
for i in range(size):
    prob += xs[i][i] == 0
```

また、**部分順回路(subtour)**を排除する必要があります。たとえばここまでの制約だけだと、`cloudstrifex`という文字列のアナグラムを探索した時に`x`→`l`→`r`→`x`→...のように、すべての文字列を回ること無く一部だけでループした結果が得られてしまいます。このような結果の対策となる制約条件を追加します。詳細は、TSPの例ですが[この記事](http://www.ie110704.net/2020/08/15/pulp%E3%81%A7%E6%95%B0%E7%90%86%E6%9C%80%E9%81%A9%E5%8C%96%E5%95%8F%E9%A1%8C%EF%BC%88tsp%E3%80%81vrp%EF%BC%89/)や[この記事](https://qiita.com/panchovie/items/6509fb54e3d53f4766aa)を読んでください。なおTSPのようなハミルトン閉路では$i,j = 0$のとき制約を追加しないのですが、今回アナグラムの探索は閉路ではない(最後の文字と最初の文字が繋がっているわけではない)ので、そのような対応はしていません。

```python
    for i in range(size):
        for j in range(size):
            if i != j:
                prob += u[i] - u[j] <= (1 - xs[i][j]) * size - 1
```

以上で準備完了です。実行してみましょう。

```bash
$ run python src/pulp_search.py
(中略)
('xrdsucelifot', 0.5606060606060607)
```

`xrdsucelifot`(クスルドスーセリフォット)という結果が得られました！

...結構大変だった割に、あまり"スコア"が高くないですね。ここまで読んだ方は気づいたかもしれませんが、今までの手法と異なりPuLPでの探索時にはレーベンシュタイン距離によるペナルティを導入できていません。このように、現実の問題はすべて一次式で表現できるわけではないのですが、一次式で表現できない制約は盛り込むことができない点が線形計画法のライブラリを使う際の弱点になります。レーベンシュタイン距離ではなく一次式で表現できるような何らかの距離を代替として導入できれば、改善できると思います。

## まとめ
`cloudstrife`に`x`を加えた文字列から得られたアナグラムをまとめると下記になります。どの結果が気に入りましたか？

手法|得られたアナグラム|(よみ)|スコア
:--|:--|:--|:--|:--
総当たり|-|-|-
N次元探索空間ランダム|Lscetfdorxui|ルスセトフドークシー|0.46212121212121204
1次元探索空間ランダム|Tslrxoiceduf|ツルルクソイセダフ|0.6363636363636364
交換近傍による局所探索|Dilecufxotsr|ディレカフクソツ|0.7272727272727273
隣接交換近傍による局所探索|Colusdtirefx|コラスドティレフクス|0.31060606060606066
PuLP|Xrdsucelifot|クスルドスーセリフォット|0.5606060606060607

本稿では、要素の並び替え問題をプログラムによって自動で解く方法を紹介しました。「ランダム法」「局所探索法」「線形計画法」を利用してアナグラムを得ることができ、またそれらの改善余地について明らかにしました。本稿ではアナグラムを題材としましたが、順列を最適化するような問題にはこれら技術の一部または全部を応用することができると思います。

冒頭に述べたとおり、私のチームではWebサイト上の掲載コンテンツ表示順序を最適化するためにこの技術を研究しており、本稿で最初から述べてきたのとほぼ同じ手順で調査・開発してきたものになります。そのため本稿自体はアナグラムの探索やTSPの解法としては非効率な手順も紹介していますが、最適化手法を一部でもなにかに応用するような時に、参考になればと思います。

もし実際にアナグラムを自動生成するときは、単なる順列ではなく要素に重複がある問題となるため、工夫をすることで探索空間を大幅に削減できると思います。また、現行のプログラムでは`Marluxia`や`Xemnas`といった文字列は正しい名前である[^organizationxiii]ものの目的関数による評価値が最大になりません。「名前らしさ」の定義をより厳密にすることで、より良いアナグラムが得られるプログラムになると思います。

[^organizationxiii]: https://kingdomhearts.fandom.com/wiki/Organization_XIII

## 参考文献
- [朝倉書店｜Pythonによる 数理最適化入門](http://www.asakura.co.jp/books/isbn/978-4-254-12895-6/)
- [最適化手法入門 \| 書籍情報 \| 株式会社 講談社サイエンティフィク](https://www.kspub.co.jp/book/detail/5170083.html)
- [局所探索法とその拡張—タブー探索法を中心として Local Search and its Extensions—Focusing on Tabu Search  野々部 宏司＊ 柳浦 睦憲＊＊](http://www.co.mi.i.nagoya-u.ac.jp/~yagiura/papers/tabu-kaisetsu-2008.pdf)
- [線形順序付け問題の解法  櫻庭　セルソ　智，柳浦　睦憲](http://www.orsj.or.jp/archive2/or57-06/or57_6_327.pdf)
- [組み合わせ最適化問題に対する近似解法](http://dse.ssi.ist.hokudai.ac.jp/~onosato/lectures/Opt19/sec09.pdf)
- [PuLP による数理最適化超入門](http://www.nct9.ne.jp/m_hiroi/light/pulp01.html#chap02)
- [巡回セールスマン問題 \| Code Craft House](http://codecrafthouse.jp/p/2016/08/traveling-salesman-problem/)
- [XIII機関 \- Wikipedia](https://ja.wikipedia.org/wiki/XIII%E6%A9%9F%E9%96%A2)
- [エースコンバット7 スカイズ・アンノウン \- Wikipedia](https://ja.wikipedia.org/wiki/%E3%82%A8%E3%83%BC%E3%82%B9%E3%82%B3%E3%83%B3%E3%83%90%E3%83%83%E3%83%887_%E3%82%B9%E3%82%AB%E3%82%A4%E3%82%BA%E3%83%BB%E3%82%A2%E3%83%B3%E3%83%8E%E3%82%A6%E3%83%B3#%E3%82%A8%E3%83%AB%E3%82%B8%E3%82%A2_2)

また、本文中の架空の英語文字列のカタカナ表記は[英語→カタカナ変換機](https://www.sljfaq.org/cgi/e2k_ja.cgi)によって機械的に生成しました。