Juliaと量子多体計算
スパースモデリングからQuantics Tensor Trainまで

品岡寛
埼玉大学
- 学術変革領域研究B「量子古典融合アルゴリズムが拓く計算物質科学」- JSTさきがけ「量子情報処理」- 基盤研究B「次元圧縮技術が拓く二粒子応答の第一原理計算」

自己紹介

物性の定量的予言がしたい → 計算物理＋第一原理計算

Fortran, Ruby, C/C++, Python, Juliaと使ってきて、今はJuliaがメイン• Python的な気楽な書き方から、Fortran的な速度優先の書き方まで可能• 学生が勉強しやすい• ALPS/CT-HYB C++による量子モンテカルロ法実装 (C++)
https://github.com/ALPSCore/CT-HYB• SparseIR.jl/sparse-ir

虚時間グリーン関数のスパースモデリング (Python/Julia)

https://github.com/SpM-lab/SparseIR.jl
https://github.com/SpM-lab/sparse-ir• DCore 動的平均場近似計算 (Python)
https://github.com/issp-center-dev/DCore• SpM

数値解析接続 (C++)

https://github.com/SpM-lab/SpM

量子古典融合アルゴリズムが拓く計算物質科学 (2023年4月〜)
https://qc-hybrid.github.io/

領域代表: 品岡寛

概要

計算物理は、桁違いに異なる長さスケールの共存を扱う
例 量子物性物理 (低温物理現象 vs 101 eV ~ 105 Kのバンド幅)
宇宙ひも (プランクスケールの細さ vs 宇宙スケールの長さ)
乱流

幅広い長さスケールを扱える表現の開発→試行錯誤＆並列計算

Quantics tensor train

連続変数関数を離散化した「テンソルネ
D=4ットワーク」表現

MSSTA.jl (仮名)スパースモデリング

グリーン関数の悪条件の積分方程式の擬似四倍精度演算
虚時間の最適基底を構成

SparseIR.jl

2009年に応用数学で提案され、ここ数年で流品岡等によって2017年に提案され、虚時間形

体力学、量子力学等への適用が始まる式の場の量子論・第一原理計算に応用が進む

発表の構成• 背景 時空データの圧縮• Quantics tensor train 連続時空の離散化とテンソルネットワーク• スパースモデリング 疑似多倍精度による積分方程式の解法• 今後のJulia応用について

情報圧縮技術がなぜ必要？ (1/2ページ)
有効少数系

量子多体系
Bare Coulomb interaction

Effective interaction
Effective potential
厳密な解法は指数的計算コスト

グランドチャレンジ

密度汎関数 (Kohm-Sham), 摂動論など• より精度の良い近似を捜す (高次の摂動論、機械学習？など)• 少数系ですら計算が難しい →次ページ

情報圧縮技術がなぜ必要？ (2/2ページ)
2つの粒子が散乱される過程
桁違いに異なるスケールの混在→単純メッシュではデータ量が発散→圧縮＆演算する技術が必要(r′′′, t′′′)(r′′, t′′)流体力学計算等でも同じような問題

有効相互作用(r, t)(r′, t′)テンソルネットワーク、スパースモデリング･･･→試行錯誤が迅速で、並列計算までカバーする
Juliaが最適t￼￼￼￼￼￼￼￼￼￼￼￼

Nt = 103, Nr = 1003 → (Nt Nr )4 = 1036

Quantics tensor train• HS et al., Physical Review X 13, 021015 (2023)• M. K. Ritter, Y. N. Fernández, M. Wallerberger, J. von Delft, HS, X. Waintal,arXiv:2303.11819

https://github.com/shinaoka/MSSTA.jl

I. V. Oseledets, Doklady Math. 80, 653 (2009)

Quantics tensor train (QTT)とは？ B. N. Khoromskij, Constr. Approx. 34, 257 (2011)
Image compression

桁違いに異なる長さスケールが混在する関数を圧縮可能

D=1

D=4D=8Bond dimension

N. Gourianov et al., Nat. Comput. Sci. 2, 30 (2022)
José I. Latorre, arXiv:quant-ph/0510031v1

Vlasov-Poisson equations for collisionless plasmas
E. Ye and N. F. G. Loureiro, arXiv:2205.11990

場の量子論の時空依存性を圧縮・演算できる！HS et al., PRX 13, 021015 (2023)

2022年8月にアイデア思いつく、2022年10月に投稿 (Juliaだから可能だった)テンソルネットワーク (1/3)
スカラーベクトルa

3階テンソル行列

aijk

aijaijiiijk

突き出している足がindexに対応。突き出している向き、図形の形・色に意味はない。

iテンソルネットワークAjx=∑jAiαBjβAij xj

共通している足に対する和を取る。

Ckα α,β β=Ai Bj Ck∑
α,β上・下付きにも意味はない。

テンソルネットワーク (2/3)
Contraction

ABCijk→ijk

実際にvirtual bondに対する和を計算することをcontractionをとるという。行列積として実行可能

テンソルネットワーク (3/3)
特異値分解

A = UΣVA=†

大きさは特異値の数UV

†

=UΣV

†

(Σ)ij = δijsi
特異値は非負。急速に減衰する場合は、行列の低ランク近似が可能

行列積表現 (matrix product state) / tensor train
≈

x1xN

≈

x1xN

⋯

x1

分解は特異値分解などで実行可能。最悪の場合、virtual bondの大きさは、指数的に大きい。
→ 多くの物理系、データにおいて、小さな次元ですむことがある。
例: 1次元量子系

離れた場所の間での「量子エンタングルエントロピー」が小さい

xN

ITensors.jl

Flatiron instituteで開発が進められている

M. Fishman, S. R. White, E. Miles
Stoudenmire, SciPost Phys.
Codebases 4 (2022)• MPSからより汎用的なテンソルネットワークまで• 階層構造- バックエンド: NDTensors.jl (型安定性を保った高速なライブラリ)- MPSコード: 動的で柔軟なコード• GPU対応 (multi GPU?)A = UΣV†==i j kU

Vi kj

https://itensor.github.io/ITensors.jl/stable/examples/ITensor.html

QTTのまとめ (1/3ページ)

I. V. Oseledets, Doklady Math. 80, 653 (2009)
B. N. Khoromskij, Constr. Approx. 34, 257 (2011)

Scalar-valued function f(x) with 0 ≤ x < 12

Binary coding: x = (0.x1x2⋯xR)2 = x1 /2 + x2 /2 + ⋯ + xR /2f(x)(0.0⋯0)2 (0.0⋯1)2(0.1⋯1)2x10f(x1, ⋯, xR) ≈D1

∑⋯xR

⋯

α1=1

≈x1

Rx1

α1DR−1∑αR−1=1(1)(2)(R)
Fx ,1α Fx1,α1α2⋯Fx ,α 1
11R R−1⋯⋯αR−1xRQTTのまとめ (2/3ページ)
Same length scalef(x, x′)≈x1⋯xRx = (0.x1x2⋯xn⋯)2 ∈ [0,1)x′ = (0.x′1x′2⋯x′n⋯)2 ∈ [0,1)x′1⋯

x′Rx1x′1⋯

xR−1 x′R−1 xR￼￼￼￼￼￼￼￼￼￼

同程度の長さスケールの変数が、強くエンタングルしている場合x′R

QTTのまとめ (3/3ページ)
指数関数f(x) = e2−xn−x1/2 −x2/2−xn/2=ee⋯e⋯x = (0.x1x2⋯xn⋯)2 ∈ [0,1)

D=1

The Sum of N exponential functions can be represented as a QTT of rank at most N.
∵ Bond dimensions are added when MPSs are added.

多項式単位行列D≤1+pM. Ali and A. Nouy, “A. Approximation Theory of Tree Tensor Networks:Tensorized Univariate Functions”, Constr Approx (2023)

f(x, y) = δx,y = δx1,y1δx2,y2⋯
00=⨂⨂

D=1⨂⨂

⨂

⋯

Representation of Continuous Functions: https://tensornetwork.org/functions/

Written by Miles Stoudenmire

場の量子論データの圧縮

HS et al., Physical Review X 13, 021015 (2023)

1P momentum space

2P vertex functions

Spectral function(a) 0.8

Ω(!)

0.6

AF insulator
Doped Mott insulator

0.4
0.2
0.0°50!

5(b)

Ω(!)

0.4

Impurity model

0.0°0.025 0.000 0.025

0.2

0.0

0.2°50!

5

Bond dimension

(c)
103
102

既存のデータを特異値分解で圧縮してみた

Multipolar susceptibility

多くの時空データが数桁に渡り圧縮可能

AF insulator
Doped Mott insulator
Impurity model

Nonequilibrium system (real-time Green’s function)101
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
Bond b

CeB6

圧縮したまま(場の量子論)計算ができるか？
Rを増やすことで、指数的な解像度向上 (計算量はRに線形)ikr̂
F(r) = dk F(k)e∫

フーリエ変換量子フーリエ変換回路＝古典コンピュータで効率的に実行可能

K. J. Woolfe et al., Quantum Inf. Comput. 17, 1 (2017), J. Chen et al.,arXiv:2210.08468v1

C(t, t′′) = dt′A(t, t′)B(t′, t′′)∫

行列積

C(t1, t′1′, ⋯, tR, t′R′) =∑t′1,⋯,t′R

A(t1, t′1, ⋯, tR, t′R)B(t′1, t′1′, ⋯, t′R, t′R′)
Tensor network contractions
t1tR⋯

A

￼

￼

￼

￼

￼

￼

t′R =t′1′t′R′

￼

￼

￼

￼

￼

￼

￼

￼

￼

￼

t′1

￼

B

￼

￼

￼
￼

￼

￼

￼

￼

￼

￼

ほぼ任意のダイアグラム計算が実行可能δδ

Juliaライブラリの実装 MSSTA.jl (multiscale space-time ansatz)
まだ、実験的なライブラリです (○ユニットテスト、△ドキュメンテーション)
https://github.com/shinaoka/MSSTA.jl• フーリエ変換などの演算を簡単に実装可能なライブラリ• バックエンドとして、ITensors.jlを利用Juliaをなぜ選んだのか？
• forループの使用が避けられない• 戦略的にJuliaで場の量子論計算のエコシステムを作る• (ITensors.jlの開発者と交流があった)

VS Codeを使ってデモしてみます！

QTTの展望• 分散並列計算、GPU対応• より複雑なテンソルネットワークへの拡張• 皆さんの周りに、桁違いの長さスケールが共存する難しい問題ありませんか？
QTTで非線型常微分方程式も解けます Cahn-Hilliard方程式?

• Quantics Tensor Cross Interpolation (TCI): 特異値分解の代替
M. K. Ritter, Y. N. Fernández, M. Wallerberger, J. von Delft, HS, X. Waintal, arXiv:2303.11819

宣伝: ICIAM2023

Quantics表現の宣伝予定！

スパースモデリング

疑似多倍精度による積分方程式の解法「虚時間」に特化した最適な基底の構成疑似多倍精度演算: MultiFloats.jl疑似四倍精度• 倍精度小数点2つの和として、浮動小数点を表現 (例: 1 + 10-100は、有効桁数が16桁程度の倍精度では表現不能)• 真の四倍精度のソフトウェアエミュレーションよりも速い (でも、倍精度演算の2桁ぐらい遅い)

Float64は桁落ち

桁落ちしていない

小数点以下31桁まで正しい

MultiFloats.jl vs DoubleFloats.jl
ヒルベルト行列

Hij = 1/(i + j − 1)

1 1/2 1/3
1/2 1/3 1/4
1/3 1/4 1/5

共に疑似四倍精度浮動小数点演算を実装• SVDは、MultiFloats.jlの方が数割速い• MultiFloats.jlは、四倍精度より高次の実装を含む

MultiFloats.jl

39 ms

DoubleFloats.jl

62 ms

cf. Float64 0.5 ms

二桁の速度差Julia v1.91, MultiFloats v1.0.3, DoubleFloats v1.2.4
on Apple M1 Pro (1 thread)虚時間形式の相関関数とは？
相関関数⟨A(r, t)B(r′, t′)⟩

違う時空点間での物理量の相関を測る。
久保公式によって、物理系の外場応答と対応。

長時間極限 | t − t′| で速く減衰しない。→数値的に圧縮しづらい

虚時間形式(松原形式とも呼ぶ)

数値的なダイアグラム計算でよく用いられる。

t = iτ

時間を仮想的に純虚数に取る
長時間極限 | τ − τ′| で指数的に減衰する。→数値的に圧縮し易い (が、情報も一部失われている)

￼￼

￼

￼

“スパースモデリング” = 虚時間での最適基底の構成

虚時間形式の相関関数における情報損失
G(τ) = − ⟨Tτ A(τ)B(0)⟩ =虚時間の量∞∫−∞dωK(τ, ω)ρ(ω)
実周波数(時間)の量

系に寄らない積分核−τωeK(τ, ω) =−βω1+eフェルミオンの場合高周波数の情報をτに伝搬しにくい

β = 1/T

虚時間の相関関数は、実時間の相関関数に比べて、高周波数の情報が落ちている。
→

積分核を使って、虚時間の最適な正規直交基底を作れるか？

積分核の特異分解による基底の構成(a) 1

ωmax = 1

HS, J. Otsuki, M. Ohzeki, K. Yoshimi, PRB 96, 035147 (2017)

l=0−410SlUl(τ)Vl(ω)

for| ω | < ωmax, 0 < τ < β = 1/T

Orthonormal intermediate-representation bases

G(τ) =∑
l=0GlUl(τ)ρ(ω) =∑
l=0−810−1210ρlVl(ω)

Gl = − Sl ρl• Dimensionless parameter Λ≡βωmax• | Gl | decay as fast as Sl ( ∵ | ρl | is bounded).
• Logarithmic growth of Sl > ϵ w.r.t Λ

Number of Ncoefficients∑

β = 100
β = 1000
β = 10000

FSlK(τ, ω) =∞(02550l75100

Matsubara frequency with tailLegendre
O(β)
Chebyshev
104
Legendre and Chebyshev
IRO( β)

10i!n3O(log β)

102
101

Intermediate representation (IR)102104106

Ø

小さな特異値に対応する特異ベクトルを倍精度で数値的に決めたいਤ 56 ඞཁͳ‫ج‬ఈͷ਺ N ͷԹ౓ґଘੑ

応用の拡がり

量子多体計算から第一原理計算まで

Analytic continuation

Migdal-Eliashberg theory

Sparse modeling (SpM)• J. Otsuki, M. Ohzeki, HS, K. Yoshimi, PRE 95, 061302(R) (2017)• [Application to lattice QCD] E. Itou and Y. Nagai, JHEP 07, 007(2020)• T. Wang, T. Nomoto, Y. Nomura, HS, J. Otsuki, T. Koretsune, and R.
Arita, PRB 102, 134503 (2020)

Nevalinna/Carathéodory formalism

Self-energy embedding calculations• J. Fei, C.-N. Yeh, E. Gull, PRL 126, 056402 (2021)• J. Fei, C.-N. Yeh, D. Zgid, and E. Gull, PRB 104, 165111 (2021)

All electron GW calculations• S. Iskakov, C.-N. Yeh, E. Gull, and D. Zgid, PRB 102, 085105 (2020)• C. N. Yeh, S. Iskakov, D. Zgid, and E. Gull, PRB 103, 195149 (2021)• K. Haule, S. Mandal, arXiv:2008.07727

FLEX• N. Witt, E. G. C. P. van Loon, T. Nomoto, R. Arita, T. O. Wehling,
PRB 103, 205148 (2021)• N. Witt, J. M. Pizarro, T. Nomoto, R. Arita, T. O. Wehling,
arXiv:2108.01121v1

Estimation of magnetic interactions• T. Nomoto, T. Koretsune and R. Arita, PRB 102, 014444 (2020)• T. Nomoto, T. Koretsune and R. Arita, PRL 125, 117204 (2020)• Y. Nomura, T. Nomoto, M. Hirayama and R. Arita, PRR 2,043144 (2020)

Extension to the two-particle level• HS, J. Otsuki, M. Ohzeki, K. Yoshimi, K. Haule, M. Wallerberger, E.
Gull, PRB 97, 205111 (2018)• HS, D. Geffroy, M. Wallerberger, J. Otsuki, K. Yoshimi, E. Gull, J.
Kuneš, SciPost Phys. 8, 012 (2020)• M. Wallerberger*, HS*, A. Kauch, PRR 3, 033168 (2021)

独自実装や、Pythonライブラリの使用も含みますSparseIR.jlM. Wallerberger, S. Badr, S. Hoshino, S. Huber, F. Kakizawa, T. Koretsune, Y. Nagai,
K. Nogaki, T. Nomoto, H. Mori, J. Otsuki, S. Ozaki, T. Plaikner, R. Sakurai, C. Vogel,
N. Witt, K. Yoshimi, H. Shinaoka, SoftwareX 21, 101266 (2023)

https://github.com/SpM-lab/SparseIR.jl• 基底関数の区分多項式 (Legendre)による展開: 特異分解→特異値分解• 疑似四倍精度演算による特異値分解 MultiFloats.jl + GenericLinearAlgebra.jl• 基底の構成が終われば、後は倍精度演算で十分Python実装 (sparse-ir)との比較
Pure Python 実装: https://github.com/SpM-lab/sparse-irnumpyの疑似四倍精度浮動小数点拡張 https://github.com/tuwien-cms/xprec
C言語 + Python
M. Wallerberger, S. Badr, S. Hoshino, S. Huber, F. Kakizawa, T. Koretsune, Y. Nagai, K. Nogaki, T. Nomoto, H.
Mori, J. Otsuki, S. Ozaki, T. Plaikner, R. Sakurai, C. Vogel, N. Witt, K. Yoshimi, H. Shinaoka, SoftwareX 21, 101266(2023)

パッケージのリリース時に、バイナリーの自動ビルド＋登録• condaとPyPIを両方メンテナンスする必要• ビルドに使うGitHub Actionsのランナーとcondaの互換性 (glibcのバージョンなど)• 普通の物理学者、学生にはメンテナンスが難しい･･･SparseIR.jlのデモンストレーション

小さな特異値まで精度良く計算

SparseIR.jlのデモンストレーション

基底関数の精度は大体10桁チュートリアル集• 典型的な計算のサンプルコード• Jupyter Book
• 手元でのビルド環境は、
VS Code + Dockerで用意
https://spm-lab.github.io/sparseir-tutorial/解説 + コード例by 永井さん

最後に: Juliaによる計算物理の展望
Juliaを使うことで、試行錯誤のスピードを上げることが出来た！
今後挑戦したい• テンソルネットワークの分散並列計算 Distributed vs MPI• Multi GPUs

研究員募集中!
- 学術変革領域研究B「量子古典融合アルゴリズムが拓く計算物質科学」(領域代表品岡: 2023〜2025年)- JST創発「2粒子レベルの量子埋め込み理論に基づく新規第一原理計算手法の開発と実証」(2024〜2026年: フェーズ1、2027〜2030年: フェーズ2)• 3年任期

Juilaを使った時空圧縮技術の開発と応用、物性物理外からの挑戦も歓迎！

• 特任助教の名称の付与可• 助教級の待遇 (1年目: 480〜530万円、2年目: 530〜580万円)• 別途、家賃補助、通勤手当あり