Maxwell方程式の表現論II(1)masslessスカラー場

The minimal representation of the conformal group and classic solutions to the wave equation
https://arxiv.org/abs/0901.2280
を読んで
Maxwell方程式の表現論
https://m-a-o.hatenablog.com/entry/20150621/p1
が正しくないことに気付いて、とりあえず、スカラー場の計算をした。
so(n+1,2)の極小表現の実現
https://vertexoperator.github.io/2021/05/05/hermite_BDI_minrep.html

論文では書かれてなかったけど、運動量空間での最低ウェイトベクトルを計算して、最低ウェイト条件の計算までをやった。スカラー場だけなので、4次元時空に限定することなく、一般の(n+1)次元時空での計算。運動量空間で計算する必要はないけど、式の形は、単純な形になる。

以前、間違えた理由は、最低ウェイトベクトルは、Lorentz不変になると思い込んでて、そうなるようにCartan部分環を選んだのだけど、Lorentz不変になるように選ぶのではダメらしい。生成子が全部コンパクトになるように選ぶのだけど、イマイチ、それで正しい理由が理解できてない。

スカラー以外の場については、今回計算したスカラー場を成分に持つ多成分場を作ればいけるはずで、前回と同じだと思うけど、まだ計算はしてない。

Lennard-Jones流体の分子シミュレーション(2)分子動力学vsメトロポリス法

以前、分子動力学の実装を作ったけど、諸条件、諸パラメータをこれと合わせることで、メトロポリス法と分子動力学で結果を比較できる。メトロポリス法を使ったことがなかったので、やってみようと思った。

分子動力学では、NVE一定の計算が、一番基本的になるのに対して、メトロポリス法では、NVT一定の計算が、一番基本的になる。メトロポリス法でも、NPT一定の計算をやる方法は考えられているし、分子動力学でも、barostatやthermostatを使えば、NVTやNPT一定の計算ができることにはなってるけど。また、メトロポリス法に不利な点として、輸送係数を計算する方法がないと思う。

全部実装して計算してから、気付いたけど、以下に書いたようなことは、Verletが、1967年の論文でやっていた。この論文は、Lennard-Jones流体の分子動力学計算を報告した最初の論文だそうで、この論文にさらっと書かれた数値積分法は、その後Verlet法と呼ばれるようになった

Computer "Experiments" on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules
https://doi.org/10.1103/PhysRev.159.98



メトロポリス法の論文は1953年に出ている。この論文は、与えられた確率分布に従う乱数列を生成するメトロポリスヘイスティングス法の最初の例が記述された論文として、物理以外の分野にとっても、重要な論文となった。

Equation of State Calculations by Fast Computing Machines
https://doi.org/10.1063/1.1699114

何故か、Wikipediaにまでページがある
Equation of State Calculations by Fast Computing Machines
https://en.wikipedia.org/wiki/Equation_of_State_Calculations_by_Fast_Computing_Machines

著者は5人いるけど、その一人Rosenbruthが書いた

Genesis of the Monte Carlo Algorithm forStatistical Mechanics
https://doi.org/10.1063/1.1632112

によると、統計力学の問題をターゲットにしようと決めて、中心となるアルゴリズムを考えたのは、Rosenbruthっぽい。また、最初は、普通に"分子動力学"をやろうと思ったが、計算機のスペックが厳しそうなので別の方法を考えたともある。論文で使われた計算機は、MANIACのようで、メトロポリスが開発を指揮したらしい。


メトロポリスらの論文で実際に計算してるのは二次元剛体球模型だけど、殆どの部分は、模型の詳細とは関係ない。ついでなので、論文を眺めた。

古典統計力学では、アンサンブル平均の計算が基本的になる。アンサンブル平均は、相空間での積分
\bar{F} = \left( \displaystyle \int F(\mathbf{q} , \mathbf{p}) \exp(-E/k_{B} T)d^{M}p d^{M}q \right) / \left( \displaystyle \int \exp(-E/k_{B} T)d^{M}p d^{M}q \right)
のこと。k_{B}ボルツマン定数Eはエネルギー(ハミルトニアン)で、相空間上の関数。Fは、相空間上の任意の関数。Tは温度。相空間の次元は2Mで、空間次元が3、粒子数がNなら、M=3Nであるが、メトロポリスは、空間次元が2で考えているので、M=2N。例えば、F=Eの場合は、アンサンブル平均は、内部エネルギーを与える。

よく出てくるアンサンブル平均の計算では、運動量pに関する積分は解析的に実行して消去できることが多いので、位置に関する積分のみが残る。例えば、ちょっと計算すれば分かるけど、内部エネルギーUは、
U = \dfrac{M k_{B} T}{2} + \left( \displaystyle \int e^{-\phi / (k_{B} T)}d^M q \right) / \left( \displaystyle \int \exp(-E/k_{B} T)d^{M}p d^{M}q \right)
のように書ける。相互作用項がない場合は、理想気体の結果として、統計力学の教科書に書いてあると思う。

運動量が動く範囲は力学的に決まるが、位置座標の範囲は、人間が決める。特に理由が無い限り、立方体領域にするのが簡単なので、普通、そうすると思う(周期境界条件を課すことが多いので、正確にはトーラス)。立方体の体積Vと、粒子数N,温度Tを外部パラメータとして指定することになるので、正しくカノニカル・アンサンブルの計算を行うことになる。

計算したいものは単純だけど、粒子数が100程度でも、M=200とかM=300なわけで、低次元で有効な数値積分アルゴリズム(台形公式やガウス求積)を使って評価するのは難しい。しかし、ポテンシャルエネルギーが大きい配置(例えば、二粒子が極めて近接してるような)は、積分に小さな寄与しかしないので、現代で言う棄却サンプリングによって、粒子の配置を\exp(-\phi /k_{B} T)に比例する確率で発生させ、このような配置に対して、F(\mathbf{q})の平均を取れば、アンサンブル平均を効率的に計算できる。

ところで、適当な統計力学の教科書を開くと、hをPlanck定数として、分配関数
Z = \dfrac{1}{N!h^{M}} \displaystyle \int \exp(-E/k_{B} T)d^{M}p d^{M}q
を計算すれば、自由エネルギーとか、エントロピーとか圧力とかの熱力学諸量は、分配関数の対数を取ったり微分したりする数学的操作で全部計算できるよ〜と書いてある(今の場合、分配関数Zは、N,V,Tの関数)。が、メトロポリス法では、分配関数そのものは計算できないので、こういう話は全然役に立たない。統計力学の教科書は破り捨てよう。

メトロポリス法で計算できるのは、アンサンブル平均なので、熱力学量をアンサンブル平均に帰着する必要がある。メトロポリスらの論文では、圧力を扱っている(論文のタイトルにもある通り、状態方程式を調べたいというのが動機だったので)。圧力は、統計力学の教科書的な定義では
P = \dfrac{1}{\beta Z} \dfrac{\partial Z}{\partial V}
である。この式を形式的に計算すると、ハミルトニアンを体積で微分した関数のアンサンブル平均という、よく分からないものが得られる。体積は、Hamiltonianの変数ではないので、素朴な意味では、この式は意味をなさない。

(メトロポリスらの論文では使われてないけど)たまに見かける議論として、体積の代わりに、Hamiltonianの位置と運動量のスケール変換(正準変換になるように、位置をr倍したら運動量は1/r倍する)を考えて、r=1での微係数を計算をすると、ビリアル圧力が得られる。正確には(空間次元を3とすれば)V=r^3 V_{0}に対して、\dfrac{\partial}{\partial V} = \dfrac{1}{3r^2 V_{0}} \dfrac{\partial}{\partial r}なので、体積と空間次元で割る必要がある。形式的な議論に見えるけど、配位空間内の積分領域でなく、相空間内の積分領域に作用する一径数(正準)変換群を考えて、この変換群で、分配関数がどう変化するのか見ていて、やってることは正当である。正準変換で、相空間の測度が不変に保たれることが重要。

ビリアル圧力は、位置のみの関数で書けるので、メトロポリス法で計算するのに適している。分子動力学でも、ビリアル圧力は、瞬時圧力を計算するのに一般的に使われている。分子動力学の場合は、マクロな圧力は、瞬時圧力の時間平均として得る。

スケール変換に限らず、もっと一般に座標の線形変換を、相空間の正準変換に拡張できる。具体的には、A \in GL(n,\mathbf{R})に対して、 \left(\begin{matrix} A & 0 \\ 0 & (A^T)^{-1} \end{matrix} \right) \in Sp(2n,\mathbf{R})を全粒子に一様に施す。A単位行列の定数倍であれば、さっきのスケール変換の状況。GL(n,\mathbf{R})には、直交変換群O(n)が含まれているが、これは、系全体を回転する(あるいは鏡に映す)だけなので、なんの影響も及ぼさない。直感的には、この線形変換は、粒子のいる立方体領域を変形する操作に相当する。

無限小線形変換で考えると、
 r_{\alpha} \mapsto r_{\alpha} + \sum_{\beta} \epsilon_{\alpha \beta} r_{\beta}
の形で、歪みテンソル\epsilon_{\alpha \beta}は対称としていい(反対称成分は、無限小回転に相当)。Hamiltonianにこの変換を施して、歪みテンソル微分すると、応力テンソル(virial stress)の成分が計算できる。virial stressは、運動量流束密度と見ることもでき、粘性係数を計算するのにも使われる。

他に得られる物理量としては、内部エネルギーの熱的揺らぎから、定積比熱が得られる(はず)。内部エネルギーやエントロピーは、直接測定できる量ではないが、定積比熱は測定できる量だけど、多原子分子だと、分子内振動の自由度からの寄与が正しく評価されない気がする(これは古典統計力学では、どうやっても計算できないので)。単原子分子でも、分子間相互作用からの寄与が、どれくらい正しく評価されるものか分からない。

エネルギーの二乗のアンサンブル平均は、\beta=1/(k_{B}T)として、以下のように計算できる
\bar{E^2} = \dfrac{M^2+2M}{4\beta^2} + \left(\dfrac{M}{\beta} \displaystyle \int \phi e^{-\beta \phi} dq + \displaystyle \int \phi^2 e^{-\beta \phi} dq \right) / \left( \displaystyle \int \exp(-\phi/k_{B} T) dq\right)

エントロピーは、相空間上の関数のアンサンブル平均として直接計算する方法はないと思う。エントロピーが得られると、ヘルムホルツの自由エネルギーF=U-TSがわかり、Z=e^{-F/k_{B}T}によって、分配関数が計算できてしまう。


以上は模型によらない一般論だけど、2次元剛体球(rigid sphereと書いてるけど、現在は、hard sphereの方が、通りがいいと思う。2次元なので、circleあるいはdiskだけど、rigid circleとかhard diskとは言わない)の計算では、もう少し、固有の工夫が要る。

剛体球模型は、粒子同士が重なると、ポテンシャルエネルギーが無限大、重ならなければ、ポテンシャルエネルギーは0というモデルなので、\exp(-\phi /k_{B} T)は、0か1しか取らない。従って、棄却サンプリングでは、重なる配置は即棄却、そうでない配置は即採用と比較的アルゴリズムが単純になる。

一方、ビリアル圧力の計算には、粒子に働く力を知る必要があり、これは、通常の場合、ポテンシャルエネルギーの微分で計算されるが、剛体球模型では、許可される配置に於いて、ポテンシャルエネルギーは常に0であって、そのまま適用できない。剛体球模型で力が働くのは、粒子同士が衝突した一瞬であって、メトロポリスの論文では、この衝突力の評価を、動径分布関数(ある粒子から一定距離にある粒子密度)から評価する式を与えている。動径分布関数も、解析的に決めることはできないが、数値的に外挿することで評価したと書いてある。

あと、メトロポリス法では、乱数生成が必要だけど、疑似乱数生成アルゴリズムに何を使ったかは書かれてない。当時は、von Neumann考案のmiddle-square methodと、Lehmerの方法というものが知られていたっぽい。Lehmer乱数は、線形合同法の特殊ケースでもある。多分、どっちかの方法を使ったんだろう。現在は、どっちも使われてないアルゴリズム


二次元剛体球模型は、実験との比較ができないtoy modelなだけでなく、計算機が速くなった今となっては、本質的じゃない部分で却って面倒が多い。というわけで、普通に、3次元のソフトポテンシャル模型を考え、ポテンシャルのモデルとしては、ポピュラーなLennard-Jonesポテンシャルを使ってメトロポリス法を実装した。希ガス気体は理想気体と変わらなすぎるので、希ガス液体(今回は希ガスアルゴン)の計算をやる。で、冒頭に書いたように、分子動力学や実測結果と比較する。

分子動力学の計算では、全粒子対で相互作用計算をすると、粒子数Nに対して、O(N^2)なので、カットオフ距離を設定して、計算量をO(N)にしている。メトロポリス法は、最初に一回だけ全粒子対の計算をするのみなので、カットオフ距離は別に必要ない。ただ、今回は、条件を揃えて計算するために、カットオフ距離を設定できるようにした。

Lennard-Jonesポテンシャルで計算するのに、カットオフを用いる場合、内部エネルギーや圧力に、カットオフで削った部分の長距離補正項を足す必要がある。g(r)を動径分布関数とすると、内部エネルギーと圧力は
U = \dfrac{3}{2}Nk_{B}T + 2\pi N \rho \displaystyle \int_{0}^{\infty} \phi(r) g(r) r^2 dr
p = \rho k_{B} T - \dfrac{2 \pi \rho^2}{3} \displaystyle \int_{0}^{\infty} \phi'(r) g(r) r^3 dr
などと書ける。\rhoは粒子数密度。長距離補正は、カットオフ距離より遠い部分の積分
 U_{LRC} = 2 \pi N \rho \displaystyle \int_{r_c}^{\infty} \phi(r) g(r) r^2 dr
 p_{LRC} = -\dfrac{2 \pi \rho^2}{3} \displaystyle \int_{r_c}^{\infty} \phi'(r) g(r) r^3 dr
という形。Lennard-Jones流体に対する補正項の計算では、g(r) \approx 1と近似するのが一般的。こうすると、積分が計算できる。圧力の長距離補正は、密度の二乗で効いてくるので、気体の時はともかく、液体では、かなり大きな影響を持ちうる。

ついでに、一般のvirial stressに対する長距離補正は
\sigma^{ij}_{LRC} = \dfrac{\rho^2}{2} \displaystyle \int_{|\mathbf{r}| \gt r_c} \frac{r_i r_j}{|\mathbf{r}|} \phi'(|\mathbf{r}|) g(|\mathbf{r}|) d\mathbf{r}
という形。
p_{LRC} = -\dfrac{1}{3}(\sigma^{xx}_{LRC} + \sigma^{yy}_{LRC} + \sigma^{zz}_{LRC})
が成り立つ。。

カットオフを入れると、ポテンシャルエネルギーが不連続になるので、微分すると、力が発散するのを気にする人もいるらしい。エネルギーが連続になるように修正したポテンシャルエネルギーを使うこともある。
cf)Numerical Experiments on the Stochastic Behavior of a Lennard-Jones Gas System
https://doi.org/10.1103/PhysRevA.8.1504

Verletは気にしなかったようなので、これは、気にしなくてもいい気がする(というか、入れてない)


アルゴンのLennard-Jonesパラメータを使って、粒子数10万、密度1407(kg / m^3)、温度85(K)あたりで、液体と思われる平衡状態を生成して、これを初期状態に、カットオフ距離を変えつつ、分子動力学とメトロポリス法で圧力を計算した結果及び圧力補正項は以下のようになった。rcutがカットオフ距離で、sigmaはLennard-Jonesパラメータのσ。通常、rcutはsigmaの2.5〜10倍程度の範囲で設定するっぽい

rcut / sigma p(atm)@MD p_LRC(atm) T(K)@MD p(atm)@MC
2.5 1.93 -312 85.26 -2.57
5.0 9.80 -39.1 85.02 11.98
7.5 11.25 -11.6 84.95 12.55
15.0 -1.45 12.90
30.0 -0.181 10.43

メトロポリス法の温度は、85(K)で、積分点を100万個サンプリングした。カットオフ距離が大きくなると、カットオフによる計算高速化の意味がなくなるので、分子動力学では、大きなカットオフはやってない。圧力の計算値は、長距離補正込みで、p_LRCは圧力の長距離補正。今は液体を扱ってることもあって、標準的な範囲のカットオフ距離だと、長距離補正の影響は、結構でかい。

Verletの論文では、rcut/sigmaは、2.5や3.3を使ったと書いてある。Verletの論文を読んだのが、この計算した後だったせいで、Lennard-Jonesポテンシャルのパラメータは、Verletの使った数値とはかすかに異なる(例えば、私はσを3.41Åにしたけど、Verletは、3.405Å)

Verletの論文で、rcut/sigmaが小さめなのは、粒子数が少ないせいもある。私が、粒子数10万を選んだのは、特に根拠があったわけでもないけど、今回くらいの条件で、長距離補正が支配的にならないようにしたいと思うと、粒子数10万で、多すぎるとも言えないように思う。


で、実験と比較したい。以下は、やや古い論文ではあるが、測定結果に基づいて、液体アルゴンの温度と圧力から密度を計算するフィッティング式が書かれている。

Density of liquid nitrogen and argon as a function of pressure and temperature
https://www.sciencedirect.com/science/article/abs/pii/0031891460900422

85(K)での測定結果はないが、この式に従って計算すると、
85(K),1(atm)での密度は1409(kg/m^3)
85(K),10(atm)での密度は1411(kg/m^3)
85.25(K),1(atm)での密度は1407(kg/m^3)
85.25(K),10(atm)での密度は1410(kg/m^3)
85.5(K),10(atm)での密度は1408(kg/m^3)
86(K),10(atm)での密度は1405(kg/m^3)
などとなるので、計算結果は、全くの的外れってことはなさそう。

最初は、(温度は85Kで)、1400(kg/m^3)、1405(kg/m^3)などで計算したが、いずれも負の圧力が実現されたので、次に1407(kg/m^3)を選んで計算をしたら、正の圧力が得られたので、それから文献調査をした。同じ温度で、1(atm)と10(atm)でも、密度が2~3(kg/m^3)くらいしか変わらない割に、そこそこの値を出していると思う。

Lennard-Jonesポテンシャル自体が、そもそも雑な近似(パラメータの決め方も正しい方法とかないし)なので、実験結果との一致がしょぼくても仕方ないけれど、もう少し、ちゃんと測定結果との比較をしてみようと思って、以下の論文の結果と対応する計算結果を出した。論文では、温度と密度と圧力が与えられてるので、温度と密度を指定して、圧力を計算する。

Measurement and correlation of the (pressure, density, temperature) relation of argon I. The homogeneous gas and liquid regions in the temperature range from 90 K to 340 K at pressures up to 12 MPa
https://doi.org/10.1006/jcht.1994.1048

T(K) density(kg/m^3) p(atm)@exp p(atm)@MD T(K)@MD p(atm)@MC
90 1381.34 9.65 25.97 90.20 26.81
100 1320.44 19.58 25.59 100.05 22.50
110 1250.17 19.76 10.49 109.68 10.42

粒子数は10万。
MDのカットオフ距離は、σの5倍。
メトロポリス法のカットオフ距離は、σの30倍で、積分点は100万点。

分子動力学で初期平衡状態を作って、NVE一定の分子動力学とメトロポリス法による計算を行った。とりあえず、圧力を高い精度で予測するには至ってない。


今回の計算に使用したコードを、gistに置いておく。
https://gist.github.com/vertexoperator/49737c13e432bf13e26c851683a6821d
パラメータとかは、mainの中でハードコードしてるのを手動で書き直して、使ってたので、使い勝手はよくない。


差し当たって、適当に検索して出てきた数値を使ってシミュレーションした。LJパラメータは、文献によって、そこそこバラツキがあるので、仮に、たまたま良いパラメータを引いて、ある範囲の条件下で、熱力学量がよく計算できても、パラメータいじれば合うよねという感じがしなくもない。なので、Lennard-Jonesパラメータにも根拠が欲しい。結局、実験的に決めるより、第一原理的に計算する方法があればベスト思われ、それには量子力学が必要だけど、随分昔にやった

量子力学はLennard-Jonesポテンシャルを再現できるか
https://m-a-o.hatenablog.com/entry/20121217/p2

を眺めると、いまいちLJパラメータの計算値と実験値は合わないと書いている。今適当に検索すると、このにっきを書いた翌年に

Theoretical Investigation of the Dispersion Interaction inArgon Dimer and Trimer
https://doi.org/10.1016/j.procs.2013.05.247

という論文が出てて、CCSD(T)/aug-cc-pVTZで計算した場合、Ar-Arの平衡分子間距離が3.75Åあたりになると書いてある。普通のLJポテンシャルで、エネルギー最小点が一致するようにパラメータを決めると、σが3.34Åくらいに相当。この論文で引用されてた

A highly accurate interatomic potential for argon
https://doi.org/10.1063/1.466051

を見ても、平衡分子間距離は、ほぼ同じ値。Verletの使ったσが、3.405Å。εについては、Fig1(a)のグラフを見ると、平衡分子間距離におけるエネルギーが、(無限遠に置いた時を基準として)約-100(/cm)=-0.0124(eV)。上のシミュレーションで使った値は、eV換算で、0.01034(eV)。

LJポテンシャルの式の形は単純すぎるので、別の関数形を使うことも考えるべきなのかもしれない。安直にパラメータを増やすのはよくないが、双極子の次は、四重極も...と考えるのは自然なことで、実際に"精密化学"界隈では、r^(-8)やr^(-10)に比例する項も計算されるのが普通らしい。適当に検索した所、London分散力の延長として、四重極を考慮した人は、1938年には存在している
Quadrupole Contributions to London's Dispersion Forces
https://doi.org/10.1063/1.1750184
古典分子動力学で、そのようなポテンシャルを使っていけない理由は別にない。

このへんの問題は、おいおい追求するとして、とりあえず、CCSD(T)/aug-cc-pVTZ計算の追試のため、Gamessに以下のような入力ファイルを食わせた所、核間距離は、3.8008895(ANGS.)で、TOTAL ENERGYは、-1054.1027355965(hartree)となった。論文より平衡分子間距離が少し大きいけど、まぁ許容範囲と思う

! rungms Ar2.inp  > Ar2.log
 $contrl runtyp=optimize CCTYP=CCSD(T) NUMGRD=.T. ISPHER=1 coord=unique $end
 $system TIMLIM=20000 MEMORY=10000000 $end
 $basis gbasis=ACCT $end
 $guess guess=Huckel $end
 $data
argon dimer
C1
 Ar  18 0.0 0.0 0.0
 Ar  18 0.0 0.0 3.7
 $end

同じようにして、希ガス二量体に対するCCSD(T)/aug-cc-pVTZレベルの計算をした結果を、まとめると、以下のようになった。

Rm(Angs.) Em(hartree) E0(hartee) (E0-Em)/k_B(K) R_nn(Angs.) 融点(K) 沸点(K) 融解熱/R(K)
Ne 3.0916867 -257.6254676649 -257.6252957319 54.293 3.13 24.56 27.07 40.29
Ar 3.8008895 -1054.1027355965 -1054.1023022117 136.853 3.76 83.80 87.30 141.92
Kr 4.1146791 -5504.4997943600 -5504.4992507992 171.644 4.01 115.79 119.93 197.25

Rm:平衡分子間距離(計算値)
Em:平衡分子間距離におけるエネルギー(計算値)
E0:核間距離を100Åにした時のエネルギー(計算値)
k_B:ボルツマン定数
R_nn:希ガス固体の最隣接原子間距離(Kittelの教科書に載ってるらしいのを孫引き)
融点:出典Wikipedia(多分、1気圧の数値?)
沸点:出典Wikipedia(多分、1気圧の数値?)
融解熱:出典Wikipedia
R:気体定数

組合せ範疇文法の漸進的構文解析

随分と昔に、CCG(combinatory categorial grammar/組み合わせ範疇文法)のパーサを作った。
https://m-a-o.hatenablog.com/entry/20160614/p3

この時の実装は、CYK法に近いもので、ボトムアップ構文解析である。

ところで、人間が自然言語を処理する時は、単語を先頭から順に処理していってると思われる。普通のCCGに於いて、例えば、I saw a cat.という文は標準的には、以下のような導出木によってparseされる。

I    saw          a          cat
--  ---------   ------      -----
NP  (S\NP)/NP    NP/N         N
                ----------------- (/-application)
                       NP
    ----------------------------- (/-application)
               S\NP
--------------------------------- (\-application)
               S

この場合、文章が最後まで提示されないと構文解析を開始できないが、英語を知ってる人は、I saw a ...まで提示された時、次に目的語が来ることを期待していて、最後まで文が提示されなくても、ある程度の構文解析を行い、かつ、文章が完成していないことも理解している。

CCGには、単純な左適用、右適用以外にも、色々な導出規則があって、特にtype raisingという規則では、最終的に得たい統語範疇がSで、途中までの統語範疇がXの場合、残りの部分の統語範疇はS\Xだろうから、統語範疇Xを、S\Xに右適用して統語範疇Sが得られるS/(S\X)に置き換えることが許される。元々は、type raisingがあると、単語の統語範疇を追加することなく、導出木を作れるケースがあるという理由で、導入された規則だと思うけど、続きの文章に対する予測を行うのに類似した作用を定式化していると考えることもできる。

type raisingを使うと、上の文章の場合、以下のように、前から順に導出木を作っていくこともできる。

   I        saw          a       cat
-------   ---------   ------    -----
  NP      (S\NP)/NP    NP/N       N
-------- (type raising)
S/(S\NP)  
-------------------- (composition)
        S/NP
----------------------------- (composition)
            S/N
--------------------------------------- (/-application)
                   S

CCGに於いて、この2つの導出木は、等価な構文解析結果を与える。

つまり、統語範疇X/YとX\Yを持つ単語は、型Y->Xを持つ関数と解釈され、上の2つの導出木は型Sを持つ値を作っている。コンビネータBをB f g = λx.f (g x)で定義すると、これは関数の合成操作で、単語wの"意味"が[ [w] ]で与えられるとすると、後者の導出木の値は以下のように計算できる。

(B (B (λf.f [[I]]) (λx.[[saw]] x)) [[a]]) [[cat]] 
=>  (B (λf.f [[I]]) (λx.[[saw]] x)) ([[a]] [[cat]])
=> (λf.f [[I]]) ([[saw]] ([[a]] [[cat]]))
=> [[saw]] ([[a]] [[cat]]) [[I]]

これは最初の導出木から得られる値と等しい。この柔軟性のために、CCGには同値な導出が無数に存在することになり、枝刈りしないと大変なことになる。追加していいCCGの導出規則は、上のような単純型付ラムダ計算の体系を破綻させないことが必要条件になるが、十分条件は良くわからない(多分、言語ごとに使っていいコンビネータのセットが多少違う?)ので、不足がないという保証はない。

細かく言えば、統語範疇も、もう少し細分化しないと、現在の標準的なCCGでは、"him are I."みたいな文章も通ってしまうので、非文を弾くということには、比較的甘い傾向がある。

前から順に処理していけるかは、実用的観点からも気になる話。機械学習を使ったNLPでも、文章は、単語列を先頭から順に処理していくことで結果を得るのが普通だし、
Unsupervised POS Induction with Word Embeddings
https://arxiv.org/abs/1503.06760
という論文では、窓幅1でword2vecで学習したベクトルからHMMを使って、品詞の教師なし学習をやって、従来法より良い結果が得られると書いている。

窓幅1なので、semanticな情報はほぼ皆無で、統語的な情報のみを含んでるだろうと感覚的には思うが、一方で、そもそも、原理的に窓幅1のword embeddingでは捉えきれない統語情報がないのか疑問に思う。もし、完全にincrementalなparserが作れるのならば、(word2vecで出来るかどうかはともかく)原理的には窓幅1で統語情報を全て得られる可能性はある。

そういうわけで、incremental parserが可能なのかというのは、興味ある問題である。

Incremental combinatory categorial grammar and its derivations
https://dl.acm.org/doi/10.5555/1964799.1964809
では、incremental combinatory categorial grammar(ICCG)と呼んでる。が、結果は完璧でない。

Incremental Derivations in CCG
https://www.aclweb.org/anthology/W12-4623.pdf

は、もう少し詳しく議論している。

3つの連続する統語範疇から一つの統語範疇Xが導出される場合、
(X/Y)/Z Z Y
Z (X/Y)\Z Y
X/Y Y/Z Z
X/Y Z Y\Z => X/Y Y/(Y\Z) Y\Z
Y/Z Z X\Y
Z Y\Z X\Y
Y (X\Y)/Z Z => X/(X\Y) (X\Y)/Z Z
Y Z (X\Y)\Z => X/(X\Y) (X\Y)/((X\Y)\Z) (X\Y)\Z
Y X/Z Z\Y => ???
などのパターンがある。いくつかは自明にincremental parseが可能で、3つはtype raisingを必要とするが、incrementalにparseできる。最後のパターンをincrementalにparseするために、論文では、Geach ruleという新たな規則を追加している(Peter Geachという人がいたらしい)。combinator birdでは、bluebirdに相当し、
A/B : f => (A/T)/(B/T) : λg.λt.f(g t)
のように意味づけすることができる。論文には書かれてないが、
A\B : f => (A\T)\(B\T) : λg.λt.f(g t)
をBackward Geach ruleと(私が勝手に)呼ぶことにする。

1つ目の統語範疇にtype raisingを使い、2つ目にGeach's ruleを使うと、
Y X/Z Z\Y => X/(X\Y) (X\Y)/(Z\Y) (Z\Y)
のような形になって、incrementalにparseできる。

後者の論文では、incrementalにparseするのが難しい例として、以下のようなthat節を挙げている。最も標準的なCCG導出木は以下の形だろう。

the woman        that       every       man          saw        laughed
---------   -------------  -------     ------     ----------    --------
   NP       (NP\NP)/(S/NP)   NP/N         N        (S\NP)/NP      S\NP
                            ------------------
                                   NP
                            ------------------
                                S/(S\NP)
                            ---------------------------------
                                          S/NP
             -------------------------------------------------
                              NP\NP
--------------------------------------------------------------
                            NP
----------------------------------------------------------------------------
                                        S

thatは、目的語を残した文(今の場合、every man sawで、統語範疇はS/NP)を受け取って、名詞句を後ろから修飾し新たな名詞句を作り出す(NP\NPの部分)ので、(NP\NP)/(S/NP)のような統語範疇は、文法解釈的に自然である。

一応、Backward Geach ruleと前者の論文では導入されているargument swap(関数の引数を置換する操作で、(S\NP)/NPを(S/NP)\NPに置き換える)を使い、更にtype raisingを一つの単語で2回やるという反則っぽいことをすると、incrementalにparseできなくはない。

the woman        that              every                   man                                 saw        
----------  -------------  ----------------      ---------------------------                 ---------- 
   NP       (NP\NP)/(S/NP)        NP/N                       N                                (S\NP)/NP 
----------
NP/(NP\NP)
-------------------------
          NP/(S/NP) 
                          ---------------------                 
                     (S/NP)/((S/NP)\(NP/N))
----------------------------------------------   ------------------------------ (type raising)
            NP/((S/NP)\(NP/N))                              NP\(NP/N)
                                                ----------------------------------------------- (type raising)
                                                ((S/NP)\(NP/N))/(((S/NP)\(NP/N))\(NP\(NP/N)))
---------------------------------------------------------------------------------------------- --------------  (argument swap)
                        NP/(((S/NP)\(NP/N))\(NP\(NP/N)))                                          (S/NP)\NP
                                                                                                ------------- (Bwd Geach rule)
                                                                                               ((S/NP)\(NP\N))\(NP\(NP/N))
------------------------------------------------------------------------------------------------------------------------
                                                      NP

こういう無理矢理っぽいことをしても、semanticsの一貫性を保証できるのがCCGの良い点ではあるが、複雑すぎるので、人間の脳が、本当に、こんなことをしてるのか(出来るのか)という点は、疑問である(私は、この導出木を見つけるのに2時間くらい考えた)

実装上は、何回までtype raisingを許すのかという問題が出てくる。全ての単語とparse途中で現れる中間統語範疇が複数回type raisingする可能性があると考えだすと、場合によっては組み合わせが膨大になって大変なことになりそう。

しかし、incrementalにparseする場合、一見、文章として完成していても、続きがあるかもしれないという問題があるので、この困難は本質的には回避不能という気もする。

例えば、He ranという単語列は、2単語で一応完結しており、構文的には、これで一つの文でありえるが、He ran quickly.という文章では、"He ran"の部分の統語範疇は、S/( (S\NP)\(S\NP) )になるべきで、He ran quickly and stopped suddenly.という文章では、"He ran"の部分には、また別の統語範疇を付与する必要がある。

incrementalにparseする場合、"He ran"という2単語が提示された段階で、どこまで文章が続くかは全く未知であって、”He ran"に付与する統語範疇として想定しなければならないパターンは、おそらく無限個ある。そうであれば、完全にincrementalなCCG parserは、一単語目で足止めされることになって、計算できない。

incremental parserが不可能であっても、任意の有限な長さの文に対して、incrementalな導出木が必ず存在するという可能性はある。現実的には、人間が話したり書いたりする文章の長さは、限りがあるし、無限のパターンを記憶することもできないので、人間の脳は、何らかのヒューリスティックで個数を制限して、incrementalな導出木を計算してるのかもしれない。

type raisingは無限のパターンを作り出すための仕組みでもあるし、また、実装上は、多相型のように扱うのがよいと思われ、いつ終わるか分からない文章で統語範疇をfixするのを先送りする役割を果たす。と考えられる。

 He                ran                     quickly
----------     -----------------        -------------
 NP                S\NP                  (S\NP)\(S\NP)
------------   -----------------
∀X.X/(X\NP)    ∀Y.Y/(Y\(S\NP))
-------------------------------- (Y:=X\NP)
     ∀X.X/((X\NP)\(S\NP))
------------------------------------------------------- (X:=S)
                          S


ついでに、1997年にLambekが提案したpregroup grammarというものを読んだ。この文法では、argument swapは自動的に入っている。
Type Grammar Revisited
https://link.springer.com/chapter/10.1007/3-540-48975-4_1

pregroup grammar
https://ncatlab.org/nlab/show/pregroup+grammar

pregroupは、半順序モノイドで、ある性質を持つleft adjointとright adjointを持つ代数構造。順序構造については、 a \leq bの代わりに、簡約的なニュアンスで、 a \to bと書かれることがある。left adjointとright adjointは a^{l} a \to 1 , 1 \to a a^{l}及びa a^{r} \to 1 , 1 \to a^{r} aを満たす。

逆元のようではあるが、a^{l} a = 1というわけではない。実際、a^{l} a = 1だと左逆元そのものであるが、初歩的な代数学の教科書の最初に書いてある通り、左逆元と右逆元は一致してしまうので条件を緩めている。群でないpregroupが存在することは、Lambekの論文やncatlabのページに例が載っている(どっちも同じ例)。明らかに、群でないpregroupは、積について非可換である。

left adjointとright adjointは、通常の逆元が持つのと類似の性質を多く持ち、
1^{l} = 1^{r} = 1
(ab)^{l} = b^{l} a^{l}
(ab)^{r} = b^{r} a^{r}
が成り立つ。また、a^{ll} = (a^{l})^{l}などと表記すると、
 a^{ll} a^{l} \to 1 \to a a^{l}
なので、a^{ll}aが等しいとは限らない。a^{rr}も同様。a^{lr}=a^{rl}=aは成立する。

pregroupはleft dualとright dualが一致しないタイプのrigid category(autonomous category)の例を与えていて、そのような圏の別の例は、Joyal and Streetの論文
AN INTRODUCTION TO TANNAKA DUALITY AND QUANTUM GROUPS
http://maths.mq.edu.au/~street/CT90Como.pdf
の§9で見ることができる。大抵の数学の文献では、left dualとright dualが同型になるような圏しか扱われない(と思う)。


pregroup grammarは、combinatoryでない範疇文法とセットで議論されることが普通なので、\-application ruleが
Y X\Y => X
の形でなく、
Y Y\X => X
と書かれている。紛らわしい上に、面倒だけど、多分、Steedmanが悪い。ここでは、組み合わせ範疇文法の話を書いてきたので、\は、組み合わせ範疇文法の記法に合わせる。

で、A/Bは、\mathrm{A} \cdot \mathrm{B}^{l}に置き換えて、A\Bは、\mathrm{B}^{r} \cdot \mathrm{A}に置き換える。後は、普通の群やモノイドのように思って計算する。合流性は成立してないので、a^{l}aa^{r}a(a^{l}a)(a^{r}a) \to a^{r}aであると同時にa^{l}(aa^{r})a \to a^{l}a \to 1でもある。このへんがparseをめんどくさくするが、CFGなどと同様、O(N^3)でparseするアルゴリズムが知られている。

例えば、"I saw a cat"は

\mathrm{NP} \cdot (\mathrm{NP}^{r} \cdot \mathrm{S} \cdot \mathrm{NP}^{l}) \cdot (\mathrm{NP} \cdot \mathrm{N}^{l}) \cdot (\mathrm{N}) \to \mathrm{S}

と計算される。他動詞"saw"の統語範疇(S\NP)/NPが\mathrm{NP}^{r} \cdot \mathrm{S} \cdot \mathrm{NP}^{l}にマップされて、argument swapを行った統語範疇(S/NP)\NPの場合と一緒になる。"I saw a ..."という途中までの文だと

\mathrm{NP} \cdot (\mathrm{NP}^{r} \cdot \mathrm{S} \cdot \mathrm{NP}^{l}) \cdot (\mathrm{NP} \cdot \mathrm{N}^{l}) \to \mathrm{S} \cdot \mathrm{N}^{l}

となる。

the woman that every man sawという句は、(式が長くなるので)"the book"と"every man"の統語範疇をNPとして、

\mathrm{NP} \cdot (\mathrm{NP}^{r} \cdot \mathrm{NP} \cdot \mathrm{NP}^{ll} \cdot \mathrm{S}^{l}) \cdot \mathrm{NP} \cdot (\mathrm{NP}^{r} \cdot \mathrm{S} \cdot \mathrm{NP}^{l}) \to \mathrm{NP}

となる。

type raisingも整合的になる。例えば、A => T/(T\A)の場合、\mathrm{T} \cdot (\mathrm{T}^{l} \cdot \mathrm{A}^{rl}) = \mathrm{T} \cdot \mathrm{T}^{l} \cdot \mathrm{A}となるが、A \to \mathrm{T} \cdot \mathrm{T}^{l} \cdot \mathrm{A}は、正しいので、問題ない。

Geach ruleとBackward Geach ruleもpregroup grammarで解釈できる。A/B => (A/T)/(B/T)は \mathrm{A}\mathrm{B}^{l} \to \mathrm{A} (\mathrm{T}^{l} \mathrm{T}^{ll}) \mathrm{B}^{l} = (\mathrm{A}\mathrm{T}^{l})(\mathrm{B}\mathrm{T}^{l})^{l}から従う。Backward Geach rule A\B => (A\T)\(B\T)も\mathrm{B}^{r}\mathrm{A} \to \mathrm{B}^{r} \mathrm{T}^{rr} \mathrm{T}^{r} \mathrm{A}から従う。

組み合わせ範疇文法のコンビネータで、pregroup grammarが処理できなさそうのもある。backward crossing composition
Y/Z:g X\Y:f => X/Z:λz.f(g z)
は、英語では必要なことになってるけど、\mathrm{Y}\mathrm{Z}^{l}\mathrm{Y}^{r}\mathrm{X}\mathrm{X}\mathrm{Z}^{l}で、\mathrm{X}\mathrm{Z}^{l}の順序が入れ替わっているので、pregroup grammarでは困りそう。英語では、backward crossing compositionは、heavy NP shiftを説明するために、よく使われる。

ただ、英語でもbackward crossing compositionの無条件な適用は、非文も生み出すと考えられているので、CCGを拡張してslash typeを導入するなどの方法(Modalized CCGとか呼ばれる)で解決を図ろうとする場合もある。

言語によっては、他にも、pregroup grammarで処理できないコンビネータが必要かもしれない。


あと、incrementalにparseすることを考えると、組み合わせ範疇文法もpregroup grammarでも、言語固有の知識を追加で積極的に使わないと効率が悪そうではある。例えば、英語で、名詞句NPが連なることはない(I, Geras, am a fixed point in time.みたくカンマ区切りで列挙するのは別)けど、組み合わせ範疇文法やpregroup grammarの枠組みだけでは、それを捨てられない。

例えば、He she theyという単語列は、英語では出現しないだろうが、pregroup grammarだったら、続く単語列の統語範疇がX \cdot \mathrm{NP}^{l} \cdot \mathrm{NP}^{l} \cdot \mathrm{NP}^{l} \cdot \mathrm{S}で、X \to 1だったら、最終的にSに簡約される。そんなことがないかどうかは、pregroup grammarは教えてくれない。組み合わせ範疇文法だと、以下のようなことができる

  He          she            they
-------   -------------  ---------------
  NP           NP             NP
-------   -------------  ---------------
S/(S\NP)  ∀X.X/(X\NP)    ∀Y.Y/(Y\NP)
------------------------
    S/((S\NP)\NP)
-----------------------------------------
         S/(((S\NP)\NP)\NP)

中世の数理科学

Wikipediaによれば、mathematicsの語源は、ギリシャ語のmáthēmaで、knowledge, study, learningを意味してたそうだ。polymath(博学)という単語には、原義が残っている。そして、ラテン語や英語に於いて、1700年頃まで、mathematicsという単語は、占星術(時々、天文学)を指していたと書かれてるけど、これは正しくない。中世ヨーロッパに於いて、「数学」という単語が、どういう分野を指してたかは、それほど固定的でもない。

以下、ヨーロッパと書いてる時、西方のラテン語文化圏を念頭に置いている。また、中世は、地域ごとに、時代区分が微妙に異なるけど、以下では、西暦800〜1600年を指す時代区分のつもりで使っている。

アラビア数学とQuadrivium

12〜13世紀のカスティーリャ王国には、トレド翻訳学派という、アラビア語文献をヨーロッパの言語(ラテン語カスティーリャ語)に翻訳する仕事をしていた学者がいた。トレド翻訳学派には、ユダヤ人が多かったとも言われ、カスティーリャ王国では、1492年にアルハンブラ勅令が出されるまで、ユダヤ系学者が活動してた節がある。

トレド翻訳学派の一人にDominicus Gundissalinus(1115~1190)という人がいて、1140年に書いたDe divisione philosophiaeの中で、"数学は、算術、幾何、音楽、天文学、scientia de aspectibus(光学)、scientia de ponderibus(science of weights)、scientia de ingeniis(science of devices)の7つのartsを含んでいる"と書いている

Mathematica quoque universalis est, quia sub ea continentur septem artes, quae sunt arithmetica, geometria, musica et astrologia, scientia de aspectibus, scientia de ponderibus, scientia de ingeniis.

ラテン語のscientiaは英語のscienceと同根の語だけど、中世には、知識とか学問のようなニュアンスで、結構広い意味だったようである。格言「知識は力なり」の原文は、scientia est potentiaらしい。

astrologiaは、直訳では、英語のastrology/占星術だけど、中世ヨーロッパの大学では、リベラル・アーツという括りで、文法学・修辞学・論理学の3学と、Quadrivium(算術、幾何、天文学、音楽)という区分をよく使っていたことを踏まえると、天文学と解する方が適切に思える。Gundissalinusの分類では、Quadriviumより広い分野を含んでるけど、この分類は、多分、アラビア科学の学問論から来ている。Gundissalinusは、トレド翻訳学派の学者だったから、別に意外ではない。


アラビアの学者al Farabi(872?〜950?)がKitáb al-Ihsa' al-'ulum("Book of enumeration of the sciences")のChapter3で、"数学"('ilm al-Ta'alim)として分類している分野は、算術、幾何、天文学、音楽以外に、後のGundissalinusが、scientia de aspectibus、scientia de ponderibus、scientia de ingeniisと呼んだ分野に相当するものがあったらしい。

アラビア語Ta'alimは、instruction/teaching/educationなどの意味だそうで、ギリシャ語のmáthēmaと似た感じっぽい(現在のアラビア語で、「数学」を意味するのは別の単語)。ラテン語にしろ、アラビア語にしろ、この分野を指すのに、あまり適切な名前とは思えないんだけど、Quadriviumにしたって、内容を反映した命名ではないので、大した理由はなかったのかもしれない。al Farabiは、これらの分野を更にtheoreticalなものとpracticalなものに分けたとよく書いてあるけど、詳細は知らない。例えば、practicalな算術は、商業とかで使われる計算ということらしい。

Quadrivium自体は、少なくとも、古代ローマ帝国末期まで遡れるらしいけど、al Farabiの分類は、どっから来たものか、よく分からない。


ついでに、De divisione philosophiaeでは、"自然科学"(scientia naturalis)として、

  • scientia de medicina(薬学)
  • scientia de iudiciis(占術?)
  • scientia de nigromantia secundam physicam(science of necromancy!?)
  • scientia de imaginibus(science of image magic?)
  • scientia de agricultura(農学)
  • scientia de navigatione(航海術)
  • scientia de speculis(反射光学?)
  • scientia de alquimia(錬金術)

の8つを挙げている。

自然科学に分類されてるscientia de speculisと、数学に分類されてるscientia de aspectibusは、現代では、どちらも"光学"の一部に相当するが、多分、ユークリッドのCatopriticaとOpticaに対応したものと思われる。Catopriticaは、鏡で反射する光の理論で、Opticaは、視覚が目から出ている視線によって生じるという考えをベースにした視覚の幾何学のようなもの。aspectibusは、aspectusの複数形で、aspectusはvisionを意味するらしいから、scientia de aspectibusは、直訳では、science of visionsとなる。

但し、プトレマイオスの本"Optica"では、両方を扱ってるっぽいし、「光学の父」と呼ばれるIbn al-Haitham(965〜1040)は、視覚が物体から目に到達した光によって生じると理解してたようなので、この区別はGundissalinusの不見識という気もする。Gundissalinusの"自然科学"には、完全にオカルトな分野も混じってるけど、いくつかは、当時のアラビア科学でも時代遅れなものだった可能性もある。


この時代には、類似の分類は色々あって、Gundissalinusのも、その一つに過ぎない。別の例として、Hugh of Saint Victor(1096〜1141)がDidascalicon de studio legendiの中で述べてる分類では、数学は、以下のように、算術、音楽、幾何、天文学の4つの知識に分かれるみたいなことが書いてあって、Quadriviumと同一と思ってよさそう。

mathematica dividitur in quattuor scientias. prima est arithmetica, quae tractat de numero, id est, de quantitate discreta per se. secunda est musica, quae tractat de proportione, id est, de quantitate discreta ad aliquid. tertia est geometria, quae tractat de spatio, id est, de quantitate continua immobili. quarta est astronomia, quae tractat de motu, id est, de quantitate continua mobili.

ここでは、4つの分野が、それぞれ異なる種類のquantitate(quantity/量)を扱うとして特徴付けされている。al Farabiの分類も、基本的には、「何らかの量を扱うのが数学だ」という思想上にあると思われる。

12世紀に於いて、mathematicaは、「定量科学」くらいのニュアンスで認識されてたのだろう。現代では、数学と実験科学は遠い印象があるけど、scientia de aspectibusの研究者としてIbn al-Haithamは、光学実験をしてたようなので、アラビア科学では、"数学"だから実験しないということは、別になかったと思われる。そもそも、アラビア天文学が、観測と計算が合うことを重視した分野だし、中世ヨーロッパの天文学者も、しばしば天文観測機器を自分で作っている。

多分、中世ヨーロッパの大学に於いて、数学とはQuadriviumのことでしたというのは、そんなに間違ってないけど、大学の外では、もっと緩いイメージで捉えられてたかもしれない。とりあえず、以下で、Gundissalinusが「数学」に分類した分野を個別に概観していく。

但し、幾何学と算術は、現代の意味でも数学の一部と見なされてるので、特に触れない。また、音楽は、よく分からないので、やっぱりpass。

数理天文学

天文学は、Quadriviumの一つではあるけど、西暦1200年頃のヨーロッパは、数学の下地が何もない状態なので、数理天文学ができるようになるのは長い時間が必要だった。江戸時代の日本も、(商人の会計や財政管理はあったにしても)数学の下地がほぼない状態から、約250年かけて和算が発展し、平方根や代数方程式の数値解法、三角関数、初歩の微分積分(極値計算法と求積法)、対数と、それらの計算法を、測量や天文学、建築(規矩術)などに応用する術を習得した。

現代の数学教育でも、小学校から(殆どが1700年までに発見された結果しか含んでいない)高校数学を学習するのに約10年かけてることを考えると、知識が整理されてない時代に、知識源にアクセスするには言語の壁がある状態で、数百年単位の時間が必要だったとしても、そういうものなのかもしれないという気がする。


以下では、中世ヨーロッパにおける、数理天文学の受容状態を簡単に追跡する。

トレド翻訳学派に、クレモナのGerard(1114〜1187)という人がいる。多くのアラビア語文献を翻訳したが、1175年に、プトレマイオスアルマゲストラテン語翻訳を行い、写本しかない時代にしては、広く回覧されたと言われる。13世紀初頭には、パリ大学にも写本があったようだ。

SacroboscoのJohannes(1195〜1256)という人は、その生涯は曖昧だが、オクスフォード大学で学んだ後、1220年代初頭からパリ大学でQuadriviumを教えるようになったと伝えられる。1230年頃、De sphaera mundi("On the sphere the World"、天球論)という本を書いて、この本は、中世ヨーロッパの大学で、天文学の標準的な"教科書"となり、1472年には、フェラーラで印刷されたらしい。ガリレオの時代でも使われてたようだ。内容の詳細は知らないけど、プトレマイオスや中世アラビア天文学の初歩的な解説を与えていたらしい。三角関数のような内容は含んでなかったぽいので、この本で、十分な計算ができるようにはならなかったと思われる。

13〜14世紀のヨーロッパで、数理天文学をやってた人たちは、情報があんまりない。13世紀のカスティーリャ王国では、アルフォンソ天文表が作られており、1270年代に完成した。計算に関わったのは、トレド翻訳学派の学者Yehuda ben MosheやIsaac ben sidという人らしい。

Levi ben Gershon(1288〜1344)は、おそらくフランス出身のユダヤ系学者らしいが、生涯は殆ど何も分かってないそうだ。1342年にローマ教皇クレメンス6世の依頼で、Sefer Milhamot Ha-Shem("The Wars of the Lord")という"科学/哲学史"の本の一部をラテン語に翻訳したとか言われている。1342年に、De sinibus, chordis et arcubus("On sine, chord, and arcs")という三角法の本を書き、ヤコブの杖という天体測量機器の発明者とされることもある。

Georg von Peuerbach(1423〜1461)は、ウィーン大学(1365年創立)で学び、天文測量機器の作成や三角関数表の計算を行い、アルフォンソ天文表の修正を試み、そして1450年代に世界最初の天文学教授となった。彼は、パドヴァ大学ボローニャ大学フェラーラ大学などで教えたこともあるようだ。Peuerbachは唐突に出現したわけではなく、彼が師事したJohannes von Gmunden(1380~1442)は、ウィーン天文学派の創始者とされ、本格的に数理天文学を開始しようと努力していたようだ。

1514年に、Viri Mathematici quos inclytum Viennense gymnasium ordine celebres habuitという本が出版されていて、ウィーン天文学派の数学者の伝記が記載されているそうである。16世紀初頭には、天文学者が、mathematici(mathematician)と見られていたことが分かる。

Peuerbachの生徒に、レギオモンタヌス(1436~1476)がいて、1464年に 三角法の著書De Triangulis omnimodusを著し、1474年に天文表Ephemerisも出版した。レギオモンタヌスの著書は、丁度、ヨーロッパで商業出版が開始した時期だったこともあり、かなり広く読まれ、影響が大きかったと思われる。

レギオモンタヌスは、1475年、教皇シクストゥス4世から、ユリウス暦の改暦相談のため、ローマに呼ばれている。グレゴリオ暦の運用を開始したのが1582年だから、改暦に至るまで100年以上かかってる。ユリウス暦への批判自体は、もっと古くからあったようである。江戸時代の日本で、西洋天文学ベースの改暦を試みた際も、何度か失敗して長い時間がかかってるので、そういうものらしい。

カスティーリャ王国出身のAbraham Zacuto(1452〜1515)は、サラマンカ大学で、天文学を勉強し、1492年に、アルハンブラ勅令でスペインを追放された後、ポルトガルジョアン2世に仕えた。1496年には、ポルトガルでもユダヤ人追放令が出され、その後の足跡は、詳しくは分かってないっぽい。

Zacutoは、1478年頃、Almanach perpetuumという本をヘブライ語で完成させた。これには、太陽の黄道座標系での予測位置以外に、赤緯表も含まれていたっぽい。Zacutoが、この本を書いた理由は知らない。レギオモンタヌスが、同時期に書いて1490年に出版されたTabulae Directionum profectionumqueにも、赤緯表が含まれてるらしいが、どれが何の表なのか、よく分からない。黄道座標での位置を赤道座標に変換するのは、別に難しくない気がする(プトレマイオスの「アルマゲスト」でも出来てた話だと思う)けど、それ以上の何かがあったのか、それとも、当時は、座標を変換するのも簡単とは思われてなかったのか、不明。

ポルトガルジョアン2世が集めた天文学者/数学者たち(しばしば、Junta dos Mathematicos/数学委員会という呼称が使われるが、当時から、そのように呼ばれてたかは不明)は、Zacutoの表をベースにしつつ、1484年にRegimento do estrolabio e do quadranteという、簡単な数表を含む航海マニュアルにまとめ、Regimento das léguas、Regimento do norte、Regimento do Solと呼ばれる遠洋航海の手法が記載されているらしい(1914年版がオープンアクセスになってるけど、殆ど何もわからない)。もう少し後に、Regimento do Cruzeiro do Sulというのも追加され、これらの手法は、江戸時代の「元和航海書」(1618年)にも、天文表と共に解説されているっぽい。

詳細はともかく、15世紀末に、ジョアン2世の雇用した数学者たちが、新しい遠洋航海術を発明して、100年以上に渡って使われたことは確からしい。コロンブスは、アメリカ大陸到達時に、レギオモンタヌスの出版したEphemerisを携行していたという逸話が残っている(これは本当かどうかは分からないが、何らかの天文表は持ってたのだろう)。ヴァスコ・ダ・ガマも、何らかの天文表を使っていたと信じられているが、どの本かについては、多少の議論があるようだ(cf. Nautical Tables for Vasco da Gama, 1497–1500?)。

そんなわけで、15世紀には、ヨーロッパも、数理天文学を習得し、航海術などに応用することができるようになった。

scientia de aspectibus

ヨーロッパで光学の研究が盛んになるのは17世紀で、理論方面では、フェルマー、スネル、デカルトホイヘンスなどの光学研究は現在まで知られている。望遠鏡や顕微鏡のような光学機器も、この時期に作られた。ヨーロッパでは、13世紀末に老眼鏡が作られたと言われ、1600年頃のヨーロッパでは、近視/遠視用共に眼鏡は一般的なものになっていたらしい。ガラス職人のレンズ加工技術も高水準になってたんだろう。

ヨーロッパの光学研究先駆者に、13世紀のポーランド人Witelo(あるいはVitello)という人がいる。彼は、1260年頃、パドヴァ大学で学び、1270年代にPerspectivaという本を書いたそうで、この本は、Ibn al-HaithamのKitāb al-Manāẓir(The Book of Optics)に基づいて書かれたと考えられている。1572年に、Friedrich Risner(1533〜1580)が、Opticae thesaurusという本を書いて、そこにal-Haithamの本のラテン語訳も含まれてたようである。この本は、ヨーロッパで広く読まれ、光学研究ブームの切っ掛けとなったっぽい。Keplerは、1604年に、"Ad Vitelonem paralipomena”という本を書いた

記録上では、最古の望遠鏡は、オランダ人Hans Lipperheyが1608年に特許を申請したものとされる。Keplerの本は、それ以前のもので、内容的には、光学というより視覚論のようである。天体観測への利用を想定してたわけではないっぽい。尤も、ティコ・ブラーエは、Camera obscuraを利用してたし、大気の屈折率の高度依存性によって、星の視高度がずれること(大気差というわかりにくい名前が付いてる)にも気付いてて、補正を行ってたらしいので、視覚論への興味も、天体観測の問題から生じてた可能性はある(が、本当かどうかは分からない)。

13世紀に生きてたらしいJordanus de Nemore(1225?〜1260?、Jordanus Nemorarius)は、その生涯が殆ど何も分かってない人で、Liber de Speculisという光学の著書があるそうだが、未出版。

Jordanusは、Liber de ponderibusという"静力学"の著書で最も知られる。また、Jordanusは、De plana speraという本で、ステレオ射影に関する命題を幾何学的に証明したようだが、動機は、天文観測機器アストロラーベにあったらしい(cf.The Mathematics of the Astrolabe and Its History)。Jordanusは、算術なんかの教科書も書いてるようだが、扱ってるテーマは、Quadriviumを超えて、al FarabiやGundissalinusが"数学"とした領域に及んでいて、アラビア科学の"数学"から影響を受けたのでないかと思われるが証拠はない。

いずれにせよ、Jordanusの光学研究が、ヨーロッパに影響を与えたことはなく、中世の間、ヨーロッパに光学研究は殆どなかったように思える。おかげで、ヨーロッパの光学研究が、アラビア科学に起源を持つことは、かなりはっきりしている。

その後、光学の知識は、望遠鏡の改良のために、かなり長い間、天文学者の教養の一部のようなものとなった。

scientia de ponderibus

ponderibusは、pondusの複数形で、pondusはweightなどの意味らしい。ポンドは、重さの単位なのでわかりやすい。従って、scientia de ponderibus('ilm al-athqal)は、直訳では、science of weightsとなる。

The Arabic Transformation of Mechanics: The Birth of the Science of Weightsに、該当するアラビア語文献が32冊挙げられている。6冊は、古代のギリシャ語文献の翻訳で、1600年以降の文献を、少なくとも8冊含み、19世紀の本まである。ギリシャ語文献は、アリストテレスのMechanica、HeronのMechanica、PappusのMechanicaを含むので、scientia de ponderibusは、Μηχανικά/Mechanicaの直接の子孫と言っていいだろう。

Euclidの"on balance"と"on heaviness and lightness"という著作もあるらしい(原題は不明)。これらの本は、梃子の原理だったり滑車だったりを扱ってて、ユークリッドの他の著書と同様、公理的なスタイルで書かれてたようで、後者は、9つの公準と6つの定理が含まれていたとも書いてある。

Mechanicaは機械学と訳されることが多いので、scientia de ponderibusも"機械学"と訳すのが妥当っぽい。ただ、機械と言っても、梃子、天秤、滑車とかなので、現在イメージする機械とは、ちょっとずれてる。mechanicaは、英語では、mechanicsになるけど、classical mechanicsやfluid mechanicsのように、力学の意味で使う今と違って、昔は、機械学の意味で使われてる。現在、機械工学は、英語だとmechanical engineeringなので、mechanicsから派生したんだと推測できる。


mechanicsの使用例として、ガリレオ(1564〜1642)が書いたRaccolta di Quelle cognizioni che à perfetto Cavaliero and Soldato si richieggono, the quali hanno dependenza dalle scienze matematiche ("完璧な騎士と兵士に求められる数理科学の知識")という文書を見る。これは、1608年、パドヴァに作られたAcadémie Deliaの教育カリキュラムに関するガリレオの提案書らしい。この文書では、算術、幾何、体積計測、機械学(scienze mecaniche)、コンパスや計測機器の使い方に加えて、砲術、軍事建築/築城術(Architettura militare)などが挙げられている。

ここのscienze mecanicheは、後に続く文章を見ると、

Cognizione delle scienze mecaniche, non solo intorno alle loro ragioni e fondamenti communi, quanto intorno a molte machine ed instrumenti particolari, insieme con la resoluzione di moltissime questioni e problemi da essa cognizione mecanica dependenti.

のように、"machine ed instrumenti"(machines and instruments)というフレーズが出てきて、様々な機械や機器に精通して、トラブルに対処できなければならないという感じのことを書いているので、"機械学"と訳すべきだろう。

ついでに、ガリレオがArchitettura militareと書いてるのは、16~19世紀に存在した築城術(fortification)という分野を指していると思われる。15~19世紀には、ヨーロッパを中心として、銃火器に対抗するために生まれた星形要塞というのが流行っていて、日本では、函館五稜郭に見られる。最古の星形要塞は、イスタンブールにあるYedikule要塞(1458年建設)とされる。

築城術は星形要塞の外形形状設計術みたいな分野っぽい。少なくとも、中世の時点では、土木工学的要素は殆どないけど、築城術は、数学の応用と見なされていたらしく、1736年の本The Young Trigonometer's Compleat Guideにも、築城術への応用が書かれている。この本は、適当に検索して見つけた古い本で、別に有名とかではないと思う。ガリレオ自身も、fortificationに関する原稿を書いている(例えばTrattato di fortificazione)。

実際のところ、中世の築城術自体には、(当時の水準で言っても)数学や物理として高度なことは何もないように見える。ただ、ウィトルウィウスのDe architectura(BC30頃?)には、幾何学日時計天文学、機械に関する記述があって、建築を監督する者は幅広い数学的素養を持つべきという思想が存在したそうである。築城術を数学の一分野のように見なす発想も、そういうところから来てるのかもしれない。


また、mechanicsの17世紀英語の文例として、John Wallis(1616〜1703)が、1640年代のロンドンにあった非公式の集まりについて書いてる文章には

Our business was ( precluding matters of theology and state affairs ) to discourse and consider of Philosophical Enquiries , and such as related thereunto ; as Physick , Anatomy , Geometry , Astronomy , Navigation , Staticks , Magnetics , Chymicks , Mechanicks , and Natural Experiments ;

とある(出典:The Royal Society 1660-1940)ようで、Staticksが静力学に対して、Mechanicksは機械学だろう。

scientia de aspectibusと違って、アラビア科学のscientia de ponderibusのヨーロッパへの影響は、はっきりしてないが、少なくとも、"mechanics"とscientia de ponderibusは、どっちも、古代ギリシャ/オリエントのMechanicaに由来する同一ルーツの分野と思われる。そして、1600年頃のヨーロッパで、mechanicsを数理科学の一つに位置付ける場合があり、Gundissalinusも、scientia de ponderibusを数学の一分野としていた。


mechanicsのニュアンスが変化した要因は、EulerのMechanica sive motus scientia analytice exposita (1736)が大きいかもしれない。その前後の力学の本を見てみると、Jakob Hermannは、 Phoronomia, sive De viribus et motibus corporum solidorum et fluidorum libri duo (1716)を書き、d'Alembertは、Traité de dynamique (1743)を書いている。Lagrangeは、Mécanique analytique (1788)を書いた。phoronomyは、kinematicsと同じ意味だそうで、kinematicsは、運動を直接扱い、運動の原因(古典力学なら力や慣性)は扱わないというニュアンスがあるとされている。

多分、dynamicsとmechanicsの差は、あんまりない。流体力学は、昔はhydrodynamicsだった(Bernoulliの1738年の本が最初で、Lambの1895年の本のタイトルにも見られる)けど、最近は、fluid mechanicsが一般的な気がする。量子力学は、quantum mechanicsなのに、量子電磁力学や量子色力学が、quantum electrodynamicsやquantum chromodynamicsなのは、どういう経緯があったのか謎だけど、後者の2つは、相互作用の原因に迫ってるから(?)。古典力学だと、重力やバネの力、摩擦力なんかを、現象論的に導入するが、その正体には関知しないので、dynamicsよりmechanicsの方が適切というニュアンスが現在はあるかもしれない。だとしても、この使い分けも、18世紀や19世紀に存在したものではないと思う。

scientia de ingeniis

scientia de ingeniis('ilm al-hiyal)は、該当文献が不明で、実態が一番よく分からない。

Jordanusのアストロラーベとステレオ射影に関する本De plana speraを、scientia de ingeniisに分類してる人もいる。Jordanusは、そう考えてたのかもしれないが、アラビア語文献ではないので、al Farabiの考えとは違うかもしれない。

イスラム法学では、hiyalとは、"本来なら違法にしか達成できない目的を法律的な工夫で合法に達成する手段"を指すらしい。現代だと節税などが例になるそうだ。おそらく、原義は"工夫"とか、その程度の意味だったと思われるが、そのhiyalと同じ単語っぽい(?)。また、ラテン語ingeniisは、ingenium(辞書によると、machine,deviceの意味があったらしい。engineと同根ぽい)の複数形で、ingeniumに"人"を意味する接尾辞torと複合すると、ingeniatorになる。ingeniatorは、英語のengineerなので、scientia de ingeniisは、engineering scienceつまり工学としてもいいだろう。


単語としてはそうでも、実際の内容は違う可能性もある。NOTES ON MECHANICS AND MATHEMATICS IN TORRICELLI AS PHYSICS MATHEMATICS RELATIONSHIPS IN THE HISTORY OF SCIENCEには、"The science of devices refers to applications of mathematics to practical use and to machine construction"と書いてある。

タイトルにingeniisが含まれるラテン語写本には、

  • De ingeniis spiritualibus
  • De decem ingeniis seu indicationibus curandorum morborum (1300年頃)
  • De aquarum conductibus et ingeniis erigendis

などがあって、1つ目は、ビザンチウムフィロンによる著作Mechanike syntaxisのPneumaticaの一部の翻訳らしい。pneumaticsは空気圧工学とか訳されるらしいが、当時は流体で駆動する"機械"(水車やポンプ、水オルガンなど)全般を指してて、現代の水理学に相当する内容も含んでた可能性もある。2つ目は、何か"医療機器"(?)の話っぽいけど、中世のラテン語に確信が持てない。3つ目は、アレクサンドリアのHeronのPneumatica(アイオロスの球という蒸気機関の一種の記載があるらしい本)のラテン語訳っぽい。訳者はトレド翻訳学派のWillem van Moerbekeらしい。

タイトルにhiyalが含まれるアラビア語文献としては、

  • Kitáb al-Ḥiyal(一般に"The Book of Ingenious devices"と訳される)
  • Kitáb fi ma'rifat al-hiyal al-handasiyya(一般に"The Book of Knowledge of Ingenious Mechanical Devices"と訳される)

があり、前者は、9〜10世紀のバグダッドにいたBanū Mūsā兄弟の著書。当時知られていた多くの発明を記載しているそうである。彼らは、数学や天文学の著書も残してるらしい。後者は、12世紀のアラブ人Al-Jazariによる著書で、カムシャフトやクランクシャフトを記載した最初の文献とされる。この中には、水力機械の記述もあるらしいので、pneumaticaの影響もあったのかもしれない。

本のタイトルを見る限り、このへんの分野が、scientia de ingeniisという気もするけど不明。


オスマン帝国の数学者/天文学者Taqi al-din(1526〜1585)が書いた機械式時計の本が2冊あって、一冊は、al-Kawākib al-durriyya fī waḍ ҁ al-bankāmāt al-dawriyya(1556or1559?)で、もう一冊は成立年不詳のKitab al-Ṭuruq al-saniyya fī al-ālāt al-rūḥāniyyaである。当時のオスマン帝国は商業出版を行ってない(活版印刷イスラムの戒律に反するという主張があったとか言われている)ので写本だと思うけど、後者はオープンアクセスになってるのを見つけた。私には微塵も読めないので、本物かどうかすら分からない。

Wikipediaによれば、後者には、Banū Mūsā兄弟とAl-Jazariへの言及があり、時計のgeometrical-mechanical structureを強調したと書いてある。つまり、Taqi al-dinは、自身の機械式時計の研究開発、Banū Mūsā兄弟の本、Al-Jazariの仕事を同分野と位置付け、数学的なものと見なしていた可能性が高い。Taqi al-dinの時計は、最小単位が5秒だったそうで、天文学者でもあったTaqi al-dinは、機械式時計を天文観測に利用した。

イスラム世界では、13世紀末に、muwaqqitという礼拝の時間を管理する仕事が誕生し、天文学の知識を持つ者が選ばれたそうである。それに関連して、ʿilm al-mīqāt(science of timekeeping)という分野が発生したらしい。どういう内容なのか知らないけど、初期の研究者として、Shams al-Din al-Khalili(1320〜1380)とIbn al-Shatir(1304〜1375)という人がいる。ヨーロッパでも、時計学/horologyという分野名が、かつては存在していたけど、ʿilm al-mīqātとの関連性は不明。Taqi al-dinの時計研究には、そういう方向からの影響もあるのかもしれないが、これも不明。

仮に、水車、オートマトン、機械式時計などの研究が、scientia de ingeniisだとすると、こっちも"機械学"に見える。個人的な感覚では、scientia de ponderibusとscientia de ingeniisの違いは、"機械"の複雑さという感じがする。scientia de ponderibusで扱う"機械"(天秤、滑車、梃子など)は、原始的な機械要素が多く、理論的解析はしやすいが、機械の印象が薄い(レバーとか投石機になると、機械感がなくもないけど)。scientia de ponderibusがmechanicsだとすれば、scientia de ingeniisは、mechanical engineering/機械工学に近い分野かもしれない。artes mechanicae(織布、建築、軍事・狩猟、鍛冶・冶金などの技術で、リベラル・アーツと対比される)との区別が曖昧でもある。


Taqi al-dinの一冊目の本では、ヨーロッパの機械式時計コレクションを調査した旨が述べられているそうなので、16世紀時点でヨーロッパの機械技術は高かったと推測される。イギリス人Richard of Wallingford(1292~1336)という人が書いたTractatus Horologii Astronomici (1327)は、機械式時計に関する(多分、世界で)最初期の文献らしい。彼は、オクスフォード大学で学んだ修道僧だったそうだ。

中世ヨーロッパでは、多くの天文時計が作られたが、1579年に、Jost Bürgi(1558~1632)が秒針を持つ時計を作ったとされる。Bürgiが、どこで教育を受けたのか、全く分かってないらしい。Bürgiは、Napierと独立に対数を発見したとも言われ、後にKeplerの計算助手になった。また、2015年になって、三角関数を計算するBürgiのアルゴリズムの詳細と証明を書いた論文が出ている(Jost Bürgi's Method for Calculating Sines (arXiv:1510.03180))。

オランダのホイヘンス(1629〜1695)は、機械式時計に対して理論的なアプローチを行い成功した。ホイヘンスは、1673年に、Horologium oscillatoriumという本を出している。タイトルは、振り子時計を意味する。ホイヘンス以前は、日差15分程度くらいの誤差があったらしいので、大体、1%くらいのずれがあったようだ(cf.機械式時計の歴史 1)

ついでに、当時の時計の分解能では、重力加速度を、自由落下時間によって高精度に決定するのは難しかったので、振り子の周期と長さから間接的に重力加速度を決定することが試みられた。最初に、この測定をしたのは、フランスのメルセンヌっぽいけど、単振り子では、振れ角に依存して周期が変わることも実験的に発見していたようだ。ホイヘンスは、サイクロイド振り子を使って、周期が二秒になる振り子(seconds pendulumと呼ばれる)の長さを測定し、この本に数値が書いてあるらしい(しかし、どこに書いてあるのか発見できなかった。ラテン語力がほしい)。

現代の単位では、この長さは、約1メートルで、18世紀末フランスのメートル法制定時に、これを基準値とする案は強く推された(cf.Why g〜π2Why does the meter beat the second? (arXiv:physics/0412078)の特に§3と4)。また、ニュートンは、月の加速度と、(ホイヘンスの測定値から算出した)重力加速度の比率が、地球半径と地球・月間距離の比の二乗に近いことを確認した。

補足)逆二乗則の推測は、等速円運動のkinematicsのみでもできる(現在の知識だと、等速円運動の軌道を座標で書いて二回微分すると、遠心加速度が分かる)。地球と月の距離が、地球半径の約60倍であることは、プトレマイオスの時代には視差の観測によって知られていた。地球半径も、紀元前から多くの推定があったが、フランスのJean Picardの1669年の測量結果では、(現在の単位に換算して約)6372(km)だった。月の公転周期は約27.3日なので、仮に完全な等速円運動をしていたとすると、遠心加速度は、(60 * 6372e3) * (2pi / (27.386400))**2=0.0027129(m/s2)と見積もられ、9.8(m/s2)は、その3612倍になる。これは、60の二乗に近い。この簡単な試算から逆二乗則に到達した人は複数いたかもしれない。また、等速円運動の時は、Keplerの第三法則から逆二乗則を導くことも容易(角速度は、公転周期に反比例し、遠心加速度は、軌道半径と角速度の二乗の積だから、Keplerの第三法則が正しければ遠心加速度は軌道半径の二乗に反比例する)。ニュートンは、この結果を楕円軌道の場合に拡張していた。重力加速度と地球半径の高精度の測定は、惑星運動と自由落下が同一の原因で生じている可能性を検証できるようにした。大雑把に言って、ここまでが1700年頃の到達点だけど、等速円運動近似に基づく試算は、十分とは言えず、もっと精密なものが求められただろう。けど、精密な検証をしようとすると、月の軌道が、地球と太陽の両方から無視できない影響を受けてるせいで、難易度が一気にあがる。この困難は、ニュートンオイラーも解決できなかった。

ホイヘンスの本と同時期の1674年に、デンマークの「数学者」Ole Roemerは、エピサイクロイド歯車に関する研究を発表した(文献未確認)。これ以降、歯車(の歯形曲線)を研究した「数学者」は多く、その中には、オイラーも含まれる。


16世紀〜17世紀前半の時点で、ヨーロッパに、機械に詳しい「数学者」が結構いたことは確かだけど、彼らが複雑な機械も数学的に解析できると思ってたのかは分からない(統一見解があったかも疑わしい)。当時は、複数分野に手を出す人が沢山いたから、単に機械にも詳しい「数学者」だったかもしれない。例えば、Gerolamo Cardanoは、三次・四次方程式の研究で知られる「数学者」だけど、本職は医者で、Cardan joint(Hooke's jointなど、色んな名前がある)というリンク機構に名前を残していて、錬金術に関する著書もある。

1695年に出版されたフランス人LahireによるTraité de mecaniqueという本の4ページには、"nouvelles machines qui ont été inventées par de tres-sçavans Mathematiciens,& par de tres-habiles ouvriers"(有能な数学者と職人によって発明された新しい機械)というフレーズがあって、17世紀末には、数学者を、機械の発明にも関与する人々と見なす場合があったっぽい。これはフランス語の本で、ヨーロッパ内でも、国によって認識の差異があった可能性は否定できない。


江戸時代の日本では細川頼直(??~1796)という人が、天文学と暦学を学んだ「天文学者」であり、『機巧図彙』(1796)という機械の本を書いている。この中には、時計や、お茶くみ人形のようなオートマトンが記述されているらしい。江戸時代のからくりも、機械式時計で使われていた機構を応用して発展したととか検索すると書いてある。第四章 実学としての和算|江戸の数学には、竹内度道(1780〜1840)という和算家が水車を設計し大工を指導してたとか書いてある。機械に関わる和算家というのが、数学者が機械に強いというヨーロッパの風潮が日本にも伝播しただけなのか、日本で独自に発展したものなのかは分からないけど、どっちだとしても、面白い事象ではある。

scientia de ingeniisは不明点が多いけど、結局は、al Farabiが10世紀に"数学"だと認識してたものを、ヨーロッパは、分裂した形で受容して、バラバラに発展させた後、再統合したというだけの話かもしれない。そして、19世紀には、また分裂した。

おまけ:16世紀の力学研究者

16世紀のmechanicsは、天文学と無縁の分野だったので、大学のカリキュラムと関係ない。そのせいか、16世紀の力学研究者は、その研究によって大学の職を得ていたわけではないし、独学の人も多い(従って「数学者」かどうかの定義が曖昧になる)。中世ヨーロッパの天文学者が、大体、大学で天文学を勉強してるのと比べると、対照的。

16世紀の力学研究者として、William Whewellの1840年の本History of the inductive sciencesに載ってる人として、タルタリア(1500?~1557)、Gerolamo Cardano(1501~1576)、Giambattista Benedetti(1530〜1590)、Michael Varro(1542〜1586)、Guidobaldo del Monte(1545~1607)、Simon Stevin(1548〜1620)、ガリレオ(1564~1642)などがいる。多分、何らかの著作が残ってる人を挙げていると思う。

ガリレオを除けば、大学で"数学"を教えた人は、いないと思う。こういう"野良数学者"が、昔からいて商業出版によって可視化されるようになったのか、この時期に実際に増えたのかは分からない。似た現象として、中世ヨーロッパの有名な魔術師や錬金術師も、この時期の人が多い。WikipediaEuropean alchemistsというリストがあるけど、誕生年で分類すると、12世紀:2人、13世紀:5~7人,14世紀:2人,15世紀:9人,16世紀:27人,17世紀:21人、18世紀:6人となっている。

ガリレオは、大学を中退して、(トスカーナ大公2代目に仕えた)宮廷数学者のOstilio Ricci(この人はタルタルアの弟子)に数学を習ったとされる。1586年に書いた論文は、La Billancetta("little balance"/小天秤)で、アルキメデスの原理の実践的利用法みたいなお話。

これが、大学で数学を教える仕事に結び付くのか不思議だけど、彼は推薦状を書いてもらうために、偉い人に面会して、イエズス会士Christopher Clavius、Guidobaldo、その伝手でトスカーナ大公(3代目)などのコネを得たらしい。最終的に、トスカーナ大公の推薦でピサ大学で教えることになったが、当時、そんな求職方法が、一般的だったわけではないだろう。

Claviusは、イエズス会の中心的教育機関ローマ学院で、幾何学天文学を教えてた学者で権威とされていた。詳細は知らないけど、イエズス会数学教育を改革したとか言われている(イエズス会学事規定1599年版)。ガリレオが異端審問にかけられた時、Claviusは既に死んでたけど、生きてたら、宗教裁判を回避できたかは定かでない。Claviusには、Opera mathematicaという著書があって、読んだわけではないけど、内容は、幾何、算術と代数、天文学の話題で構成されているっぽい。アストロラーベに関する著書Astrolabiumは、漢訳『渾蓋通憲図説』が存在するという事実から、よく知られている。

Guidobaldoは、何者かよく分からないけど、Wikipediaによれば、裕福な家の出身で、1564年のBattle of Saint Gotthardに参加した後、Federico Commandino(1509~1575、ヘレニズム時代のギリシャ語文献を何冊もラテン語翻訳した)に師事して、数学を学んだらしい。1588年からは、トスカーナ大公国の"築城監督官"の役職を得ていた。ガリレオがピサ大学をクビになった時も、Guidobaldoの援助があって、パドヴァ大学の職を得たらしい。

そんなわけで、ガリレオも野良数学者みたいな人に師事して、苦労して、大学で仕事を得たっぽい。

Simon Stevinの力学上の仕事は、1586年に出版したDe beghinselen der weeghconst(The principles of weighing)とDe Beghinselen des Waterwichts(The principles of water balance)の2冊の本。Stevinも、独学だったそうだが、やはり誰に何を学んだかは分かっていない。

A first mathematics curriculum : Stevin’s Instruction for engineers (1600)によると、1600年、Simon Stevinは、ライデン大学内で、"Duytsche Mathematique"というコースを開くよう要請されたそうだ。

名前は、数学だけど、目的は、military engineerの教育。カリキュラムは、fortificationがメインのようで、それに必要な算術や幾何学、測量も含まれる。機械の使い方なんかは、記述がないっぽい(?)し、砲術も含まれてない。建築機械については、実地で学ぶ感じだったのかもしれない。デカルトやOtto von Guericke(1602~1686)も、ここで勉強したことがあるらしい。

おまけ2:中世以降の砲工兵教育

17世紀にはヨーロッパ諸国で、砲兵隊、工兵隊が専門分化し、17〜18世紀には、西ユーラシア諸国で、砲工兵学校が作られた。17世紀にも学校は存在したが、正式名称も定かでないものが多い。18世紀のものを、いくつか挙げておく。

  • 1720年のスペインで、Real Academia Militar de Matemáticas y Fortificación(直訳は、Royal military academy of mathematics and fortification)
  • 1739年にサルデーニャ王国トリノで、Regie Scuole di Artiglieria e Fortificazione(直訳では、Royal School of artillery and fortification)
  • 1741年のイギリスで、WoolwichのRoyal Military Academy(王立陸軍士官学校)
  • 1748年のフランスで、École royale du génie de Mézières(メジエール工兵学校)
  • 1773年のオスマン帝国で、Mühendishâne-i Bahrî-i Hümâyun(Imperial Naval Engineering School)
  • 1795年のオスマン帝国で、Mühendishâne-i Berrî-i Hümâyun(Imperial School of Military Engineering)

Regie Scuole di Artiglieria e Fortificazioneは、若い頃のLagrangeが教師をやってたことがあり、メジエール工兵学校では、Mongeが教えていたことがある。Mongeは、Monge-Ampere方程式に名前が残ってるけど、エコール・ポリテクニークの教育カリキュラム決定に大きく関与した人でもある。

The Teaching of Mathematics at the Royal Military Academy: Evolution in Continuityに、1792年のWoolwich Royal military academyのカリキュラムが列挙されているが、初等的な内容以外に、力学(mechanics)、流率法(fluxions,微積分)、水理学(hydrostatics and hydraulics)、pneumatics、砲術(Gunnery)などが含まれている。

ここのmechanicsは、力学とする方が適切だと思う。詳細は、以下のように書いてある。

Mechanics: Including Motions equable and variable; Forces constant, variable and percussive; Gravity, Sound and Distances; Inclined Planes; Projectiles; Pratical Gunnery; Pendulums; Centres of Gravity; Percussion, Oscillation, and Gyration; Ballistic Pendulum, &c.;

出典は、1892年のRecords of the Royal Military Academy, 1741-1892という本らしいけど、1792年の文言を、どれくらい正確に収録してるのかは知らない。機械工学の方は、1846年には、イギリスで、スチーブンソンが発起人となってInstitution of Mechanical Engineers(英国機械学会)が作られているので、19世紀には、mechanicsとmechanical engineeringが区別された。

18世紀当時は、水理学や弾道学も、「数学者」の研究対象だった。Johann Bernoulliは、1743年に"Hydraulics"という本を出している。弾道学も、例えば、Eulerが研究を行っている。Eulerは、論文を書く以外に、イギリスの砲工兵Benjamin Robinsによる弾道学の本を翻訳して、詳細な注釈を付けている。BernoulliやEulerが関与していることから推測される通り、これらの分野は、当時としては、かなり数学的なものではある。

‘Mixed’ mathematics in engineering education in Spain: Pedro Lucuce's course at the Barcelona Royal Military Academy of Mathematics in the eighteenth centuryは、スペインのRoyal military academyのカリキュラムについて書いている。おそらく18世紀半ば頃のもので、Curso Mathematico para la Instruccio ́n de los Militares(軍事教育の数学コース)は、大きく8つの主題から構成され、算術、ユークリッド幾何、実用幾何(平面三角法や対数、測量など)、築城術、砲術、cosmography(celestial sphere/天球と書いてるのは、天文学?他に、地理学や海図作成など)、静力学(staticsと書いてるが、内容は、"the motion of heavy bodies"、機械学/machinery、水理学、光学など)、土木建築になっている。各主題は、概ね同程度の分量だったっぽい。まだ微積分はないけど、オイラーのMechanicaが出て10年足らずの時期のカリキュラムなので、まだ十分に、その有用性が知れ渡ってはいなかったのだろう。

mixed mathematicsは、応用数学の古い呼び方で、1800年頃までは、よく使われたらしい。The Evolution of the Term "Mixed Mathematics" には、mixed mathematicsという用語を最初に使ったのは、フランシス・ベーコンだと書いてある。

18世紀は、理論物理学が形成され、そのための数学的道具として、微分方程式論が開発された時代でもある。 一方、d'Alembertは、幾何、算術、代数、微積分学を純粋数学に分類し、力学、天文学、光学、確率論をmixed mathematicsに分類したそうだ。これは、確率論を除けば、現代の認識に近い。確率論は、現在は、数学の一分野と思われてるけど、ヒルベルトの第6問題「物理学の公理化」で、物理学と確率論の公理化をあげてるので、1900年頃でも、微妙な立ち位置にいたのだろう。

経験から完全に切り離して、機械的な論理操作のみで数学ができるという発想が、いつから存在してたのか分からない。少なくとも、算術の公理は19世紀末になるまで存在してない。純粋数学は、日常感覚に任せてた部分を整理する方向に進み、mixed mathematicsに分類されてた分野は、自然科学と融合して、物理学になった。この仕切り直しが、本当に良かったのかは分からないけど、とりあえず、名前を変えてくれなかったせいで、昔の話をする時は、ややこしくなった。

遺伝子の起源

生物が持つ遺伝子は、翻訳されて、数百〜数千アミノ酸のタンパク質になるのが一般的である。数百〜数千残基の長大なタンパク質or遺伝子が、ランダムにアミノ酸ヌクレオチドの重合を繰り返しただけで出現する確率は、普通に考えると極めて低い。ヒトの遺伝子で最も巨大なものの一つはtitin遺伝子で、363個のエクソンで構成され、最大のisoformは36000残基近くになるようである。

一般的に、生物の持つタンパク質は、複数の共通するドメインからなることが知られていて、ヒトの遺伝子数は2〜3万とされている(選択的スプライシングがあるので翻訳されるタンパク質の種類は10万種類くらいになると考えられている)が、1000個程度のドメインの組み合わせでバリエーションを増やしているらしい。また、真核生物、真正細菌古細菌間でも、共通するドメインの利用が確認されているとされる。

cf)Evidence for PDZ domains in bacteria, yeast, and plants.
https://dx.doi.org/10.1002%2Fpro.5560060225

SMART: Domain List
http://smart.embl-heidelberg.de/smart/domain_table.cgi
を見ると、1400以上のドメインが列挙されている。

多くのドメインは、数十〜数百残基程度のアミノ酸からなり(例として、PDZドメインは、90残基弱で、ホメオドメインは60残基程度、SAMドメインは70残基程度、SH2ドメインは100残基程度、pairedドメインは100残基以上、MAMドメインは170残基程度...)、固有の立体構造を取り、単体でも何らかの機能を有し、ドメインの組み合わせによって、新規遺伝子が作られることもあったと考えられている。が、これでも、まだちょっと長い。

例えば、1Lの水中で、毎秒1mM相当の新規タンパクが生成する(つまり、毎秒6.02×10^20個の新規タンパクが生成される)ような(法外な)イベントがあったとしても、一億年かけても、最大で10^36個オーダーのタンパク質しか試せない。現在の海水(その量は、10^21リットルのオーダーだそうだ)を全部使っても、20の100乗よりは大分小さく、ランダムな重合によって、100残基のタンパクを全部網羅するには到底及ばない。

というわけで、生物が遺伝子ガチャに成功している理由が説明されるべきで、以下のように、いくつかの可能性が考えられる。
①全ての遺伝子は、単一の遺伝子に由来する
②ランダムに生成したアミノ酸配列でも高確率で何らかの活性を持っている:例えば、ドメインの配列も高度に保存された(連続してるとは限らない)短い配列と、それ以外の部分からなることが多々あり、保存された配列以外の部分は、割と何でもいいとか
ドメインより小さい構造単位が存在する


現在、多くの遺伝子が、遺伝子重複によって、既存の遺伝子を雛形にして生成し、重複遺伝子は元になった遺伝子と独立な変異を経ることで、しばしば機能を変化させてきたと考えられている。ヒトに多くの相同遺伝子(パラログ)が含まれていることが分かっているし、大腸菌酵母も、パラログ遺伝子を持っていて、これらの事実は遺伝子重複によって説明できる。

1970年に出版された大野進の本Evolution by gene duplicationには、“In a strict sense, nothing in evolution is created de novo. Each new gene must have arisen from an already existing gene.”と書かれていて、これは①の考え方に相当する。当時は知られてなかった遺伝子の水平伝播も、今ではありふれた現象と考えられてるけど、既存の遺伝子を雛形として新しい遺伝子を獲得するという点では共通している。

進化の歴史上、このような方法によってのみ遺伝子の獲得が起きたとすれば、確率的な困難は、最初の遺伝子が誕生した時に一回だけ存在したことになり、一回だけなら、奇跡が起きることもあるだろう。但し、このような奇跡が必要だったなら、地球以外の星に生命体が存在する可能性は極めて低いことになる。


遺伝子重複は、古くから予想されていたことではあるけど、最近は、個別の遺伝子について、かなり詳細な進化の軌跡を辿ることができるようになっている。

遺伝子重複の顕著な例は、ヒトのLオプシン(赤オプシン)とMオプシン(緑オプシン)で、この2つの遺伝子はX染色体上に並んで存在する。SオプシンはX染色体上になく、分岐したのはもっと古いと思われる(塩基配列上も、最も違いが大きい)。現在、原猿や有胎盤類は、SオプシンとL/Mオプシンの2色性色覚のものが多いらしく、真猿類の中だと、新世界ザルのオスは(一部を除いて)2色型色覚で、旧世界ザル、類人猿、ヒトは3色型色覚を持つ。新世界ザルは、3〜4000年前に真猿類の一部がユーラシア大陸から南米に移住して分岐し、アフリカに残った真猿類の一種で、L/Mオプシンの遺伝子重複が起き、この種は、旧世界ザルと類人猿の共通祖先になったと考えられている。新世界ザルの祖先が、海を渡った方法は謎だけど、当時は今より大陸間は近かったと考えられている。

と、よく書いてあるのだけど、ヒトの場合、X染色体上にオプシンを3つ持つ人がいるらしい。
OPN1LW
https://www.ncbi.nlm.nih.gov/nuccore/NM_020061
OPN1MW
https://www.ncbi.nlm.nih.gov/nuccore/NM_000513
OPN1MW2
https://www.ncbi.nlm.nih.gov/nuccore/NM_001048181
OPN1MW3
https://www.ncbi.nlm.nih.gov/nuccore/NM_001330067

OPN1MW2とOPN1MW3は、目で見ると同じもののようなので、ヒトのX染色体上には、OPN1LWとOPN1MW、OPN1MW2の3つのオプシンがある。OPN1MW2は、ある人とない人がいるらしい(?)から、かなり最近になって出現した遺伝子と考えられる(数千年とか数万年とか)。OPN1MWとOPN1MW2は、アミノ酸配列が一致してるが、一部の人は、これによって4色型色覚となるようなので、多分、なんか多型があるのだろう。

チンパンジーには、(私が探した限り)2つしかないっぽい。ヒトでもチンパンジーでもXq28領域に遺伝子がある。
OPN1LW領域(チンパンジー)
https://www.ncbi.nlm.nih.gov/nuccore/CU467112.1?report=fasta&from=144100&to=156324
OPN1MW領域(チンパンジー)
https://www.ncbi.nlm.nih.gov/nuccore/CU467112.1?report=fasta&from=186659&to=199065
ヒトXq28領域
https://www.ncbi.nlm.nih.gov/nuccore/AC092402.2

どーでもいいけど、
チンパンジーX染色体
https://www.ncbi.nlm.nih.gov/nuccore/1348310337
ゴリラX染色体
https://www.ncbi.nlm.nih.gov/nuccore/NC_044625.1
を見ても、OPN1LW/OPN1MWに相同な遺伝子が一個しか見つからないので、アセンブルをミスってるような気がする。

ヒトのOPN1LWとOPN1MWは、共に364aaで、アミノ酸13個の違いしかない。通常、オプシンは、一つの視細胞に一種類しか発現してないらしいけど、4色型色覚の人で、そのへんの発現制御がどうなってるのかは謎。

OPN1MWとOPN1LWの場合は、上流に共通のエンハンサー領域LCRというのがあって、一方の遺伝子のプロモーターを(ランダムに)選択し、一旦行った選択は変化しないらしい(どうやって安定化するのかは分かってないらしい)。ゼブラフィッシュでは、RH2-LCRという領域が4つの下流にあるRH2遺伝子の発現制御に関与しているそうなので、ヒトでも、OPN1MW2が追加されたところで、同種の機構が働いているのかもしれないけど不明要素が多い

OPN1LW,OPN1MW,OPN1MW2とOPN1SWは、エクソンの数が合致していて、個々のエクソンのサイズもよく似ているから、これらの遺伝子が単一の遺伝子に由来するというのは尤もらしい。ヒトのオプシンファミリー遺伝子には、他にも、OPN2/RHOやメラノプシン(OPN4)などの複数の遺伝子があるが、これらではエクソン数が違っている。

OPN1以外のオプシンファミリー遺伝子も、元を辿れば、単一の遺伝子に起源があると信じられている。更に、オプシンは、GPCR(Gタンパク共役型受容体)で、最初のオプシンは、他のGPCR遺伝子から遺伝子重複と変異による新機能獲得(neofunctionalization)で生成された遺伝子かもしれない。ヒトには1000近いGPCR遺伝子があるが、277個のGPCRの系統樹を作成した論文(PMID 12429062)によれば、オプシンファミリーに最も近いのはメラトニン受容体のようである。カイメンには、メラトニン受容体もオプシン遺伝子もないらしいが、また別の論文(PMID23112152)には、平板動物はオプシンを持つと書いてある。

原核生物古細菌も、バクテリオロドプシンやSchizorhodopsin、Heliorhodopsinと呼ばれる遺伝子を持つが、これらはGPCRではないようだ。レチナールと結合し発色団としている点は共通し、また7回膜貫通型受容体ではあるらしいから、全てのGPCRの起源という可能性はなくもない(?)(7回膜貫通型受容体という用語は、GPCRとほぼ同義で使われる)


他の例として、ホタルのルシフェラーゼ遺伝子は、遺伝子重複とneofunctionalizationで出現し、ついでに、その後、更に遺伝子重複が起きたと考えられるらしい。
cf)発光生物学 2つの最新トピックス
https://doi.org/10.1271/kagakutoseibutsu.57.409


単一遺伝子の重複だけでなく、ゲノム単位の重複も起きると考えられている。Hox gene clusterは、ショウジョウバエでは、8遺伝子からなる単一のクラスターだけど、多くの脊椎動物では、4つのクラスターがあり、それぞれのクラスターは10個前後の遺伝子からなる(進化の中で、一部が脱落したり増えたりしてるけど、並び順は保存されてる)。ヒトやマウスではHoxa,Hoxb,Hoxc,Hoxdと名前がついてて、異なる染色体上(ヒトではそれぞれ7,17,12,2番染色体)にあるけど、これはゲノム単位の重複が2回起きて生じたという仮説がある。クラスター単位で重複して転座しただけという可能性もあるけど、ParaHoxクラスターとかも4つあるようだ。

他にも、哺乳類のFOXP遺伝子は4つあるけど、ハエには一つしかない。2回ゲノム重複があったと考えると、丁度いいケースが結構ある。2010年代後半までの報告では、脊椎動物に近縁のナメクジウオやホヤが単一のHoxクラスターしか持たない一方、無顎類の生き残りである円口類のHox gene clusterやParaHox gene clusterは複数あるものの、どっちも2つとか、両方4つずつというわけではなく、ゲノム重複の起きた時期は、まだ完全には確定されていない。
cf) Chromosome evolution at the origin of the ancestral vertebrate genome
https://www.biorxiv.org/content/10.1101/253104v1
cf) Lampreys, the jawless vertebrates, contain only two ParaHox gene clusters
https://www.pnas.org/content/114/34/9146

ゲノム重複は魚類や両生類でも起こると考えられ、有性生殖を行うフナ類や金魚は、52対104本の染色体を持ち、1400万年前にゲノム重複があったと2019年に報告された(他の事例は、アフリカツメガエル)。また、多くの真骨魚類が、最低8つのHox gene clusterを持つらしく、数億年前に、初期の真骨魚類でゲノム重複が起きた可能性も考えられている。節足動物では、一部のクモやサソリに2つのHox gene clusterが見られるようである。

哺乳類だと、仮に4倍体が生じても、子孫を残せる可能性は絶望的と思われてるけど、21世紀になって、哺乳類でも、現在知られる中で最大の染色体数を誇るメンドサビスカーチャネズミ(学名:Tympanoctomys barrerae)及び近縁のゴールデンビスカーチャネズミ(学名:Pipanacoctomys aureus)というマイナーな齧歯類が、4倍体の可能性が議論されている。
Discovery of tetraploidy in a mammal
https://www.nature.com/articles/43815

Molecular cytogenetics discards polyploidy in mammals
https://doi.org/10.1016/j.ygeno.2004.12.004

Evolution of the Largest Mammalian Genome
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5569995/




ヒトの持つ遺伝子の多くが、重複遺伝子として生まれて進化してきたことは確からしいが、全部の遺伝子がそうだと証明されてるわけではない。21世紀になってから、既存遺伝子を雛形としない遺伝子(de novo gene)の生成があるという議論がされるようになった。orphan geneという言葉もあるけど、orphan geneは、遺伝子の水平伝播で獲得した遺伝子でも、そう見えることがあるのに対して、de novo geneは、DNAのnon-coding sequenceから生成されたものと定義される。

de novo gene birthという名前が付いてるけど、ヒト遺伝子の中で、de novo geneだと特定されてるものは今の所ないし、de novo gene birthのメカニズムもよく分かってない。

de novo geneに関する初期の論文として、2006〜7年頃の以下の論文が挙げられる。

Evidence for de novo evolution of testis-expressed genes in the Drosophila yakuba/Drosophila erecta clade
https://doi.org/10.1534/genetics.106.069245

論文中に、We refer to this class of orphan genes as “de novo,” to suggest the possibility that they may derive from ancestrally noncoding sequence.と書かれていて、de novoという接頭辞を使おうと提案したのは彼らっぽい。

彼らは、分子生物学でよく使われるモデル生物Drosophila melanogasterの近縁種Drosophilia yakubaの精巣で発現しているmRNAを解析し、D. melanogaster, D. erecta, D. ananassaeなどの近縁種に相同遺伝子の存在が確認できないD. yakuba固有の遺伝子を探索した。これらの種は、Drosophila melanogaster species subgroupもしくは、Drosophila melanogaster species groupに分類されている。species subgroupとかspecies groupはハエ業界以外ではあんまり使われてないっぽいけれど、属と種の中間の分類階級らしく、とりあえず、Drosophilia melanogasterに最も近縁な種たちを調べたということ。

Wikipediaによると、Drosophilia yakubaに最も近縁なのが、D. santomeaで、次が、D. teissieriとなっている。
Drosophila melanogaster species subgroup
https://en.wikipedia.org/wiki/Drosophila_melanogaster_species_subgroup

こうして見つかった遺伝子は、2〜4個の比較的少数のエクソンからなると予測されている。論文では、これらの遺伝子が、本当に、non-coding sequenceから出現したものであることが証明されてるわけではないし、タンパク質に翻訳されてない可能性もあるが、これらの遺伝子がタンパク質をコードしているならば、そのサイズは、100aa前後で、かなり小さく、突然、巨大な遺伝子が出現する可能性は低いという直感とは整合的である。


2012年の以下の論文では、酵母のORFの内、遺伝子としてアノテーションされてないものを調べている。
Proto-genes and de novo gene birth
https://www.nature.com/articles/nature11184

ゲノム配列からのORFのアノテーションでは、最初に、300塩基を閾値として取捨選択されていて、約6000の遺伝子に対して、約261000のアノテーションされてないORFがあったそうである(過去形だけど、多分、今も状況はそんなに変ってない)(個別に機能が同定された遺伝子の中には、閾値より小さいものもあると思われる)。論文では、長さが30塩基以上のアノテーションされてないORFについてタンパク質に翻訳されてるかどうかを調べた所、108000のORFの内、1139個について翻訳されている証拠を得たとしている("We found that 1,139 of 108,000 ORFs0 show such evidence of translation")。

割合としては1%だけど、そもそもプロモーターがなければ転写もされないので、少ないとも言い切れない(どれくらい転写されてるのかは分からないが)。実験方法としては、リボソームプロファイリングと呼ばれてる方法で、2009年に提案され、次世代シーケンサーを使うので、20世紀には出来なかった実験。但し、翻訳されているというだけで、それが何かの機能を果たしているかは分からない。

"weakly conserved annotated ORF"とアノテーションなしのORFでは、翻訳が確認されるものが1256個、300塩基以下の長さのものは1133個ある(Fig4a及び、aupplementary information fig8 middle)らしいので、100アミノ酸以下の割合が、かなり高い。


2018年の以下の論文は、ちょっと毛色が違って、以前から知られていた遺伝子がde novo geneであると主張している。
De Novo Gene Evolution of Antifreeze Glycoproteins in Codfishes Revealed by Whole Genome Sequence Data
https://dx.doi.org/10.1093%2Fmolbev%2Fmsx311

不凍タンパク質と呼ばれる凍結防止の機能を持つ遺伝子は、様々な生物で見つかっていて、独立に何度も獲得されたと考えられている。論文では、タイセイヨウダラ(Gadus morhua)やコダラ(Melanogrammus aeglefinus)が持つ不凍タンパク質AFGPが、1300〜1800万年前にde novo gene birthで獲得したと推論している。2019年の論文

Molecular mechanism and history of non-sense to sense evolution of antifreeze glycoprotein gene in northern gadids
https://doi.org/10.1073/pnas.1817138116

も、同様の結論だけど、こっちの方が読みやすい(結論だけ見たければ,Fig4)。Fig2(Boreogadus saidaという種のafgp遺伝子)を見ると、3つのエクソンからなり、冒頭2つのエクソンはシグナルペプチドで、3つ目の最も大きな本体は、TAAが繰り返し並ぶ特異な配列をしている。afgp遺伝子は、そこそこ大きいけど、例外的ケースじゃないかと思う。


2019年の以下の論文は、現存する自然界の遺伝子の話ではなく、ランダムに生成した人工的なsmall ORF(長さは10~50aa)によって、大腸菌抗生物質耐性を獲得しうるという実験。
De Novo Emergence of Peptides That Confer Antibiotic Resistance
https://doi.org/10.1128/mbio.00837-19
10の9乗に満たない種類の配列から、3つの候補遺伝子が得られたそうだ。あくまで、抗生物質耐性の有無でスクリーニングしているだけなので、大腸菌にとって有用な機能を果たす遺伝子は、もっと多く含まれている可能性もある。

ヒトに、de novo geneが、どれくらいあるかは分からないが、現在のところ、小さい遺伝子は見過ごされてる可能性も高いので、意外と多い可能性もある。酵母の例を見ると、既知の遺伝子6000に対して、アノテーションされてないが翻訳されるORFが1000以上あり、これらの内、どれくらいが何らかの機能を果たしているのか分からないが、全遺伝子の一割くらい、de novo geneという可能性もある。



ドメインより小さい構造単位が存在するかどうかについて、以下のような論文を書いてる人たちがいる。

Centripetal modules and ancient introns
https://doi.org/10.1016/S0378-1119(99)00292-9

タンパク質の3次元構造に基づいて、各アミノ酸残基に対して定義される"F-value"(とは呼んでないが、勝手に命名した。周辺アミノ酸のα炭素との距離の二乗の平均)の極小点を境界として、centripetal moduleというものを定義して、moduleの境界は、フェーズ0のエクソンイントロン境界(エクソンの始まりや終わりが、コドンの途中になってないような境界)と、よく相関すると主張する論文。

centripetal moduleの概念は、1987年に提案されたものらしいが、この論文では、module境界を計算するプログラムを書いて、統計的に相関があるのを確認したということらしい。プログラムの詳細は説明されてないし、公開されてもいない。

彼らは、ドメインについては特に何も言ってないけど、centripetal moduleは、直感的には、空間的にコンパクトにまとまっている単位に分割しようという発想で、ドメインも、しばしばコンパクトにまとまっているとされているので、ドメイン境界を含むものになるだろうと予想される。


色々疑問はあるけど、例えば、ヒトのピルビン酸キナーゼM2(PKM2)で、エクソンイントロン境界とドメイン境界に、何か関わりがあるか見てみる。

PKM2は、グルコースを2つのピルビン酸に分解する解糖系に於いて、ホスホエノールピルビン酸からADPにリン酸を転移し、ピルビン酸を生成する過程を触媒する酵素である。ヒトは、複数のピルビン酸キナーゼを持つが、その内の一つ。

PKM2(PKM isoform a) domains
http://atlasgeneticsoncology.org/Genes/GC_PKM.html

によれば、PKM2は、4つのドメインからなり、A-domainは、配列上は2つに分断されていて、N-domain (aa 1-43), the A-domain (aa 44-116 and 219-389), the B-domain (aa 117-218) and the C-domain (aa 390-531)らしい。

pkm2遺伝子のエクソンは、
Homo sapiens pyruvate kinase M1/2 (PKM), transcript variant 1, mRNA
https://www.ncbi.nlm.nih.gov/nuccore/NM_002654.6
によれば、11個のエクソンからなるようで、その境界は、タンパク質のドメイン境界と、特に関係してそうには見えない。


PKM2は、解糖系に関わる酵素なので、その進化的起源は、物凄く古い可能性がある。

ヒトPKM2の構造は、
Functional cross-talk between allosteric effects of activating and inhibiting ligands underlies PKM2 regulation
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6636998/
のFig1(A)で見られ、

Evolutionary plasticity in the allosteric regulator-binding site of pyruvate kinase isoform PykA from Pseudomonas aeruginosa
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6802521/
のFig3Bに緑膿菌のピルビン酸キナーゼPykAの構造があるが、目で見る限り、ヒトPKM2に似ている。緑膿菌PykAのA-domainには、高度に保存された配列MVARGDLGVEがあり、触媒活性に重要な役割を果たすと書いてある。ヒトのPKM2のA-domainには、似た配列MVARGDLGIEが含まれている(触媒活性に関与するのかは知らない)


それはそれとして、PKM2のcentripetal moduleを調べることにする。構造データとして、
https://www.rcsb.org/structure/1T5A
を使用して、計算した。横軸は、残基で、縦軸は、F-value。kは、周辺アミノ酸何残基までを計算に入れるか決めるウインドウ幅(の半分。論文では、k=40~90がいいとか書いてるので、同じような範囲で計算した)。
f:id:m-a-o:20201213154649p:plain
は、平滑化なしの値をプロットしたもので、既知のドメイン境界も、重ねて書いてある。

PKM2のドメインは、目で見ても、何となく分かるので当然のことではあるが、ドメイン境界が極小点になってそうだとは分かる。この例では、極小点の中でも、ドメイン境界になっているものは、F-valueがウインドウ幅kの選び方に依存してないようである。いつでも、こんな風にドメイン境界が採れるという保証はない。

生データだと、極小点が多すぎるので、8残基移動平均をとって平滑化して、極小点を抽出してみた。この例だと、あまり、エクソンイントロン境界と相関してそうに見えない。"+"が付いてるのは、フェーズ0で、"-"が付いてるのは、フェーズ0じゃないもの。
f:id:m-a-o:20201213154709p:plain
平滑化したグラフに、ドメイン境界を重ねてみる
f:id:m-a-o:20201213154712p:plain

もうちょっとencourageされる結果だったら、他の例も見てみようと思ったけど、大変微妙な結果。今の所、10〜20残基程度の短い構成単位が存在するという証拠はなさそうに思える。


描画に使用したコード

import matplotlib.pyplot as plt
import numpy as np
from scipy import signal

def readcif(filepath):
   ret = []
   nprev = None
   for line in open(filepath):
       line = line.strip()
       if len(line)<5:continue
       if line[:5]!="ATOM ":continue
       ls = line.strip().split()
       if len(ls)<10:continue
       if ls[0]!="ATOM":continue
       if ls[3]!="CA":continue  #-- not alpha carbon
       if ls[4]!="." and ls[4]!="A":continue
       if ls[6]!="." and ls[6]!="A":continue
       aa = ls[5]
       n = int(ls[8])
       x,y,z=float(ls[10]),float(ls[11]),float(ls[12])
       assert(nprev==None or nprev+1==n),(nprev,n)
       ret.append( (aa,x,y,z) )
       nprev=n
   return ret


def getseq(r):
   dic = {'Ala': 'A', 'Asx': 'B', 'Cys': 'C', 'Asp': 'D', 'Glu': 'E', 'Phe': 'F', 'Gly': 'G', 'His': 'H', 'Ile': 'I', 'Lys': 'K', 'Leu': 'L', 'Met': 'M', 'Asn': 'N', 'Pro': 'P', 'Gln': 'Q', 'Arg': 'R', 'Ser': 'S', 'Thr': 'T', 'Sec': 'U', 'Val': 'V', 'Trp': 'W', 'Xaa': 'X', 'Tyr': 'Y', 'Glx': 'Z'}
   def normalize(s):
      s= s.lower()
      return s[0].upper()+s[1:]
   return "".join([dic[normalize(x[0])] for x in r])



def getF(r , k=50):
   ret = []
   for i,(aa,x,y,z) in enumerate(r):
       j1 = max(0 , i-k)
       j2 = min(len(r)-1 , i+k) + 1
       assert(j2>j1),(j1,j2)
       Fval = sum([(x-x2)**2+(y-y2)**2+(z-z2)**2 for(_,x2,y2,z2) in r[j1:j2]])/(j2-j1)
       ret.append(Fval)
   return ret


def phase(n):
  if n%3==0:return "+"
  else:return "-"


if __name__=="__main__":
  r = readcif("1t5a.cif")
  Nlen = 531  #-- PKM2 has 531aa
  offset = Nlen - len(r) + 1
  xtick = "module"
  smoothing = True
  for k,c in [(40,'r'),(50,'brown'),(60,'grey'),(70,'g'),(80,'b')]:
    d = getF(r,k=k)
    x = np.linspace(offset ,Nlen, len(d))
    y = np.array(d)
    if not smoothing:
      plt.plot(x,y,c,label='k={0}'.format(k))
    else:
      num=8  #移動平均の個数
      b=np.ones(num)/num
      y2=np.convolve(y, b, mode='same')#移動平均
      plt.plot(x,y2,c,label='k={0}'.format(k))
      minpos = signal.argrelmin(y2, order=3)
      plt.plot(x[minpos],y2[minpos],c, marker='o',linestyle='None')
    if xtick=="exon":
       exon_ticks=[((c-89)/3,str(n+3)+phase(c-89)) for (n,c) in enumerate([243,335,467,654,925,1076,1229,1396,1578])]
       plt.xticks([c[0] for c in exon_ticks],[c[1] for c in exon_ticks])
    elif xtick=="module":
       plt.xticks([44,117,219,390])
    plt.yticks([])
    plt.grid(True)
    plt.xlim(0,Nlen)
  plt.legend()
  plt.show()

対数共起頻度のSVD

対数共起頻度を用いた四項類推:word2vecとPMI との比較
https://doi.org/10.11517/pjsai.JSAI2020.0_4Rin177

という報告を読んだ。四項類推は、queen-woman+man ≒ kingみたいなタイプの類推問題。word2vecは、この手の問題によく答えられるということになっている。

四項類推の精度については、通常、最も類似度の高い単語を1〜数個取ってきて、妥当性を確認するという方法が取られ、類似度の値そのものは見てないことが多いけど、ちょっと類似度の値を見てみる。
Pretrained Embeddings
https://wikipedia2vec.github.io/wikipedia2vec/pretrained/

で公開されているword2vecで学習したenwiki_20180420_300dに含まれるデータでは、kingとqueenのコサイン類似度は、元々、0.63ある。kingもqueenも、最上級の支配者を表すという点では、共通していて同じような文脈で出現することが多いと考えられるので、それ自体は問題ない。しかし、king-man+womanとqueenのコサイン類似度は、0.65になるが、単にkingとqueenが、最初から類似度の高い単語であることが効いてるだけのように思える。類似度が上がってるように見えるのも多分偶然の産物で、queen-woman+manとkingのコサイン類似度は、0.60になって、むしろ小さくなる。

しかし、四項類推について全く意味のある情報を含んでいないとも言えない。(king,queen),(man,woman),(father,mother),(god,goddess),(boy,girl),(boys,girls),(son,daughter),(prince,princess),(waiter,waitress)という9組について、差のベクトルの平均を取って、(actor,actress)の差ベクトルとのコサイン類似度を計算すると、0.82という高い値になる。個別の差ベクトルとのコサイン類似度は、0.43〜0.69までバラバラであるが、平均化することで類似度が上がるので、"性差"を表現するベクトルが存在している可能性もある。

actorとactressのコサイン類似度は、元々、0.78ほどあり高いのだけど、actor-"性差"ベクトルとactressのコサイン類似度は、0.92とかなので、0.4~0.7くらいの胡散臭い数値を、うろうろしてるのとは違うようにも思える。個別のペアに関しては、性差以外の部分でも、用法の差があるのかもしれない(例えば、日本語で、"ボーイ"を給仕的な意味に使うことがあり、英語圏でも、この用例は存在するっぽいが、girlで、これに相当する使い方は、多分しないっぽいとか、そういった感じのこと)が、平均化することで、こうした微細な違いが消失するのかもしれない。

gloveの場合も、githubレポジトリのREADMEにリンクがあるglove.42B.300dで実験すると、同じようなことが起きる。例えば、上と同様にして"性差"ベクトルを計算して、actor-"性差"ベクトルとactressのコサイン類似度を見ると、0.89程度になっている。he-"性差"ベクトルとsheのコサイン類似度は、0.93以上ある。
stanfordnlp/GloVe
https://github.com/stanfordnlp/GloVe

gloveに同梱されてるdemo.shでは、小さなコーパスtext8をダウンロードして学習するようになっているけど、text8で学習したベクトルを使うと、同じことをやっても成績は悪いので、大きなコーパスを使うことは重要かもしれない。text8は、Wikipediaを元に作成したコーパスから100MBを抜粋したデータで、約1700万tokensを含む。glove.42B.300dは、420億トークンからなるコーパスで学習しているらしい。Wikipedia+αから作成したコーパス(60億トークン)で学習したglove.6B.300dでも、学習精度は、glove.42B.300dと大差ないように思える。

これらのコーパスは、通常の本などと比べて桁違いに大きく、例えば、King James版英訳聖書は、延べ80万tokens弱(unique wordsは約5400単語)らしい。
単語の分散表現を使った教師なし単語翻訳
https://m-a-o.hatenablog.com/entry/2019/04/01/221230
の末尾で、King James版英訳聖書で学習した分散表現がイマイチっぽいと書いたけど、本1〜2冊程度で学習できる分散表現は、(少なくとも、現在のアルゴリズムでは)やはりカスなんじゃないかと思う


で、冒頭の報告は、対数共起頻度の(trunated) SVDを使っても、四項類推問題で、それなりの成績を収めるということらしい。多くの点で、以下の2015年の論文を参照している。
Improving Distributional Similarity with Lessons Learned from Word Embeddings
https://www.aclweb.org/anthology/Q15-1016/

論文には、(LSAとかで古くから行われてたような)old-styleの"count-based"な手法でも、prediction-basedな手法(≒neuralなんとか)に近い性能を実現することができる的なことが書いてある。対数共起頻度も、count-based modelではあるけど、以下の論文では実験されてないので、それをやったのが冒頭の報告ということ。

対数共起頻度をSVDするだけというのは、理屈上は簡単なので自分でも試してみることにした。上にも書いた通り、コーパスサイズが小さいとゴミになる可能性があるので、報告と同じく、英語Wikipedia全文を利用することにする(少なくとも、word2vecやgloveは、これで、そこそこの性能が出てるらしいし)。報告では、メモリの都合が〜と言い訳しつつ、なんか回りくどいことをしているが、メモリを潤沢に利用できる環境を用意して直接計算した。

Wikipediaの前処理には、fasttextのレポジトリに含まれるwikifil.plを使用して、以下のように、コーパスを作成した。

wget https://raw.githubusercontent.com/facebookresearch/fastText/master/wikifil.pl
curl -L -O http://dumps.wikimedia.org/enwiki/latest/enwiki-latest-pages-articles.xml.bz2
bzcat enwiki-latest-pages-articles.xml.bz2 | perl wikifil.pl > enwiki_all.txt

次に、共起頻度の計算には、gloveで使われているcooccurを流用。デフォルトの設定では、distance-weightingという距離による重み付けが入ってたりするけど、オプションで切ることができる。
GloVe
https://github.com/stanfordnlp/GloVe

gloveのレポジトリに含まれてるdemo.shを改変して、以下のようなシェルスクリプトを使用した。メモリの設定は64GBとかにしてるので、環境に合わせて変更しないとメモリを使い切るかもしれない。shuffleとgloveは、やらなくてもいいけど、ついで。

#!/bin/bash
set -e

make

CORPUS=enwiki_all.txt
VOCAB_FILE=vocab.txt
COOCCURRENCE_FILE=cooccurrence.bin
COOCCURRENCE_SHUF_FILE=cooccurrence.shuf.bin
BUILDDIR=build
SAVE_FILE=glove.raw.300d
VERBOSE=2
MEMORY=64.0
VOCAB_MIN_COUNT=100
VECTOR_SIZE=300
MAX_ITER=15
WINDOW_SIZE=10
BINARY=2
NUM_THREADS=24
X_MAX=10

echo "$ $BUILDDIR/vocab_count -min-count $VOCAB_MIN_COUNT -verbose $VERBOSE < $CORPUS > $VOCAB_FILE"
$BUILDDIR/vocab_count -min-count $VOCAB_MIN_COUNT -verbose $VERBOSE < $CORPUS > $VOCAB_FILE
echo "$ $BUILDDIR/cooccur -memory $MEMORY -vocab-file $VOCAB_FILE -verbose $VERBOSE -window-size $WINDOW_SIZE -symmetric 1 -distance-weighting 0 < $CORPUS > $COOCCURRENCE_FILE"
$BUILDDIR/cooccur -memory $MEMORY -vocab-file $VOCAB_FILE -verbose $VERBOSE -window-size $WINDOW_SIZE -symmetric 1 -distance-weighting 0 < $CORPUS > $COOCCURRENCE_FILE
echo "$ $BUILDDIR/shuffle -memory $MEMORY -verbose $VERBOSE < $COOCCURRENCE_FILE > $COOCCURRENCE_SHUF_FILE"
$BUILDDIR/shuffle -memory $MEMORY -verbose $VERBOSE < $COOCCURRENCE_FILE > $COOCCURRENCE_SHUF_FILE
echo "$ $BUILDDIR/glove -save-file $SAVE_FILE -threads $NUM_THREADS -input-file $COOCCURRENCE_SHUF_FILE -x-max $X_MAX -iter $MAX_ITER -vector-size $VECTOR_SIZE -binary $BINARY -vocab-file $VOCAB_FILE -verbose $VERBOSE"
$BUILDDIR/glove -save-file $SAVE_FILE -threads $NUM_THREADS -input-file $COOCCURRENCE_SHUF_FILE -x-max $X_MAX -iter $MAX_ITER -vector-size $VECTOR_SIZE -binary $BINARY -vocab-file $VOCAB_FILE -verbose $VERBOSE

無事計算を終えると、cooccurrence.binに共起頻度が入っている。(symmetricオプションを有効にしたため)対角成分が、2倍になってたりするけど、どうせ対数を取るので、大きな問題はないと考えて、そのままにしてある

私がやったところ、共起頻度行列の次元は40万弱で、非ゼロ要素の数は、17.8億弱になった。非ゼロ要素の半分近くが1なので、対数を取ると0になるけど、それでも、疎行列を作るだけで10GB以上のメモリが要る(cooccurrence.binのサイズは28GB強)。疎行列のSVDには、RedSVDを使った(RedSVDの実装が正しいのかチェックしてないけど信じることにする)。使用した実装が、Eigenに依存してたので、Eigenも必要。
Eigen
http://eigen.tuxfamily.org/index.php?title=Main_Page

header-only version of RedSVD
https://github.com/ntessore/redsvd-h

実際のプログラムは、以下の通り。UとSを別々に出力してるけど、冒頭の報告では、UとSの平方根を掛けたもの(rootSVDと呼んでる)が一番成績がよいらしい。これは、後から計算する。U単体を、word vectorと見たものは、uSVDと呼んでいる。

#include <cmath>
#include <vector>
#include <iostream>
#include <Eigen/Sparse>
#include <RedSVD/RedSVD.h>

typedef struct cooccur_rec {
    int word1;
    int word2;
    double val;
} CREC;

int main(){
    std::vector< Eigen::Triplet<float> > triplets;
    FILE *fin = fopen("cooccurrence.bin","rb");
    fseek(fin , 0 , SEEK_END);
    size_t sz = (size_t)ftell( fin );
    triplets.reserve( sz / sizeof(CREC) );
    fseek(fin , 0 , SEEK_SET);
    int D = 0;
    while(!feof(fin)){
        CREC rec;
        fread(&rec, sizeof(CREC), 1, fin);
        if ( D < rec.word1 ) D = rec.word1;
        if ( D < rec.word2 ) D = rec.word2;
        if( rec.val > 1.0) {
            triplets.push_back( Eigen::Triplet<float>(rec.word1-1 , rec.word2-1 , logf((float)rec.val)) );
        }
    }
    fclose(fin);
    
    Eigen::SparseMatrix<float> M;
    M.resize(D,D);
    M.setFromTriplets(triplets.begin() , triplets.end());
    triplets.resize(0);
    triplets.shrink_to_fit();

    RedSVD::RedSVD<Eigen::SparseMatrix<float>> svd(M,300);
    auto U = svd.matrixU();
    auto S = svd.singularValues();
    std::cout << "S:" << std::endl << S << std::endl;
    std::cout << "U:" << std::endl << U << std::endl;
    return 0;
}

計算は、gloveより大分早かった(gloveは複数スレッド使って頑張ってるらしいのに数時間、rootSVDやuSVDは1スレッドで10分くらい)。後は、出力を適当に整形すれば、分散表現が得られる。で、四項類推とか以前に、そもそも、(king,queen),(man,woman),(father,mother),(god,goddess),(boy,girl),(boys,girls),(son,daughter),(prince,princess),(waiter,waitress)という9組で、単語のコサイン類似度がどうなってるのか見ると、以下のようになった。

名詞1 名詞2 rootSVD uSVD word2vec glove.42B.300d glove.text8
king queen 0.69 0.40 0.63 0.76 0.52
man woman 0.70 0.30 0.69 0.80 0.57
father mother 0.89 0.69 0.81 0.82 0.65
god goddess 0.64 0.45 0.61 0.47 0.42
boy girl 0.85 0.48 0.77 0.83 0.53
boys girls 0.84 0.65 0.86 0.86 0.54
son daughter 0.86 0.67 0.84 0.82 0.64
prince princess 0.79 0.64 0.68 0.66 0.50
waiter waitress 0.85 0.65 0.72 0.78 0.11

word2vecは、上の方で触れたenwiki_20180420_300dを使った時の数値。glove.text8は、gloveに同梱されてるdemo.shの学習結果で、上に書いたようにコーパスが小さい。雰囲気としては、uSVDやglove.text8は、他よりコサイン類似度が小さい傾向にあるように見える。

rootSVDでは、king-man+womanとqueenのコサイン類似度は、0.64で下がる。king-male+femaleとqueenのコサイン類似度も、0.67で、元々のkingとqueenの類似度が効いてるだけのように見える。性差ベクトルを計算して、king-"性差"ベクトルとqueenのコサイン類似度は、0.75になり改善はする。

網羅的なテストはしてないけど、rootSVDが、四項類推課題に対して、そこそこの性能を持つという報告は正しそうに見える。ただ、それは、元々、類似度の高い単語を取ってるだけでないかという疑惑はある。報告を見る限り、skip-gramよりは多少精度が悪そうなのだけど、共起頻度のみを使う場合、共起単語の出現位置に関する情報(the kingとking theという単語列を区別できない)は完全に捨てられていて、入力の時点で、情報が少ないので、当然という気はする。SVDベースの方法で、このような情報を取り込むには、前方共起頻度と後方共起頻度を別にカウントするとか、もっと詳細には、窓幅4以下での共起と窓幅8以下での共起を分けるとか、複数のcontextを同時に使う方法が思いつくけど、実際有効かどうかは不明

NVIDIA Thrustで分子動力学(1)単原子分子集団のMD

MDは時々使うけど、自分で実装したことがなかったので、一回くらいは、やっておくべき気がした。

実装
https://gist.github.com/vertexoperator/1d2b4c6b71138d59484665cb524143ac

最初なので、相互作用はLennard-Jonesポテンシャルのみで、周期境界条件、NVE(粒子数、体積、総エネルギー)一定という簡単なやつ。粒子の軌跡を追跡するだけでは詰まらないので、任意時点での温度と圧力を計算できるようにしてある。

普通に、全粒子間のLJ相互作用を全部計算すると、粒子数Nに対して、O(N^2)の計算量が必要になるので、分子動力学では、適当なカットオフ距離を決めて、カットオフ距離内の粒子との相互作用のみ計算することが多い(バッキンガムポテンシャルを使う場合も同様の方法が使えると思うけど、クーロン相互作用を計算する場合は、LJ相互作用より減衰しにくいので、ももう少し工夫がいる)。この場合、計算全領域をカットオフ距離幅で分割することで、計算量をO(N)に減らせる。

実装に難しいところはないけど、生のCUDAカーネルを直接書かないという目標は、一箇所だけ、いい方法が思いつかなくて達成できなかった(sortされた整数の配列がある時に、各整数の開始indexを計算する処理で、scan_by_key→unique_by_key→scan→gatherとかで出来るけど、これだけのことに、4つも関数を呼びたくない)。力の計算も、なんか、thrustの正しい使い方じゃない気もするけど、一応正しく、それなりの性能で動いてる。

時間積分は、velocity Verletだけど、今の所、知られてる時間積分では、どの方法であっても、長時間計算すると、エネルギーが発散するという問題は避けられないらしい。エネルギーを保存する時間積分法が欲しいところではある。

カットオフ距離2nm、タイムステップ1fsで、うちのマシンで適当に性能を測定すると、
倍精度、N=100万:1ステップ20msくらい
倍精度、N=1000万:1ステップ140msくらい
単精度、N=100万:1ステップ10msくらい
単精度、N=1000万:1ステップ70msくらい
という感じ。計算は、常温・約8気圧で、カットオフ距離は2nmと、やや大きめにしたものの、単精度と倍精度で、性能比が、ほぼ2:1なことから、残念な感じが察せられる。

力の計算の演算数を数えると、一粒子あたり
・3個の除算(最適化すれば、3個の乗算)と3個のfloor
・8〜11個の加減乗算が近傍グリッドにいる粒子の個数回
・(25個の加減乗算+1個の除算)がカットオフ距離内にいる粒子の個数回
となる。

100万粒子の計算では、系の一辺の長さを、約167nmに設定しているので、粒子が一様に分布していれば、近傍グリッドにいる粒子の個数は、平均46個。カットオフ距離内にいる粒子は(自分自身を除いて)6個強になる。使ったGPUの倍精度演算性能は420Gflopsと公表されていて、100万粒子に対して、420Gflops出れば、これくらいの演算をやるのに、2msあれば十分というくらい。

近傍グリッドは27個あるので、平均すると、1グリッドあたり、1〜2個の粒子しか入ってない。より計算需要があるだろう超高圧の気体や液体になれば、粒子の密度は、数千倍にもなって、一つのグリッドにいる粒子数が増大するので、多分、性能は向上する。希薄気体であっても、全ペアの計算をするよりは、ずっとマシ(100万粒子に対して、全ペアの計算をすれば、10Tflopsの演算性能があっても、1ステップ1秒以上は必要になるはず)なので、こんなんでもいいだろう。


一番はまったのは、初期位置を、単純に一様分布で生成すると、かなり近い位置に粒子が落ちることが結構あって、そうすると、運動エネルギーが普通でも、総エネルギーが異常に大きくなって困る。こういう希薄な気体でも、こういう問題が起きるのは、意外だった。

もう一つ、本来は初期化でやるべき気がすることとして、系の運動量及び角運動量を0にすることがあるが、やっていない。(初期位置を決めれば)運動量と角運動量が0という条件は、速度に関する6個の線形な拘束条件になるので、この拘束条件が6x3N行列Cで書けるとすると、I-C^{T}(CC^{T})^{-1}Cを掛ければ、運動量と角運動量が0の部分空間に落とせる(射影すると、速度分布は、Maxwell-Boltzmann分布からは、ずれると思われる)。LJポテンシャルは方向に依存しない二体相互作用なので、運動量と角運動量は(原理的には)保存するはずである。


NVE条件で計算すると、温度や圧力は別に保存量じゃないので一定値にはならないが、系が平衡状態にあれば、(相転移点の近傍にいない限りは)ある平均値の周りで変動するものと思われる。日常的な感覚に基づく推測でも、外部とエネルギーのやり取りがない孤立系で、突然、温度や圧力が不安定化することは、あんまりないと思われる。けれど、所与の温度や圧力が、どのようなNVE条件で成立するか事前に知ることは、通常、そんなに簡単じゃない。普通、計算したい物理量は、特定の温度や圧力条件下における値なので、これは、少し不便。

なので、分子動力学でも、温度や圧力を制御する方法が色々と考えられている。これらの方法の多くは、1980年台に提案されたものらしい。Nose-Hoover thermostat, Andersen barostat, Berendsen thermostatなどがよく使われるが、実装は、別に難しくない。

何らかのthermostatを使ったNVT計算を、統計力学に於けるカノニカル・アンサンブルの実現と見なすようなことは理論的根拠がないので、これらのbarostat/thermostatは、所与の温度、圧力に対応する密度とエネルギーを発見するために使われる(この値が間違ってないことは理論的に保証されないと思うので、NVE条件による"production run"で、温度や圧力はチェックすべきである)。

生体高分子の計算をやる場合は、NVE一定条件でやると、タンパク質のフォールディングによって温度が変わることがある(十分多くの水分子を用意すれば、系の温度変化は極めて小さくなるだろうが、計算量が増える)ようなので、production runを、NVE条件でやらずに、上に書いたようなthermostat/barostatを使っているのを見ることがあるかもしれない。


実用的な分子動力学に必要な一般化として、多原子分子の計算もある。(古典)分子動力学では、電子の自由度を無視しているので、多原子分子の扱いは、色々無理があるけど、電子の自由度を考慮するには、量子論的な計算が必要になって計算量が爆発的に増えるので、古典分子動力学は、色んな業界で使われている。

電子の自由度を無視した結果として、分子内の相互作用は、二体相互作用の形ではなく、多体相互作用の形でモデル化されることも、しばしばある。こういう相互作用を許すと、運動量や角運動量の保存則が壊れる気がする(そもそも、原子核だけでは、これらの量は保存しないので問題ないとも言える)けど、計算結果は、それなりに有用なものとして扱われる。分子内の相互作用の計算は、実装としては、そんなに難しいものでもない。

分子間の相互作用は、Lennard-Jonesポテンシャルに加えて、Coulomb力が考慮されることが多い。Coulomb力は、LJポテンシャルより、減衰しにくいので、LJポテンシャルのように単純に小さなカットオフを使って遠くの粒子からの影響を切り捨てるわけにはいかないし、それどころか、Coulombポテンシャルは周期境界条件を考慮して、周期対称性を持つ形にしなければならない。つまり、1/r^2ではなく、\sum_{\mathbf{k} \in L} 1/(\mathbf{r} + \mathbf{k})^2の形を取る(Lは周期格子)

Coulombポテンシャルを考慮した計算も、O(N^2)で良ければ、定義どおりに実装してみることは、特に難しくない。計算量を減らすアルゴリズムとしては、particle Mesh Ewald法というものが、よく使われて、O(N log N)になる。最近は、fast multipole method(FMM)を使った実装が増えてきていて、こっちは、O(N)になる。

cf)FMMについては、2020年に実装しましたというだけの論文が出てる程度には枯れてない話題
A GPU-Accelerated Fast Multipole Method for GROMACS: Performance and Accuracy
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7660746/

分子内と分子間の相互作用をモデル化する方法は多分色々提案されているが、特に、水分子については、専用のモデルがいくつも作られている。有名なものとしては、SPC,SPC/E, TIP4P,TIP5Pなどがある。Wikipedia
水モデル
https://ja.wikipedia.org/wiki/%E6%B0%B4%E3%83%A2%E3%83%87%E3%83%AB
を読めば、大体、どういうものかはわかる。

次は、とりあえず、水の分子動力学計算を実装して、いくつかの物理量を計算して、水モデルによる違いでも見ることにする。