2021-12-21

Googleの量子超越性をスパコンで300秒で計算した論文を解説してみる

2021-12-03

OpenJijを使ったスピングラスの数値計算

この記事は2021年物理学アドベントカレンダーの3日目です。

2021年のアドベントカレンダーということなので2021年のノーベル物理学賞の一つのスピングラスを題材にしたいと思います。

さっそく本題のOpenJijを使った数値計算について知りたい方はここに飛んでください。

パリージさんはスピングラスという模型の解析手法についての功績が評価されての受賞となりました。スピングラスは物理の統計力学における複雑な磁性体の数理モデルを表す系ですが、その後この解析手法は今の深層学習の元となるニューラルネットの解析など様々な複雑系の解析への応用がなされており物理を飛び出して非常に多くの分野に影響を与えています。

日本物理学会のノーベル賞解説に大関先生が解説を書かれているのでぜひ読んでみてください（https://www.jps.or.jp/public/2021nobel1.php）。

この記事ではそのスピングラスの簡単な解説と手軽に使えるOpenJijというイジングモデルの数値計算ツールを使って実際にスピングラスの数値計算をしてノーベル賞の公式ページにあるPDF(https://www.nobelprize.org/uploads/2021/10/sciback_fy_en_21.pdf)にある図を再現することを目標にします。

読者には学部の統計力学でイジングモデルについて一応学んだことがあるくらいを想定します。

イジングモデルと相転移

簡単にイジングモデルについて紹介しておきます。

イジングモデルは磁石を表すモデルで小さな磁石（スピン）の集まりを考えてシンプルな相互作用のみを考えます。

$$E = -\sum_{i, j}J_{i, j}\sigma_i \sigma_j,~\sigma_i \in {-1, 1}$$

イジングモデルのエネルギーを見るとわかるように相互作用係数$J_{ij}$が正の場合は各スピン同士が同符号の場合がエネルギーが下がります。この$J_{ij}$は結合しているスピンの揃いやすさを表しています。すべての$J_{ij}$が正の場合、ボルツマン分布

$$p(\{\sigma_i\}_{i=1}^{N}) = \exp\left(-\beta E(\{\sigma_i\}_{i=1}^{N})\right)/Z$$

に従ってスピン配列$\{\sigma_i\}$を生成して、その各サンプルを使って以下の磁化

$$m = \frac{1}{N}\left\langle \sum_i\sigma_i\right\rangle$$

を温度$T$を横軸にプロットするとある温度より低いところで急激に$m$が有限の値を持つようになります。これが磁性体の相転移を表しています。$m$は各スピンの揃い具合を表していて、温度が低いところでは多くのスピンが同じ向きを向くことで大きな磁化を作ります。

スピングラス

では$J_{i, j}$がもっと複雑な値を持っていたらどうでしょうか。特に$i, j$によって全く異なる値を持っている時にうまく統計力学的な解析によって相転移のようなマクロな物理現象を見出すことはできるでしょうか。

簡単な強磁性イジングモデル（$J_{ij}$がすべて正）の場合は系全体を特徴づける量として磁化が導入されましたが、$J_{i, j}$がバラバラな値を持っていると、エネルギーの低いスピン配列も同様にバラバラになっていると考えられます。そのため磁化のように単純にスピンの値を足し算してしまうと磁化は0になってしまうため温度が高い時と温度が低い時で区別がつきません。

このため$J_{i, j}$がランダムな分布の場合は磁化ではなく、別の特徴づけが必要となります。その特徴を見出す解析手法がレプリカ法です。

レプリカ法

系の自由エネルギーを解析すれば、各温度での系の安定な状態を見つけることができます。自由エネルギーは分配関数の対数を取れば計算できるのでした。なので分配関数の対数を計算することを考えます。これはランダムな相互作用の時でも同じです。

しかし対数を取ってもこれ以上計算できないので対数と極限に関する以下の恒等式を用います。

$$\ln Z = \lim_{n\rightarrow 0}\frac{Z^n-1}{n}$$

こうするうことで分配関数の対数ではなく$Z^n$を計算すれば良いことになりました。$Z^n$は$\ln Z$に比べると計算が楽になるのでこのまま$Z^n$を計算していきます。

$n$は極限を取っているので実数なのですが、一旦整数だと思うことにします。そうすると$Z^n$は

$$Z^n = \left(\sum_{\sigma^{(1)}} \exp \left[-\beta H(\sigma^{(1)})\right]\right) \left(\sum_{\sigma^{(2)}} \exp\left[-\beta H(\sigma^{(2)})\right]\right) \cdots $$

という風にかけて独立な同じ系を$n$個用意した系の分配関数とみなすことができます。このように$Z^n$の$n$をいったん整数だとみなして$Z^n$を計算してその結果の$n$を$n \rightarrow 0$へと外挿します。これをレプリカ法と呼びます。

急に$n$を整数だと思ってその後外挿するなんて乱暴だと思われるかもしれませんが、これが実際これまで多くの有用な結果を導きだしてきています。こちらの数学的な正当性についてや条件などについては詳しくないので別の文献を参照していただければと思います。

このあとレプリカ法では上記の手続きにおいて$n$を整数だと思って$n$個のレプリカに対する分配関数を計算していくのですが、この記事では理論的な計算ではなく数値的にレプリカを用意してスピングラスの秩序を見ていきたいと思います。

Sherrington-Kirkpatrick 模型

スピングラスの計算をするときに理論的にも扱いやすいSherrington-Kirkpatrick 模型（SK模型）を題材にします。SK模型は全結合のランダム相互作用のある系で以下のハミルトニアンで記述されます。

$$H =- \sum_{i < j}J_{ij}\sigma_i\sigma_j - h\sum_i\sigma_i$$

ここで相互作用係数$J_{ij}$は平均が$J_0/N$で分散が$J^2/N$のガウス分布から生成されているとします。

$$J_{ij} \sim \mathcal N(J_0/N, J^2/N)$$

この系のマクロな特徴づけを以下の各レプリカ$\alpha, \beta$の間の内積で行うことにします。

$$\sum_i \langle\sigma_i^\alpha\rangle \langle\sigma_i^\beta\rangle$$

$\sigma_i^\alpha \sigma_j^\beta$はそれぞれレプリカ$\alpha$の$i$番目のスピンとレプリカ$\beta$の$j$番目のスピンを表しています。$\langle \cdot \rangle$はボルツマン分布による平均値です。さらにガウス分布から生成している$J_{ij}$についての平均を$[\cdot ]$で記述することにすることにして以下の$q_{\alpha, \beta}$を定義します。

$$q_{\alpha, \beta} = \left[\sum_i \langle\sigma_i^\alpha\rangle \langle\sigma_i^\beta\rangle\right]$$

この$q_{\alpha, \beta}$がスピングラスを特徴づけるスピングラス秩序変数となります。磁化 $m = \sum_i \langle\sigma_i\rangle$が有限の値を持つときスピングラス変数も0ではない量になりますが、逆にスピングラスのようにスピン配列はバラバラでも配列が時間的に固まってしまっている場合は磁化$m$は0になりますが、スピングラス秩序変数は有限の値を持ちます。

この$q_{\alpha, \beta}$はレプリカ法による分配関数の計算の際に積分変数として表れてきて、スピングラスの理論的な解析のキーとなる変数なのですが、$q_{\alpha,\beta}$は人工的に導入したレプリカの番号$\alpha, \beta$に依存していますが、人工的に導入したパラメータなので、依存しなくてもよい気がします。なのでまず理論的に解析する一歩としてはこの$\alpha, \beta$の依存性について無視できること（レプリカ対称性 $q_{\alpha, \beta} = q$）を仮定して計算を進めます。このようにして$\beta, J_0$というマクロなパラメータだけで$m, q$といった秩序変数の解を求めることができます。このような解をレプリカ対称解と呼びますが、レプリカ対称性を仮定したことで、エントロピーが負になるようなパラメータ領域があり、物理的に正しくない結果を導いてしまいます。今回パリージさんがノーベル賞を受賞したのはこの$q_{\alpha, \beta}$に関して、解析的に計算を進めることができつつ、レプリカ対称解の問題を解決する$\alpha, \beta$の依存性、つまりレプリカ対称性の破り方を発明した功績が評価されています。

OpenJijを使ったモンテカルロ計算

ではさっそくスピングラス秩序変数をOpenJijで計算していきましょう。

OpenJijはイジングモデルを使ったモンテカルロシミュレーションを気軽に行うことができるOSSライブラリです。

pip install openjij

で気軽にインストールすることができます。

OpenJijではシミュレーテッドアニーリングを実行する関数がついていますが、その温度スケジュールを任意に操作することができます。そのため温度を一定に設定することで一定温度からのモンテカルロサンプリングを行うことが可能です。その機能を使ってイジングモデルのモンテカルロサンプリングを行ってスピングラス秩序変数の計算に必要な平均値を計算して、各温度で$q$の分布がどうなるかを確認してみましょう。

SK模型の用意

まずSK模型のイジングモデルを生成する関数を用意します。OpenJijではイジングモデルの相互作用係数を$i, j$をキーとして相互作用の大きさを値とした辞書型で受け付けるのでその形でSK模型を生成します。

import numpy as np

def sk_model(N: int) -> Dict[Tuple[int, int], float]:
    J = {}
    for i in range(N):
        for j in range(i+1, N-1):
            J[i, j] = -1*np.random.normal(0, 10/N)
    return J

ここでは縦磁場無しで、平均値0, 分散が$1/N$のガウス分布から相互作用係数を生成します。

OpenJijを使った一定温度でのモンテカルロシミュレーション

ではさっそくOpenJijを使ってイジングモデルのモンテカルロシミュレーションを行います。

OpenJijでは温度スケジュールを逆温度とその温度でのモンテカルロステップ数の組み合わせの配列で指定します。具体的な実装を見ていきましょう。

import openjij as oj

# beta=5で100モンテカルロステップ実行
beta = 5
mcs = 100
beta_schedule = [[beta, mcs]]
# [[beta, mcs], [beta2, mcs2]]とすると途中で温度が切り替わる
# OpenJijのシミュレーテッドアニーリングはこのスケジュール機能で温度を下げていくスケジュールを生成してアニーリングを行っている

num_reads = 10
sampler = oj.SASampler()
# num_reads を指定することで独立にnum_reads個のスピン配列をサンプルしてくる
response = sampler.sample_ising(h={}, J=sk_model(N=10), schedule=beta_schedule, num_reads=num_reads)

print(response.record["sample"])

これで完了です。帰ってきた10個のスピン配列がプリントされます。SKモデルから生成してサンプリングした結果なのでバラバラな結果が出力されると思います。

ここで注意しないといけないのがOpenJijの引数で指定するJとhは以下のハミルトニアンにおけるJ, hです。

$$H = \sum_{i, j}J_{ij}\sigma_i\sigma_j + \sum_i h_i \sigma_i$$

物理の教科書などの定義とはJ, h の符号が逆であることに注意してください。

まずここまで出来たら sk_model の平均値を少し変更したり、betaを変更してみたりして結果が変わることを確かめてみてください。

スピングラス秩序変数の計算

ではさっそくSK模型のスピングラス秩序変数を計算しましょう。

from typing import Dict, Tuple
import numpy as np
import openjij as oj
import matplotlib.pyplot as plt
from openjij.sampler import response

def sk_model(N: int) -> Dict[Tuple[int, int], float]:
    J = {}
    for i in range(N):
        for j in range(i+1, N-1):
            J[i, j] = -1*np.random.normal(0, 10/N)
    return J

num_system = 100    # Jの平均のためのサンプリング回数
N = 50              # スピンの数
num_replica = 50   # レプリカの数
mcs = 500           # モンテカルロステップ数
beta_list = 1/np.array([0.1, 0.5, 1, 2, 5])  # 掃引する逆温度
q_list = []
for beta in beta_list:
    q_beta = []
    qab_beta = []
    for _ in range(num_system):  # Jに関する平均のためのループ
        # スピングラスは最初から低温にしてしまうと全然動かないことがあるので少しアニールする
        beta_schedule = [[beta/10, 100], [beta/2, 100], [beta, mcs]]
        # SK模型の相互作用係数
        Jsk = sk_model(N = N)
        # OpenJijを使ったサンプリング
        sampler = oj.SASampler()
        response = sampler.sample_ising(h={}, J=Jsk, schedule=beta_schedule, num_reads=num_replica)
        # num_reaplica個のレプリカの結果を取得 samples.shape == (num_replica, N)
        samples: np.ndarray = response.record["sample"]
        # レプリカ間の全組み合わせの内積を計算
        q = 1/N*np.einsum("ai, bi -> ab", samples, samples)
        q_beta.append(q)
    q_list.append(np.array(q_beta).flatten())

beta_index = 0
plt.title("T = {}".format(1/beta_list[beta_index]))
plt.hist(q_list[beta_index], bins=20)
plt.show()

beta_index を変えることで各温度でのスピングラス秩序変数$q$の分布をみることができます。

例えば温度が低いところだと

f:id:jijtech:20211203203505p:plain — T=0.1でのサンプリングしたqのヒストグラム

のように両端にピークがあってそれ以外は平らになっている（スピングラスは多くの多谷構造を持つのでスピングラス相ではスピングラス秩序パラメータの分布がフラットに近くなる）ことが見えます。

少しずつ温度を上げていくと f:id:jijtech:20211203203655p:plain

f:id:jijtech:20211203203715p:plain

f:id:jijtech:20211203203743p:plain

このように分布が q=0に近くなっていきます。これは温度揺らぎによってスピンがバラバラに動き、レプリカ間のオーバーラップがほぼなくなっていく様子を表しています。

今回は古典スピン系のスピングラスの数値計算をOpenJijで行う方法を紹介しました。OpenJijでは横磁場を入れた量子イジングモデルの量子モンテカルロシミュレーションも行うことが可能なので今度は量子系のシミュレーション方法についても解説したいと思います。

参考文献

スピングラスに関する理論を知りたい方は以下の書籍がおすすめです。

西森 - スピングラス理論と情報統計力学

2021-06-25

熱浴法、模擬焼き鈍し法と種々のアルゴリズム、及びそれらの収束定理

熱浴法で更新する場合の、イジング模型と呼ばれる $\{-1,+1\}^N$ を定義域とする多項式函数の最小に近い解を求めるための、有名な古典アルゴリズムである模擬焼き鈍し法の説明およびその収束定理（Geman-Geman）を、原論文を和訳し、さらに証明もコンパクトなものに纏めました。

さらに、一般の測度空間の可算直積を基礎の空間とする、ある種の自然な条件を充たす離散時間マルコフ過程についての収束性についての一般定理も、私が証明を構成したものを纏めました。（調べた限り、これは新しい結果なのではと思います。）

基礎知識は、位相、測度、及び $L^p$ 空間の完備性程度を理解しておけば理解できると思います。確率論は、用語の定義程度で充分です。

ブログに直接書くと長すぎて数ページに及んでしまうので、以下に一般公開されたグーグルドライブのファイルへのリンクを貼るという形式にさせて頂きました。

それでは、興味のある方は是非ご覧になって下さい。

https://drive.google.com/file/d/1BNn9ZX3yELT7hA2NhE4sFudGG26UfxxL/view?usp=sharing

2021-06-14

単位円上の有理点のなす群と、p=a²+b²の形に表せる素数

代数幾何

単位円上の有理点のなす群と、 $p=a^2+b^2$ の形に表せる素数

本記事では、単位円上の有理点のなす群 $\rm SO(\mathbb Q)$ を考察します。単位円上の有理点は、「ピタゴラスの三つ組」、つまり

a^2+b^2=c^2

を充たす $(a,b,c) \in \mathbb Z^3$ の全体の、定数倍を同一視した商集合に一対一対応しています。

まず、結構知られていることだとは思いますが、単位円上の任意の有理点がパラメータを用いて表示できることから復習していきましょう。

点 $(-1,0)$ を通る、傾き $t$ の直線は

y=t(x+1)

と表せます。単位円

x^2+y^2=1

との交点を求めるためには、二つの式を連立させて解けば良いです。つまり

t^2(x+1)^2+x^2=1

整理して

(t^2+1)x^2+2t^2x+t^2-1~0

だから、結局、交点は

(x,y)=\left(\frac{1-t^2}{1+t^2},\frac{2t}{1+t^2}\right)

と表示できます。これらに $(-1,0)$ を付け加えれば、単位円上の有理点の全体は

\left\{ \left(\frac{1-t^2}{1+t^2},\frac{2t}{1+t^2}\right)|t \in \mathbb Q \right\} \cup\{(-1,0)\}

のような形になっています。

さて、当然ながら、単位円上の点には「角度」が対応していて、単位円上の点の回転に対応する「群」を考えることができます。つまり、演算を

(s,t)\cdot(u,v) \mapsto(su-tv,sv+tu)

により定義するということです。有理点に限らない、全ての単位円上の点のなす群は $\rm SO(2)$ と呼ばれるものですが、その中で有理点のなす部分群を考察してみましょう。

$t=a/b$ と既約分数で表示すれば

\left(\frac{1-t^2}{1+t^2},\frac{2t}{1+t^2}\right) =\left(\frac{a^2-b^2}{a^2+b^2},\frac{2ab}{a^2+b^2}\right)

の形になります。言い換えれば、ピタゴラスの三つ組はいつも

(a^2-b^2,2ab,a^2+b^2)

の形をしているという訳です。虚数単位 $\sqrt {-1}$ を使って、対応する複素数

\frac{a^2-b^2}{a^2+b^2}+\frac{2ab}{a^2+b^2}\sqrt{-1}

を考えます。単位円上の点のなす群を「良く知る」ためには、どんな群と同型になっているか、言い換えれば群の「生成系」を知りたくなってきます。（有理点は可算個ですから、高々可算個の群の直和または直積で表されるでしょう。）その為には、「どのような条件の下で $\rm SO(\mathbb Q)$ の元は分解（二つの元の積で表されるか）できるか？」を考えるのが重要です。

そこで、ガウス整数環

\mathbb Z[\sqrt{-1}]:=\{m+n\sqrt{-1}|(m,n) \in \mathbb Z ^2\}

を考えます。平面の上で、原点 $O=(0,0)$ と $(m,n)$ を結ぶ直線と $y$ 軸に平行な直線 $x=1$ は点 $(1,n/m)$ で交わりますが、上で考えた

\frac{n}{m} \mapsto \left(\frac{n^2-m^2}{n^2+m^2},\frac{2nm}{n^2+m^2}\right)

なる写像と合わせれば、ガウス整数環の乗法群 $\mathbb Z[\sqrt{-1}]^\times(:=\mathbb Z[\sqrt{-1}]/\{0\})$ から ${\rm SO(\mathbb Q)}$ への写像 $f$

m+n\sqrt{-1} \mapsto \left(\frac{n^2-m^2}{n^2+m^2},\frac{2nm}{n^2+m^2}\right)

を得ます。これが全射になっていることは定義から明らかです。さらに

f(1+0\sqrt{-1})=(1,0)

が成り立っていて、 $f$ は対応する角度が２倍のものに写し、 $\mathbb Z[\sqrt{-1}]^\times$ 及び $\rm SO(\mathbb Q)$ のそれぞれの群演算も回転を意味しているので

f((m+n\sqrt{-1})(m^\prime+n^\prime\sqrt{-1}))=f(m+n\sqrt{-1})\cdot f(m^\prime+n^\prime\sqrt{-1})

が成り立つことは直感的にも明らかでしょう。つまり、 $f$ はアーベル群の全射準同型になっているということです。

ところで、ガウス整数環 $\mathbb Z[\sqrt{-1}]$ は代数学の分野で良く調べられている対象です。 $\mathbb Z[\sqrt{-1}]$ の既約元に関する基本的事実を証明なしで引用しておきます。例えば「https://pomb.org/pdf/2015algebra.pdf」が参考になるでしょう。

定理： $\mathbb Z[\sqrt{-1}]$ の既約元は、以下のもののみからなる。

$(1)$ 有理整数（ $\mathbb Z$ に於ける素数）に対し、 $p=2$ または $p\equiv 1(\mod 4)$ ならば、 $p=m^2+n^2$ となるような $m,n \in \mathbb N$ が一意的に存在し、 $m+n\sqrt{-1}$ は既約である。

$(2)$ 有理整数（ $\mathbb Z$ に於ける素数）に対し、 $p\equiv 3(\mod 4)$ ならば $p$ は $\mathbb Z[\sqrt{-1}]$ に於いても既約である。

例えば $(1)$ に関しては

$2=1^2+1^2$

$5=1^2+2^2$

$13=2^2+3^2,\cdots$

となっていて、 $(2)$ の例を挙げるとすれば

$3,7,11,19,\cdots$ となっています。

この結果を、単位円上の有理点のなす群 ${\rm SO}(\mathbb Q)$ に適用してみましょう。

$\mathbb Z[\sqrt{-1}]$ は一意分解整域なので、素因数分解の一意性が成り立つことにより、 $\mathbb Z[\sqrt{-1}]^\times$ の任意の元は上のような既約元たちの積で書けます。準同型 $f$ は全射準同型なので、やはり $\rm SO(\mathbb Q)$ の元はそれら既約元たちの像の積で書けます。（一意性は、この時点では保証されません。）言い換えれば、 $\mathbb Z[\sqrt{-1}]^\times$ の既約元の像が $\rm SO(\mathbb Q)$ の生成系になっているということです。

$(1)$ の形の既約元の $\rm SO(\mathbb Q)$ に於ける像は

f(m+n\sqrt{-1})=\left(\frac{n^2-m^2}{p},\frac{2nm}{p}\right)

です。

$(2)$ の形の既約元の像は $f(p)=(1,0)$ となっています。

さて、これらより、 ${\rm SO}(\mathbb Q)$ の生成系として

\left\{\left(\frac{n^2-m^2}{p},\frac{2nm}{p}\right)|p=2,p\equiv 1(\mod 4)\right\}

がとれることが分かりました。ここで、式に現れる $p$ は素数を動き、 $n>m>0$ とします。

例えば $p=2$ の場合、 $2=(1+\sqrt{-1})(1-\sqrt{-1})$ より $n=1,m=1$ ととれて

\left(\frac{n^2-m^2}{p},\frac{2nm}{p}\right)=(1,0)

のようになっており、これは $\rm SO(\mathbb Q)$ の位数 $4$ の部分巡回群を生成し、その中には $\rm SO(\mathbb Q)$ の単位元 $(1,0)$ も含まれます。

これら生成系の間に関係はあるのでしょうか？つまり、もし自明な関係しかないとすれば、準同型定理により、 $\rm SO(\mathbb Q)$ は生成元のそれぞれが生成する部分群の直和になっていることが言えます。これに関しては、実際、自明な関係しかないことが言えます。

補題：

\left\{\left(\frac{n^2-m^2}{p},\frac{2nm}{p}\right)|p=2,p\equiv 1(\mod 4)\right\}

は非自明な関係を持たない。

証明：

ある $a_p \in \mathbb Z$ があって

\prod_{p \equiv 1(\mod 4)} a_p \left(\frac{n^2-m^2}{p},\frac{2nm}{p}\right)=(1,0)

を充たすとすれば、 $\mathbb Z[\sqrt{-1}]$ に引き戻して考えると

\prod_{p \equiv 1(\mod 4)} a_p f(n+m\sqrt{-1})=\prod_{p \equiv 1(\mod 4)} f((n+m\sqrt{-1})^{a_p})=(1,0)

となる。即ち、ある $k,k^\prime \in \mathbb Z$ があって

k\prod_{p \equiv 1(\mod 4)} (n+m\sqrt{-1})^{a_p}=k^\prime

でなければならない。故に、素元分解の存在と一意性により、 $\forall a_p,a_p=0$ である。

$\square$

上の証明から、 $p=2$ でないならば

\left(\frac{n^2-m^2}{p},\frac{2nm}{p}\right)

の形の元は何倍しても $(1,0)$ になることがない、つまりこれにより生成される部分巡回群の位数は無限であることも分かりました。以上により、 $p=2$ の場合の

\left(0,1\right)

により生成される部分群 $G_2 \simeq \mathbb Z/4\mathbb Z$ 、及び $p\equiv 1(\mod 4)$ なる素数の場合の

\left(\frac{n^2-m^2}{p},\frac{2nm}{p}\right)

により生成される部分群 $G \simeq \mathbb Z$ の直和、即ち

{\rm SO(\mathbb Q)} \simeq \bigoplus_{p=2,p\equiv 1(\mod 4)}G_p

となることが分かりました。

これは、単位円上の有理点のなす群 $\rm SO(\mathbb Q)$ と $p=a^2+b^2$ と平方数の和で表される素数についての、とても面白い関係と言えるのではないでしょうか？

ついでにですが、実は、ガウスの整数環 $\mathbb Z[i]$ は体 $\mathbb Q(i)$ の中で $\mathbb Z$ 係数多項式の根になるようなものを全て集めたものになっています。一般に、体 $K$ が有理数体 $\mathbb Q$ を含み、 $\mathbb Q$ 上のベクトル空間として有限次元の場合（これを、 $K$ は $\mathbb Q$ 有限次拡大であると言います）、 $K$ の元で、ある $\mathbb Z$ 係数多項式の根になっているようなものの全体（これを、 $\mathbb Z$ の $K$ に於ける整閉包と言います。）は「代数的整数環」と呼ばれるもので、デデキント整域というある特別な環になっていて、デデキント整域に於いては普通の整数環 $\mathbb Z$ における素因数分解のある種の一般化である「イデアルの素イデアル分解」ができるという性質を持っています。ここではこれ以上詳しいことはあまり述べませんが、もし興味を持った方がおられれば、「代数的整数論」方面の文献を見るといいかも知れません。

また、このように、ある体の上で定義された、多項式により表される図形の点の間に群構造が入るようなものは、代数幾何や数論幾何でよく研究されている対象でもあります。有名どころとしては楕円曲線があり、その一般化として群多様体やヤコビ多様体、アーベル多様体と言ったものがあります。個人的には、純粋に代数的な議論をやるよりは、この記事のような幾何的直感の働く議論の方が好きなのですが、皆さんはどうでしょうか？もし興味を持たれたら、宣伝と言う訳ではありませんが、永井保成先生の「代数幾何学入門：代数学の基礎を出発点として」はお薦めの本です。

参考：

https://pomb.org/pdf/2015algebra.pdf

https://www.maa.org/sites/default/files/pdf/upload_library/22/Allendoerfer/1997/0025570x.di021195.02p0087x.pdf

雪江明彦「代数学２　環と体とガロア理論」

2021-04-10

Nimのインストール

NimはPythonのような文法を持つ静的型付け言語です。実際使ってみるとPythonの皮を被ったCって感じです。

まだ日本語の記事が少ない（2021/4）ので、こちらでインストール方法を紹介します。

choosenim のインストール

Nimをインストール標準的な方法はNimのバージョン管理を行う choosenin を使います。Pythonのpyenvのようなツールです。

choosenim: https://github.com/dom96/choosenim

choosenimのREADMEにある手順でインストールしていきます。Mac, Linux では以下のようにcurlで直接インストールします。

curl https://nim-lang.org/choosenim/init.sh -sSf | sh

インストールした後は~/.nimble/binにパスを通せと言われるので

export PATH="/$HOME/.nimble/bin/:$PATH"

を~/.zshrcまたは~/.bashrcに追記してください。

source ~/.zshrc

(bashを使っているなら source ~/.bashrc) で変更を反映させると choosenim コマンドが使えるはずです。

choosenim -V

でバージョンが出力されれば成功です。

そしてchoosenimを使ってNimの安定バージョンをインストールします。

choosenim stable

これでNimの安定バージョンがインストールできるはずです。

nim -V

と

nimble -V

の両方がうまく動けばNimがインストールできたことになります。

2021-03-27

セグレ埋め込みと射影多様体の積（２）～代数幾何学～

代数幾何

セグレ埋め込みと射影多様体の積（１）～代数幾何学～

https://jijtech.hatenablog.com/entry/2020/11/13/105951

から続きます。

射影多様体の積

いよいよ、射影代数多様体の積を定義します。任意の準射影多様体 $X \sub \mathbb P^n$ 及び $Y\sub \mathbb P^m$ についてセグレ埋め込みによって直積集合としての $X \times Y$ に位相構造を入れたもの $X \times Y \sub \mathbb P^n \times \mathbb P^m$ が、射影多様体の積です。まず、 $\mathbb P^n,\mathbb P^m$ の閉集合 $Z,W$ に対し $Z \times W$ は閉集合であることを示しましょう。

証明：

Z\times W=\cap_i Z(f_i)\times\cap_j Z(g_j)=\cap_{i,j}Z(f_i) \times Z(g_j)=\cap_{i,j}(Z(f_i)\times \mathbb P^m \cap \mathbb P^n\times Z(g_i))

が成り立つので、 $Z(f) \times \mathbb P^m$ という形の集合が閉集合であることを示せば充分である。しかし、 $f(x_0,...,x_n)=0$ と $\forall j,f(x_0y_j,...,x_ny_j)=0$ が同値であることにさえ注意すれば明らかに、 $\mathbb P^N$ の閉集合 $\square$

また、 $X,Y$ が局所閉なら $X\times Y$ が局所閉であることも容易に確認できます。

証明：

任意の準射影多様体は局所閉だから、閉集合 $T_1,T_2,S_1,S_2$ を用いて

X \times Y=(T_1-T_2)\times (S_1-S_2)=T_1\times S_1-T_2\times S_1\cup T_1 \times S_2

と書ける。

$\square$

次に、局所閉集合に関する節で述べたように、局所閉集合の閉包が既約ならば、これは準射影多様体であることが言えます。よって、閉包を考えることで $X,Y$ が既約の時に $X\times Y$ が既約であることを言えば良いことが分かります。

命題：

位相空間 $X$ は、既約なら連結である。

証明：

背理法により示す。 $X$ が連結でないと仮定すると、 $X$ に真に含まれるような開かつ閉なる集合 $A$ が存在する。 $X-A$ は空でなく、開かつ閉であって $X=A\cup(X-A)$ となる。

命題：

局所既約（任意の $x \in X$ に対し、既約な開近傍がとれること。既約集合の開部分集合は既約だから、これは $x$ の各近傍が既約な近傍を包むことと同値である）な位相空間 $X$ に於いて、既約成分は開である。

証明：

既約成分 $\cal I$ に対し $x \in \cal I$ を任意にとると、 $x$ の既約な開近傍 $U_x$ が $X$ の中でとれる。 $U_x \cap \cal I$ は空でなく、 $U_x$ の中では既約成分である。これは以下のようにして分かる。 $U_x \cap {\cal I}$ を真に含む既約集合が $U_x$ の中にあるとする。 $U_x \sub \cal {I}$ ならば何も言うことはないが、そうでないならば閉集合 ${\cal I}-U_x$ は空でない。その閉包 $V$ は ${\cal I}\not \subset V$ で ${\cal I}=({\cal I}-U_x)\cup V$ 、つまり $\cal I$ が既約でないことになる。これは矛盾。従って、 $U_x\cap {\cal I}= U_x$ 即ち $U_x \sub \cal I$ を得る。これは $\cal I$ が開であることを示す。

$\square$

命題：

局所既約な位相空間 $X$ に於いて、連結成分であることと既約成分であることは同値である。

証明：

$X$ の連結成分 $\cal I$ は、全ての $\cal I$ と交わる開かつ閉である集合に含まれる。よって特に $\cal I$ と交わる既約成分をとると、 $X$ が局所既約だからそれは開かつ閉であることにより、それに含まれる（適当に $x \in {\cal I}$ をとり、それを含む既約成分をとれば良い。その存在は既に示している）。既約なら連結であるから、 $\cal I$ の極大性により $\cal I$ は既約でなければならない。故に、連結成分は既約成分である。逆はこれより明らか。

$\square$

命題：

位相空間 $X$ に於いて、既約な開被覆 $\{U_\lambda\}_{\lambda \in \Lambda}$ があってそれらの共通部分が空でないならば、 $X$ は既約である。

証明：

既約なら連結であるから、連結な集合族 $\{M_\lambda\}_{\lambda \in \Lambda}$ の共通部分 $\cap_\lambda M_\lambda$ が空でないなら合併 $\cup_\lambda M_\lambda$ も連結であることにより、 $X=\cup_\lambda U_\lambda$ は連結である。また、 $X$ は局所既約でもあるから前命題により既約である。

$\square$

では最後に、既約集合 $X \sub \mathbb P^n,Y \sub {\mathbb P}^m$ のセグレ埋め込み $X \times Y$ が既約であることを示しましょう。

証明：

$H_i:x_i=0,H^\prime_j:y_i=0$ とし、一般性を失わずに、 $\forall i,X \not \sub H_i$ 及び $\forall j, Y \not \sub H^\prime_j$ を仮定する。セグレ埋め込みの中で $(\mathbb P^n-H_i)\times (\mathbb P^m-H^\prime_j)$ は開集合だが、これは $\mathbb A^{n+m}$ と同型。従って $U_{ij}:=(X-H_i)\times (Y-H^\prime_j)$ は、アファイン多様体の時の議論により既約。つまり、 $X \times Y$ のアファイン開被覆がとれている。さらに、共通部分 $U:=\cap U_{ij}=(X-\cup _iH_i)\times (Y-\cup_jH^\prime_j)$ は空でないので $X \times Y$ は既約である。

$\square$

参考：

[1]http://home.p07.itscom.net/strmdrf/set07.htm

[2]ハーツホーン-代数幾何学１

[3]https://arxiv.org/pdf/1411.5901.pdf

2021-02-26

航空機荷物搭載最適化問題

D-Wave 数理最適化量子アニーリング

この記事の概要

Pilon et al., 2021, "Aircraft Loading Optimization - QUBO models under multiple constraints"を読み、その理解を深めるために内容をまとめたものです。

この記事の概要
Aircraft Loading Optimization
モデル化
不等式制約の計算法
最適化計算の実装
数値実験結果
- 古典ソルバー
- 量子ソルバー
結言
参考文献
文責

Aircraft Loading Optimization

航空機に搭載できる荷物には限りがあります。荷物コンテナをどのように航空機内に配置すれば、安全かつ運べる荷物重量が最大となる輸送が行えるかを最適化するのがこの最適化問題の目的です。

モデル化

モデル化にあたり以下のようなバイナリ変数を用います。

$\displaystyle{ p_{i, j} = \left\{ \begin{array}{ll} 1 & (コンテナiを位置jに配置するとき) \\ 0 & (それ以外) \end{array} \right. \tag{1} }$

添字 $i=1, \dots, n$ は利用可能コンテナの番号、 $j=1, \dots, N$ は荷物を置くことができる位置を表します。

目的関数

コンテナ $i$ は重量 $m_i(>0)$ を持ちます。またコンテナサイズとしてT1(中くらいのサイズで1の広さを占有), T2(小さいサイズで1/2の広さを占有), T3(大きいサイズで2の広さを占有)の3種類を考えます。今、最大化したいのは運ぶ荷物の総重量としましょう。すると、目的関数はこれらを用いて

$\displaystyle{ f^\mathrm{obj}(\mathbf{p}) = - \sum_{i=1}^{n} \sum_{j=1}^{N} t_i m_i p_{i, j} \tag{2} }$

のように表現することができます。ここで

$\displaystyle{ t_i = \left\{ \begin{array}{ll} 1 & (\mathrm{T1}, \mathrm{T2}のとき) \\ \frac{1}{2} & (\mathrm{T3}のとき) \end{array} \right. \tag{3} }$

です。例えばあるT3コンテナの重さが2Mでこれを1つ搭載することは、重さMのT1コンテナを2つ搭載することと等価です。よってT3をこのように区別しています。

制約

Payload Limits (荷重制限)

No-overlaps (コンテナにおける位置に対する制約)

1つの場所におけるコンテナは中/大サイズのコンテナを1つ、または小サイズのコンテナなら2つまでです。この制約を数理モデルで表現すると

$\displaystyle{ \sum_{i=1}^n d_i p_{i, j} \leq 1 \tag{4} }$

ここで

$\displaystyle{ d_i = \left\{ \begin{array}{ll} 1 & (\mathrm{T1, T3}のとき) \\ \frac{1}{2} & (\mathrm{T2}のとき) \end{array} \right. \tag{5} }$

は位置の占有を表現するための係数です。この制約をQUBOで定式化すると以下のようになります。

$\displaystyle{ f^{\bar{O}}_j = P^{\bar{O}} \left(\sum_{i=1}^n d_i p_{i, j} + \sum_k c_k v_k -1 \right)^2 \tag{6} }$

ここで $v_k$ はスラック変数、そして $c_k$ はその係数です。もし $j$ という同じ位置に不適切なコンテナ量が割り当てられたら、これが正の値となります。

No-duplication (1つのコンテナは1つの位置にしか置けない制約)

小/中サイズのコンテナは1つの位置にしか置けません。大サイズのコンテナは隣り合う2つの位置に置けますが、これもある1つの位置に限られます。これを数理モデルで表現すると

$\displaystyle{ t_i \sum_{j=1}^N p_{i, j} \leq 1 \tag{7} }$

となります。よってこの制約をQUBOで定式化すると以下のようになります。

$\displaystyle{ f_i^{\bar{D}} = P^{\bar{D}} \left(t_i \sum_{j=1}^N p_{i, j} + \sum_k c_k v_k -1 \right)^2 \tag{8} }$

Contiguity for big containers (大サイズコンテナに対する制約)

(7)式において、大サイズのコンテナに関しては連続する位置を占有させる必要があります。よって(7)式を拡張する必要があります。もしコンテナ $i$ が大サイズであるとすると

$\displaystyle{ \sum_{j=1}^{N-1} p_{i,j } p_{i, j+1} = \frac{1}{2}\sum_{j=1}^N p_{i, j} \tag{9} }$

が成り立ちます。このコンテナを置くことが決まっているのであれば右辺は1、どこにも置かない場合には0です。よってこれをQUBOで定式化すると以下のようになります。

$\displaystyle{ f_i^{C} = P^{C} \left( \frac{1}{2} \sum_{j=1}^N p_{i, j} - \sum_{j=1}^{N-1} p_{i, j} p_{i, j+1} \right) \tag{10} }$

この制約式の形はすでに2次式になっているので、これ以上の変形の仕様がありません。しかし $P^C$ と $P^{\bar{D}}$ の値を適切に取ることで、 $f^C_i + f^{\bar{D}}_i$ の値を負にならないようにすることができます。 $\sum_{j=1}^N p_{i, j} = 1$ のときは $f^C_i > 0$ となるのは明かなので、ここからは $\sum_{j=1}^N p_{i, j} >2$ のときに $f_i^C>0$ となる場合を考えましょう。

$\sum_{j=1}^N p_{i, j} = m$ のとき、 $f_i^C$ が最も小さくなる最悪の場合として考えられるのは $\sum_{j=1}^{N-1} p_{i, j} p_{i, j+1} = m-1$ です。また $f_i^{\bar{D}}$ が最も小さくなるときを考えると

$\displaystyle{ f_i^C + f_i^{\bar{D}} = P^C \left( \frac{1}{2} m - (m-1) \right) + P^{\bar{D}} \left(\frac{1}{2} m-1 \right)^2 = \frac{1}{2} P^C (2-m) + \frac{1}{4} P^{\bar{D}} (2-m)^2 }$

これが正の値となるには

$\displaystyle{ P^{\bar{D}} > \frac{2(m-2)}{(m-2)^2} P^{C} }$

であれば良いことがわかります。 $\frac{2(m-2)}{(m-2)^2} \ (m>2)$ は $m=3$ のとき最大値2を取ることがわかるので、

$\displaystyle{ P^{\bar{D}} > 2 P^C \tag{11} }$

であれば常に $f^C_i + f^{\bar{D}}_i > 0$ となります。

最大積載量以下でなければならない制約

飛行機に搭載できる荷物の最大総重量を $W_p$ と書くと、この数理モデルは

$\displaystyle{ \sum_{i=1}^n \sum_{j=1}^N t_i m_i p_{i, j} \leq W_p \tag{12} }$

となります。これをQUBOで定式化すると

$\displaystyle{ f^W = P^W \left( \sum_{i=1}^n \sum_{j=1}^N t_i m_i p_{i, j} + \sum_k c_k v_k -W_p\right)^2 \tag{13} }$

です。

Center of Gravity Limits (重心制約)

荷物コンテナの配置が悪いと、航空機は姿勢制御に燃料を費やし最適な輸送となりません。よってここでは搭載した荷物の重心を計算し、それに対する制約を課します。ここではコンテナを配置する位置を

$\displaystyle{ x_j = \frac{L}{N} \left( j - \frac{N}{2} \right) - \frac{L}{2N} \tag{14} }$

ターゲットとする重心位置との一致

重心がこの位置に一致していることが望ましいとする値を $x_{cg}^t$ とすると

$\displaystyle{ \hat{x}_{cg} = \frac{\bar{x}_{cg}}{\underline{x}_{cg}} = \frac{\sum_{i=1}^n \sum_{j=1}^N t_i m_i p_{i, j} x_j + W_e x_{cg}^e}{\sum_{i=1}^n \sum_{j=1}^N t_i m_i p_{i, j} + W_e x_{cg}^e} = x_{cg}^t \ \Longrightarrow \ \bar{x}_{cg} = x_{cg}^t \underline{x}_{cg} }$

ここで $W_e$ は荷物を搭載していないときの航空機の重量、そして $x_{cg}^e$ は航空機の重心です。これをQUBOで定式化すると

$\displaystyle{ f^{C_t} = P^{C_t} \left( \bar{x}_{cg} - x_{cg}^t \underline{x}_{cg} \right)^2 \tag{15} }$

となります。

上限と下限の制約

先ほどのターゲットにピタリと一致していなくとも、荷物を搭載したときの重心が $[x_{cg}^{min}, x_{cg}^{max}$ ]の範囲に収まってほしいというのを制約として課します。下限に対する制約は

$\displaystyle{ f^{C_l} = P^{C_l} \left( \bar{x}_{cg} - x_{cg}^{min} \underline{x}_{cg} - \sum_{k} c_k v_k \right)^2 \tag{16} }$

そして上限に対する制約は

$\displaystyle{ f^{C_u} = P^{C_u} \left( \bar{x}_{cg} - x_{cg}^{max} \underline{x}_{cg} - \sum_{k} c_k v_k \right)^2 \tag{17} }$

です。 $P^{C_l}, P^{C_u}$ は $P^{C_t}$ よりも大きくなければいけません。(16), (17)式は(15)式に比べて強い制約でなければならないからです。この論文では $P^{C_l} = P^{C_u} \simeq 10 P^{C_t}$ としています。

Shear Limits (シアに対する制約)

重心位置が同じ2通りの置き方があるとして、一方は重心の周りに重たいものを配置する置き方、もう一方は重心から離れた位置に重たいものを配置する置き方だとします。後者の方法では航空機体に大きなシアが発生し、安全な航行ではなくなる可能性があります。よってここではシアの大きさに対する制約を考えましょう。

線形・対称な場合を考えるとシアの最大値は

$\displaystyle{ S^{max}(x) = \left\{ \begin{array}{ll} S_0^{max} \frac{L + 2x}{L} & (x<0) \\ S_0^{max} \frac{L - 2x}{L} & (x>0) \end{array} \right. \tag{18} }$

と書かれます。ここで $S_0^{max}$ は適当な定数、 $L > 0$ は航空機の長さです。そして $u = 1, \dots, N$ における位置は

$\displaystyle{ x_u = \frac{L}{N} \left( u- \frac{N}{2} \right) \tag{19} }$

そして $u$ の位置に働く、左側からのシアは

$\displaystyle{ \hat{S}^l (u) = \sum_{i=1}^n \sum_{j=1}^u t_i m_i p_{i, j} \quad (u: x_u < 0) \tag{20} }$

同様に $u$ の位置に働く、右側からのシアは

$\displaystyle{ \hat{S}^r (u) = \sum_{i=1}^n \sum_{j=u+1}^N t_i m_i p_{i, j} \quad (u: x_u > 0) \tag{21} }$

わかりやすくするため、ここでは $N$ が偶数の場合を考えます。

左側のシア

$\displaystyle{ \hat{S}^l (u) \leq S^{max} (x_u) \quad (u = 1, \dots, N/2) \tag{22} }$

これをQUBOで定式化すると

$\displaystyle{ f_{x_u}^{S_l} = P^{S_l} \left( \hat{S}^l (u)+ \sum_k c_k v_k - S^{max} (x_u) \right)^2 \tag{23} }$

右側のシア

$\displaystyle{ \hat{S}^r (u) \leq S^{max} (x_u) \quad (u = N/2, \dots, N-1) \tag{22} }$

これをQUBOで定式化すると

$\displaystyle{ f_{x_u}^{S_r} = P^{S_r} \left( \hat{S}^r (u)+ \sum_k c_k v_k - S^{max} (x_u) \right)^2 \tag{23} }$

Total QUBO

これまでの制約と目的関数をまとめると、この最適化問題を解くための全ハミルトニアンは以下のようになります。

$\displaystyle{ H^{tot} = f^{\mathrm{obj}} + f^W + \sum_{i=1}^n f_i^{\bar{D}} + \sum_{i:iがT3の場合}f_i^{C} \sum_{j=1}^N f_j^{\bar{O}} \\ + f^{C_t} + f^{C_l} + f^{C_u} \\ + \sum_{u=1}^{N/2} f_{x_u}^{S_l} + \sum_{u=N/2}^{N-1} f_{x_u}^{S_r} \tag{24} }$

不等式制約の計算法

例えば $(\mathrm{left \ side}) \leq u$ のような不等式制約では、Log encodingスラック変数を用いて

$\displaystyle{ (\mathrm{left \ side}) + \sum_{k=0}^{r_u} 2^k v_k = u }$

のように等式制約に直してから最適化計算を行います。同様に $(\mathrm{left \ side}) \geq l$ のような不等式制約では

$\displaystyle{ (\mathrm{left \ side}) - \sum_{k=0}^{r_l} 2^k v_k = l }$

のようにして等式制約にします。

最適化計算の実装

この論文では以下の3つの手法を用いて最適化計算を実装し、数値実験を行っています。

QBsolv: D-Waveが公開している、Tabu-Searchに基づいたメタヒューリスティック最適化アルゴリズム。古典ソルバーとして使用。
D-Wave 2000Q: D-Waveが開発した量子アニーリングハードウェアD-Wave 2000Q_5のキメラグラフに問題を埋め込み、最適化計算を行う手法。量子ソルバーとして使用。
Exact solver: 全ての状態を計算しその中から一番エネルギーの低い状態を探すという、力まかせ(brute-force)ソルバーによる計算手法。

数値実験結果

論文では3つの場合に対してExact solver, QBSolv, D-Wave 2000Q手法を用いて計算を行い、それぞれのベンチマークを取得しました。

Case: PL (Only Payload Limits)
Case: PL + CL (Payload + Center of Gravity Limits)
Case: PL + CL + SL (Payload + Center of Gravity + Shear Limits)

それではベンチマーク結果を見ていきましょう。

古典ソルバー

Average Time-to-Solution

QBsolvでこの最適化問題を解いたときのAverage Time-to-Solution (TTS)を示します。制約が増えるに連れてスラック変数の数が増えるため、計算時間が増えていると解釈できます。

f:id:jijinc_nakamura:20210225041846p:plain — 500回の試行における平均TTS.

成功確率

以下は解の成功確率を示したものです。CaseがPL+CLの場合、PLが満たされている解が出た確率は98.0%、CLが満たされている解が出た確率は100%というように結果を見ます。

f:id:jijinc_nakamura:20210225042121p:plain — モデル別、制約を満たしている成功確率。

下図はPLとPL+CLのときの、重心の出現分布です。PLのみ(左パネル)では重心の位置が $x_{cg}^{min}$ (黄色線)よりも小さな値になっている解が出現していますが、PL+CL(右パネル)ではそれが解消されています。

f:id:jijinc_nakamura:20210225042329p:plain — 重心位置の分布

シア制約

以下の表はシアのエラー割合です。0 errorsは全ての位置がシア制約を満たした解が出現した割合、1 errorsはシア制約を満たせなかった位置が1つあるような解が出現した割合などです。

f:id:jijinc_nakamura:20210225042500p:plain — シア制約に終えるエラー割合

シア制約(SL)を導入することでシアの大きさが規定範囲内に収まっているような解が出現する割合が若干増えているように見えます。しかし、PL+CLとPL+CL+SLのerrors=1列を見ると、PL+CL+SLの方が値が大きくなっています。これは重心制約とシア制約が拮抗しているためと考えられますが、さらなる調査が必要としています。

量子ソルバー

以下の表は厳密計算(exact), 古典ソルバー(qbsolv), 量子ソルバー(dwave)でのベンチマーク比較です。

f:id:jijinc_nakamura:20210225042639p:plain — 3つのソルバーにおける、平均TTS, 実行可能解の出現率, 最適解の出現率の比較

qbsolvは100%の確率で実行可能解を算出し、そのうち33%は最適解を導いています。dwaveは80%で実行可能解を算出し、最適解はわずか4%です。これはD-Wave 2000Qにノイズが多く存在するためと考えられます。

さらに下の表は、得られた解から算出された荷物の総重量の最大値と平均値です。qbsolvに比べてdwaveが少ない重量の解が現れていることがわかります。

f:id:jijinc_nakamura:20210225042829p:plain — 複数回試行したときの最大重量と平均重量の比較

結言

航空機に積む荷物の重量や重心・シアに関しての制約の定式化の方法と、それら制約を満たしながらも積載量を最大化する問題を解く方法をご紹介いたしました。
またこの問題において最適解を見つける速度はD-Wave量子アニーラー上で量子ソルバーを用いたときの方が高速ですが、安定して良い解を出せるのは古典ソルバーであることもわかりました。

参考文献

文責

中村翔、株式会社 Jij
Sho K. NAKAMURA, Jij Inc.
憂いの篩 -Pensieve-

※ 疑問点や間違いなどございましたら、お気軽にブログ上でご質問ください。

イジングモデルと相転移

スピングラス

レプリカ法

Sherrington-Kirkpatrick 模型

OpenJijを使ったモンテカルロ計算

SK模型の用意

OpenJijを使った一定温度でのモンテカルロシミュレーション

スピングラス秩序変数の計算

参考文献

単位円上の有理点のなす群と、p=a^2+b^2の形に表せる素数

choosenim のインストール

射影多様体の積

この記事の概要

Aircraft Loading Optimization

モデル化

目的関数

制約

Payload Limits (荷重制限)

Center of Gravity Limits (重心制約)

Shear Limits (シアに対する制約)

Total QUBO

不等式制約の計算法

最適化計算の実装

数値実験結果

古典ソルバー

Average Time-to-Solution

成功確率

シア制約

量子ソルバー

結言

参考文献

文責

単位円上の有理点のなす群と、 $p=a^2+b^2$ の形に表せる素数