マス・フォア・インダストリ研 究 No.27

数学と物理における
Juliaの活用

編集

富谷
横山
永井
寺崎

昭夫
俊一
佑紀
敏志

Institute of Mathematics for Industry
Kyushu University

About the Mathematics for Industry Research
The Mathematics for Industry Research was founded on the occasion of the certification
of the Institute of Mathematics for Industry (IMI), established in April 2011, as a MEXT
Joint Usage/Research Center – the Joint Research Center for Advanced and
Fundamental Mathematics for Industry – by the Ministry of Education, Culture, Sports,
Science and Technology (MEXT) in April 2013. This series publishes mainly proceedings
of workshops and conferences on Mathematics for Industry (MfI). Each volume includes
surveys and reviews of MfI from new viewpoints as well as up-to-date research studies
to support the development of MfI.
October 2022
Kenji Kajiwara
Director
Institute of Mathematics for Industry

Use of Julia in Mathematics and Physics
Mathematics for Industry Research No.26, Institute of Mathematics for Industry, Kyushu University
ISSN 2188-286X
Editors: Akio Tomiya, Shun’ichi Yokoyama, Yuki Nagai, Satoshi Terasaki
Date of issue: 18 January 2024
Publisher:
Institute of Mathematics for Industry, Kyushu University
Motooka 744, Nishi-ku, Fukuoka, 819-0395, JAPAN
Tel +81-(0)92-802-4402, Fax +81-(0)92-802-4405
URL https://www.imi.kyushu-u.ac.jp/
Printed by
Kijima Printing, Inc.
Shirogane 2-9-6, Chuo-ku, Fukuoka, 810-0012, Japan
TEL +81-(0)92-531-7102 FAX +81-(0)92-524-4411

数学と物理における Julia の活用

編集

：

富谷

昭夫

横山

俊一

永井

佑紀

寺崎

敏志

巻頭言:
Julia 言語 は 2018 年にバージョン 1 が公開されたオープンソースの科学技術計算言語で
あり、一言でまとめると Fortran の様に高速でかつ Python の様に生産性の高い言語であ
る。とはいえ、Fortran や Python、Mathematica 言語や MATLAB 言語など科学技術計算
が可能なプログラミング言語はすでに数多ある。なぜ今 Julia 言語なのだろうか。Julia
言語は後発であるがゆえ、様々な言語の良いところを集め、また反省点が生かされてい
る。後発ゆえの利点は数多いが例えば JIT (Just in time)コンパイルシステムは、Python で
も Numba や JAX などを通じて実装が進んでいるが Julia 言語では言語の開発段階から想
定されている。また多重ディスパッチによる関数の多重定義は、様々なデータ型に対し
ての操作を自然な実装で実現することができる枠組みとなっている。LLVM を用いたコ
ンパイルも後発ゆえの選択だろうと思う。加えて開発コミュニティがユーザーに近いた
めユーザーの声が届きやすく、多くのパッケージ開発もユーザーコミュニティで進んで
いる。Julia 言語の難点といえば、様々な分野において活用が始まっているが、2023 年
時点では、Python や C/C++ の様に専門家が各機関にいるということもなく、最新の情
報はいくつかの書籍かネットを参照するしかないのが現状である。そこで本研究会は、
産学を含め、数学・物理学の最前線で研究を行っている研究者に講演を頂いた。数学・
物理学を含んだ広い分野における様々な利用例が全国的に共有され、議論も活発に行わ
れた。本研究会に参加した方、資料を参照した方の中から Julia 言語を用いた数学・物
理学の研究が生まれれば至極幸せである。
富谷昭夫

－i－

材料科学においては，近年，急速に必要となっている幾何・代数の材料科学への適応に関
わる研究の更なる加速が期待される．本研究会はその礎・足場として期待に応えるものとな
ったと考え，ここに報告をする．
組織委員代表 富谷昭夫
2024年１月18日

組織委員

－ii－

富谷昭夫

大阪国際工科専門職大学

横山俊一

東京都立大学

永井佑紀

日本原子力研究開発機構

寺崎敏志

AtelierArith

Contents
巻頭言 ・ ・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・

i

Julia 入門

後藤俊介（来栖川電算）, 寺崎敏志 (AtelierArith)・・・・・・・・・・・・・

1

BasicBSpline.jl で始める B-spline
堀川由人（大阪大学）・・・・・・・・・・・・・・・・・・・・・・・・

89

応用数学の概念を用いた物理モデルをいかに Julia でプログラミングを行うか
降籏大介（大阪大学）・・・・・・・・・・・・・・・・・・・・・・・・ 123
Julia と量子多体計算 : スパースモデリングから Quantics Tensor Train まで
品岡寛（埼玉大学）・・・・・・・・・・・・・・・・・・・・・・・・・ 133
数論における Julia の援用
横山俊一（東京都立大学）・・・・・・・・・・・・・・・・・・・・・・ 151
Julia による科学技術計算 : 大規模並列計算について
永井佑紀（日本原子力研究開発機構）・・・・・・・・・・・・・・・・・ 169
素粒子物理における Julia の活用〜格子 QCD の大規模計算に向けて〜
富谷昭夫（大阪国際工科専門職大学）
・・・・・・・・・・・・・・・・・ 195

－iii－

概要:
題目: Julia 入門
講演者:後藤俊介 (来栖川電算) , 寺崎敏志 (AtelierArith)
概要: ライトユースにも使いたい、けれどもそれで高速性を犠牲にしたくない。総合的
に開発できるものにしたい。超初心者にも習得は容易でありながら、ハッカーの満足に
も応えられるものがほしい。」そういった願いからプログラミング言語 Julia が誕生し
ました。直近では v1.9 がリリースされ、長年課題であったパッケージの読み込み時間
の短縮など、開発体験の改善が行われました。さらに和書の入門書も充実してきてお
り、学習を始めるには絶好のタイミングといえるでしょう。本チュートリアル講演では
Julia の入門として 1. インストールと基本的な使い方, 2. 型と多重ディスパッチ, 3. スレ
ッド並列の基礎 の 3 つのトピックを解説します。この機会に Julia に入門し、新しい体
験を始めてみませんか。
題目: BasicBSpline.jl で始める B-spline
講演者: 堀川由人（大阪大学）
概要: B-spline は滑らかさを保証した区分多項式を作るための道具であり、幾何形状表
現や数値解析など広範な分野で活用されている。Julia における B-spline のパッケージは
幾つか存在するが、その中でも BasicBSpline.jl を本講演では扱う。BasicBSpline.jl パッケ
ージの使用方法を交えながら B-spline に関する理論を解説し、その応用例についても触
れる。
題目:応用数学の概念を用いた物理モデルをいかに Julia でプログラミングを行うか
講演者: 降籏大介（大阪大学）
概要: 物理上の相分離現象に対するモデル方程式としては Cahn-Hilliard 方程式が有名か
つ有用だが，この数値解析は計算量的に高コストである．そこでこの現象のごく初期過
程を除いて状況を記述するようなシンプルなモデルを考案したい. この際，流体的挙動
を模するために空間の Voronoi 分割という応用数学的手法をモデルに取り入れるのだ
が，これまでのコンピュータ言語ではこうしたモデルの数値解析にはそれなりに苦労す
る面が多々あったのが現状である．これに対し，Julia を用いると比較的容易に数値解
析が行えることを示すことで，微分方程式の数値解析に対して Julia がどのように用い
られるかを紹介したい．
題目: Julia と量子多体計算：スパースモデリングから Quantics Tensor Train まで

－v－

講演者: 品岡寛（埼玉大学）
概要:本講演では、Julia を用いた情報圧縮技術の量子多体計算への応用を紹介する。 前
半では、虚時間グリーン関数の「スパースモデリング」技術 [1]を実装した SparseIR.jl
[2] を紹介する。本ライブラリでは、疑似四倍精度浮動小数点演算による特異値分解を
活用し、悪条件の積分方程式を解くことで、虚時間グリーン関数の最適な正規直交基底
を構成する。これにより、虚時間形式に基づく場の量子論計算を大幅に高速化すること
ができる [2]。後半では、桁違いに異なる複数の長さスケールが共存する関数を圧縮可
能な quantics (quantized) tensor train (QTT) [3]を紹介する。近年、QTT は、流体力学[4]や
量子多体計算[5]など、物理学の様々な分野で応用されている。テンソルネットワー
ク、QTT を解説した後、品岡等が開発する実験的なライブラリのデモンストレーション
を行う。このライブラリを用いていることで、QTT によって圧縮されたデータのフーリ
エ変換、畳み込みなどの演算が簡便に行える。量子物理学を超えた将来の応用について
も議論したい。
[1] H. Shinaoka, J. Otsuki, M. Ohzeki, K. Yoshimi, PRB 96, 035147 (2017); H. Shinaoka et al.,
SciPost Physics Lecture Notes, 063 (2022). [2] M. Wallerberger, S. Badr, …, H. Shinaoka,
SoftwareX 21, 101266 (2023). [3] I. V. Oseledets, Doklady Math. 80, 653 (2009); B. N.
Khoromskij, Constr. Approx. 34, 257 (2011). [4] N. Gourianov et al., Nat. Comput. Sci. 2, 30
(2022). [5] H. Shinaoka et al., Physical Review X 13, 021015 (2023).
題目: 数論における Julia の援用
講演者: 横山俊一 (東京都立大学)
概要: Julia 言語は現在急速に普及しているが、数学分野においては比較的幾何学・解析
学の分野においての活用例が多く、代数学の分野においてはまだまだ少ないのが現状で
ある。しかしながら近年では AbstractAlgebra.jl などに代表される、非常に汎用性の高い
パッケージが充実してきている。本講演では代数系分野における Julia 活用の可能性に
ついて触れた後、とくに数論における Julia native の数式処理システム開発プロジェクト
NemoCas / OSCAR について紹介する。
題目: Julia による科学技術計算:大規模並列計算について
講演者: 永井佑紀 (日本原子力研究開発機構)
概要: Julia を使った科学技術計算の方法について具体例を交えながら解説を行う。特
に、現在 Fortran を用いて数値計算を行ってきている方が Julia に期待しているであろう
点について言及する予定である。Fortran は、高速でありかつスーパーコンピュータ上
での大規模並列計算が可能である。一方で、機械学習に関するパッケージはほとんどな

－vi－

く、Fortran 使用者にとって、既存の数値計算と機械学習を組み合わせるのは非常に難
しい。そこで、Julia である。Julia は高速で動作し、MPI 並列計算も Fortran（と同程度
かそれ以上に楽に）実装できる。そして、最新の機械学習アルゴリズムを簡単に実装で
きる。本講演では、Julia による大規模並列計算の方法と、機械学習の方法について述
べる予定である。
題目: 素粒子物理における Julia の活用～格子 QCD の大規模計算に向けて～
講演者: 富谷昭夫 (大阪国際工科専門職大学)
概要: 格子 QCD の数値計算周りの話と Julia の利用について話します。

－vii－

チュートリアル講演
チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

この資料について
数学と物理にるJuliaの活⽤ でのチュートリアル講演資料
piever/Remark.jl を使って Markdown 原稿をWライドとして表⽰

更新⽇: 2023-07-10T03:59:43.607

本⽇の資料は AtelierArith/julia_tutorial に;2．

2 / 87

－1－

ここでの Julia とは
プログラミング⾔語の⼀つ
2012 年 Why We Created Julia に発表. 開発は 2009 年ごろから
Jeff Bezanson, Stefan Karpinski, Viral B. Shah, and Alan Edelman
In short, because we are greedy.

冒頭で
個⼈的に好きな箇所

とる9に既存のプログラミング⾔語の良いところを取り⼊れた⾔語

We want the speed of C with the dynamism of Ruby
We want something as usable for general programming as Python
Something that is dirt simple to learn, yet keeps the most serious hackers happy.
(Did we mention it should be as fast as C?)

要するに
⾼級⾔語の9に使い7すく計算機の上で⾼速に動作する
速く動くRードを早くを書くこと できる
3 / 87

Why We Use Julia, 10 Years Later
年は 10 周年記念 �� Why We Use Julia, 10 Years Later が公開
リポジトリで管理されている
この原稿 にみんなが集まって書いたもの
例: miguelraz さんの Pull Request など
Julia というプログラミング⾔語に出会ったきっかけが紹介されている. 例えば以下の9うなもの:
Why We Created Julia を読んでソースをビルドした
同僚，職場，講義で出会った
C/C++/Fortran は難しい..., Python は簡単に使えるけれど遅い...
気づいたら特定のj^Qージのメンテナになっていた
2022

JuliaLang/www.julialang.org

4 / 87

－2－

Star History
Star History
40.0k

JuliaLang/julia

30.0k

sratS buHtiG

20.0k

10.0k

2012

2014

2016

2018

2020

Date

2022
star-history.com

5 / 87

直近の動き
がリリース
TTFX 問題が改善される !!!
可視化ツールの使⽤が捗る！
デバッグまわりのツールが増えてきた
鈍器(褒め⾔葉)扱いの 実践Julia⼊⾨， Juliaプログラミング⼤全が登場してきた
い>A%本が今年は出版される:しい
Juila 1.9

6 / 87

－3－

時代は Julia なのでは？

そこでチュートリアル講演です
概要: ライトユースにも使いたいけれどもそれで⾼速性を犠牲にしたくない。総合的に開発できるものにしたい。超初
⼼者にも習得は容易でりながらhッカーの満⾜にも応えられるものがほしい。そいた願い らpグラwン
グ⾔語 Julia が誕⽣しました。直近では v1.9 がリリースされ⻑年課題でたjッQージの読み込み時間の短縮など
開発体験の改善が⾏われました。さらに和書の⼊⾨書も充実してきており学習を始めるには絶好の\イwングといえ
るでし8。本チュートリアル講演では Julia の⼊⾨として
1. インストールと基本的な使い⽅ (<– このスライcの担当箇所)
2. 型と多重ディスjッチ
3. スッc並列の基礎
3 つのトピッOを解説します。この機会に Julia に⼊⾨し新しい体験を始めてみまA 。
7 / 87

Julia

のインストール
8 / 87

－4－

Julia

のインストール

要するに julia というコマンドが利⽤できれば良い
$ julia

_
_
_ _(_)_ | Documentation: https://docs.julialang.org
(_) | (_) (_) |
_ _ _| |_ __ _ | Type "?" for help, "]?" for Pkg help.
| | | | | | |/ _` | |
| | |_| | | | (_| | | Version 1.9.1 (2023-06-07)
_/ |\__'_|_|_|\__'_| | Ofﬁcial https://julialang.org/ release
|__/
|
julia>

すでに使える⼈はスキップする

素直な⽅法
Julia

の公式サイトから⼊⼿ https://julialang.org/downloads/ 移動し各⾃の環境に合わせて導⼊

9 / 87

それはそうだが
任意の⼈間が「各⾃の環境に合わせて導⼊してください」がで れ*苦労し%い
各⾃の環境が何とか, 環境変数とかjス@通すとか
バージョンが上が<度に⼿動でインストールしたいか?

10 / 87

－5－

それはそうだが
任意の⼈間が「各⾃の環境に合わせて導⼊してください」がで れ*苦労し%い
各⾃の環境が何とか, 環境変数とかjス@通すとか
バージョンが上が<度に⼿動でインストールしたいか?

そA%%たに Juliaup
Juliaup - The Julia installer and version multiple…
multiple…

11 / 87

Juliaup - Julia version manager

を使う⽅法 (1)

Windows
パワーシェルを開く．
PS> winget install julia -s msstore

Windows

アプリケーシ~ ら も⼊⼿ができる

Julia in the Windows Store | …

12 / 87

－6－

Juliaup - Julia version manager
juliaup/deploy/shellscript /juliaup-init.sh

を使う⽅法 (2)

を叩いている．

Mac and Linux
$ curl -fsSL https://install.julialang.org | sh
$ source ~/.bashrc
$ julia --version
julia version 1.9.1

--yes

オプションを使えば途中の対話操作を省略し進めることができる. command-line-arguments を39．

$ curl -fsSL https://install.julialang.org | sh -s -- --yes

13 / 87

動作確認
インストールができたと仮定して進める． julia というコマンドを使うことができるか確認をする:
$ julia --version
julia version 1.9.1

何も指定し%い場合 REPL(Read-Eval-Print Loop) が起動する:
$ julia

_
_
_ _(_)_ | Documentation: https://docs.julialang.org
(_) | (_) (_) |
_ _ _| |_ __ _ | Type "?" for help, "]?" for Pkg help.
| | | | | | |/ _` | |
| | |_| | | | (_| | | Version 1.9.1 (2023-06-07)
_/ |\__'_|_|_|\__'_| | Ofﬁcial https://julialang.org/ release
|__/
|
julia>

14 / 87

－7－

Appendix: Julia

⾃体のアップデート

というアナウンスが出た
経由で⼊れていると次の9う%通知が出る

Julia v1.9.2 has been released
Juliaup

$ julia
The latest version of Julia in the `1.9` channel is 1.9.2+0.x64.apple.darwin14. You currently have `1.9.1+0.x64.apple.darwin14` installed. Run:
juliaup update
to install Julia 1.9.2+0.x64.apple.darwin14 and update the `1.9` channel to that version.
$ juliaup update
Installing Julia 1.9.2+0.x64.apple.darwin14
Downloading: [>
] 1.35 MiB/114.18 MiB eta: 14m

15 / 87

Appendix: Juliaup

以外の選択肢

abelsiqueira/jill
jill - Julia Installer 4 Linux - Light
johnnychen94/jill.py
pip3 install jill && jill install
Docker
docker run -it --rm julia:1.9.1

を⽤いて

の9&てR `dを起動する

16 / 87

－8－

参考資料
系が現時点で最新安定版. もし 1.10 がースされると 1.9 系のサポート（バグ修正など）はしなくなる．
特定のバージョンを⻑く使いたい場合は LTS (1.6 系) を使うと良い
1.9

Julia’s Release Process

17 / 87

Julia

を動かす (REPL)
18 / 87

－9－

なぜ REPL で動かすのか?
もちろんファイルにコードを記述し次のよに動かしても良い．
$ julia script.jl

は JIT コンパイルによって動作する⾔語
するとにパ^Qージのロード関数のコンパイルが⾏われるそのためのコストは無視できない
コンパイル結果を使い回し効率よく作業する必要がある
試⾏錯誤時は REPL の上で作業をするのがよく⾏われる
REPL の上でで学べることが多い
Python と異なりインaントに関してYンU`ィブではないのでコードをコmーペーストにより⾃由に実⾏
させることができる
もちろん VS Code, Jupyter, Pluto な$の^]な環境でも良い
Julia

julia script.jl

19 / 87

VS Code

を使いたい場合

のホームページ
エクステンションの導⼊
いくつか便利なショートカットNー
ファイルにフォーカスを当てて Shift + Enter を押す
Alt-J Alt-O で REPL を開く
VS Code

Pluto.jl

使いたい場合

Pluto.jl

⼊⾨

JuliaTokai

で話た勉強会資料
20 / 87

－10－

Julia

を動かす (REPL)

versioninfo()

を使うと詳細な情報を得ることができる．バグレポートなどに添付すると良い．

$ julia

_
_
_ _(_)_ | Documentation: https://docs.julialang.org
(_) | (_) (_) |
_ _ _| |_ __ _ | Type "?" for help, "]?" for Pkg help.
| | | | | | |/ _` | |
| | |_| | | | (_| | | Version 1.9.1 (2023-06-07)
_/ |\__'_|_|_|\__'_| | Ofﬁcial https://julialang.org/ release
|__/
|
julia> versioninfo()
Julia Version 1.9.1
Commit 147bdf428cd (2023-06-07 08:27 UTC)
Platform Info:
OS: macOS (x86_64-apple-darwin22.4.0)
CPU: 16 × Intel(R) Core(TM) i9-9900K CPU @ 3.60GHz
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-14.0.6 (ORCJIT, skylake)
Threads: 1 on 16 virtual cores
Environment:

Julia

21 / 87

を動かす (REPL)
22 / 87

－11－

REPL

を使う (Julian mode)

julia> 1 + 1 #
2

簡単な算数

julia> println("Hello World")
Hello World
julia> print("Hello"); print(" "); print("World") #
Hello World
julia>

julia>

セミコロンで繋げることもできる

の部分はプロンプトと呼ばれる. julia> の部分も含めてコピペしても REPL 側が適切に処理する．

julia> 1 == 2
true #
false
julia> println("Hello")
World #
Hello

実際は なで実際に動かすと false になる
実際は が出るはず

と異なりインデントに関してYンU`ィブではない
コードを適にコピペして⾃由に実⾏させることができる
Python

23 / 87

Example
正の整数 a, b を選ぶ．gcd(a, b) が 1 となる確率が 1/ζ(2) = 6/π2 となるとを数値計算で確認する様⼦
julia> function calcπ(N)
cnt = 0
for a ∈ 1:N # ∈ \in + <tab>
for b in 1:N # `in`

は

と⼊⼒する. `in` でも良い
の代わりに `=` と書いても良い
とすることもでる

if gcd(a, b) == 1 # gcd(a, b) == 1 && (cnt += 1)
cnt += 1
end
end
end
prob = cnt / N / N
return √(6/prob) # √ \sqrt + <tab>
end
calcπ (generic function with 1 method)
julia> @time calcπ(10^4)
1.791594 seconds
3.141534239016629 #

は

円周率に近い値

上と同じコードを Python で書いてみると⾯⽩いかも．少なくとも数⼗秒かかる．)

(

－12－

24 / 87

REPL

を使う (Help mode)

例えば x⁽¹⁾, xₖ, θ̂, ⊻ のような⽂字も使う#ができる
どう⼊⼒するべきなのか? REPL が教えてくれる!!!
ここで を押す

julia> #
?
help?> ?θ̂
"θ̂" can be typed by \theta<tab>\hat<tab>

通常の Julian mode で呼び出す#もできる．
julia> (Base.Docs.doc)((Base.Docs.Binding)(Main, :θ̂))

上記のRードの代わりに下記のようにて @doc マクロを使う#もできる：
julia> @doc θ̂

25 / 87

寄り道(マクロについて)
マクロは呼ばれる式の断⽚を受け取り新しい式を⽣成する.
例えば @doc θ̂ は (Base.Docs.doc)((Base.Docs.Binding)(Main, :θ̂)) を⽣成する. その様⼦は次の9うにして確認できる:
julia> @macroexpand @doc θ̂
:((Base.Docs.doc)((Base.Docs.Binding)(Main, :θ̂)))

⼊⼒の仕⽅の他に，関数の使い⽅を知ることもできる:
help?> ndims
ndims(A::AbstractArray) -> Integer
Return the number of dimensions of A.
See also: size, axes.
Examples
≡≡≡≡≡≡≡≡≡≡
julia> A = ﬁll(1, (3,4,5));
julia> ndims(A)
3

26 / 87

－13－

寄り道(マクロについて)
1

の 原始 3 乗根 ω = exp(2π −1/3) をそのままコードに落とし込むと次の9にな<:
は虚数単位

julia> ω = exp(im * 2π/3) # im
-0.4999999999999998 + 0.8660254037844387im
julia> @show ω
ω = -0.4999999999999998 + 0.8660254037844387im
-0.4999999999999998 + 0.8660254037844387im
julia> ω ^ 3 # 3
1
0.9999999999999998 - 6.106226635438361e-16im
julia> @assert ω ^ 3 ≈ 1 # "≈" can be typed by \approx<tab>

乗するとほぼほぼ になる．

マクロを使わないと次の9なコードを書いてい<とにな<:
julia> ω = exp(im * 2π/3)
-0.4999999999999998 + 0.8660254037844387im
julia> println("ω = ", repr(begin dummyvariable = ω end)); dummyvariable
ω = -0.4999999999999998 + 0.8660254037844387im
-0.4999999999999998 + 0.8660254037844387im
julia> if ω ^ 3 ≈ 1
nothing
else
Base.throw(Base.AssertionError("ω ^ 3 ≈ 1"))
d

27 / 87

寄り道(フォントの問題)
⽂字化けする場合は⼊⼒に⽤いるフォントを変更するとよい． 例えば次の候補がある:
cormullion/juliamono
tohgarashi/JuGeM
yuru7/juisee
miiton/Cica

例えば macOS の場合 ~/Library/Fonts/ 直下に .ttf 拡張⼦のフBEルを配置する．
三⾓関数", title="⽇本語", fontfamily="JuGeM-Regular")

julia> using Plots; plot(sin, label="

28 / 87

－14－

Pkg mode
パッケージマネージ{ 付属している．
julia> using Pkg; Pkg.add("Example") # JuliaLang/Example.jl
julia> using Example; Example.hello("World")
julia> hello("Example")
julia> Pkg.rm("Example") #

をインストールする

アンインストール

Pkg mode

を⽤いて次のようにして書くこともできる:

を⼊⼒する

julia> # ]
(@v1.9) pkg> add Example
(@v1.9) pkg> ^C # Ctrl C
julia> #

元に戻る

と を同時に押す

29 / 87

寄り道(Example パッケージについて)
上記のコードは実際に動作する例. Example.jl の実装は概ね次のように%っている. 詳細は こちら を参照せよ.
module Example
export hello, domath
hello(who::String) = "Hello, $who"
domath(x::Number) = x + 5
end

export hello

と宣⾔しているので using Example によって hello を即時使うことができる．

julia> using Example
julia> hello("World")
"Hello, World"
julia> domath(3)
8

30 / 87

－15－

名前空間の話
によってユーザが name を使うことができる
便利である⼀⽅「この関数はどこで定義しているのか+っと⾒わからない」問題がある
export name

julia> using Example: hello
julia> hello("World") # hello
julia> domath(3) #

がロードされていることがわかる
これはできない
こように使う

julia> Example.domath(3) #

julia> using Example: Example
julia> import Example
julia> Example.hello("World")

その他下記のようにしてk トを得ることができる
julia> @doc hello("world")
julia> @which hello("world")
julia> @less hello("world")

31 / 87

TerminalClock.jl
Example.jl

は簡単すぎるかもしれない？

julia> using TerminalClock; clock()

JuliaPackaging/Preferences.jl

によって時計の⽂字盤を制御する# " る:

$ julia -q
julia> using TerminalClock
julia> tomlﬁle = joinpath(dirname(pathof(TerminalClock)), "dials", "UnicodeBox.toml");
julia> TerminalClock.set_dials(tomlﬁle)
julia> exit() # REPL
$ julia -q
julia> using TerminalClock; clock()

を再起動

32 / 87

－16－

Search mode
で⼊⼒した過去の履歴を検索できる. ~/.julia/logs/repl_history.jl に残っている
各々 Ctrl-R, Ctrl-S でSーチができる
⼊⼒途中の式 + Ctrl-P で以前⼊⼒した6のを補完してくれる
REPL

OhMyREPL.jl

を使う場合

KristofferC/OhMyREPL.jl

を使うと直感的に探すことができる:

julia> using OhMyREPL
julia> 1 + 1 #
julia> # Ctrl-R

シンタックスハイライトされる
で直感的な⼊⼒履歴を参照でる

33 / 87

Julia

で開発する際のワークフロー
34 / 87

－17－

Julia

で開発する際のワークフロー

の機能を触れたので Julia での開発ワークフローについて触れる
知る##で⼈⽣が豊かになる（か6）
REPL

35 / 87

REPL

周り学習リソース
2022 Workshop

Revise.jl

Julia REPL Mastery Workshop | Miguel Raz Guzmá…
Guzmá…

の使い⽅v1.9 で追加された機能に関しての解説も

Julia: Work ow tips, and what is new in v1.9 | 2nd D…
D…

36 / 87

－18－

Project.toml/Manifest.toml
⼤抵のプログラムは何かしらのパッケージ(ライブラリ)に依存
ユーザは何をインストールすればいいかを知りたい
Project.toml に!使⽤するパッケージを記述
Manifest.toml はより詳細な依存関係の情報を格納 (Project.toml から⾃動⽣成される)
Pkg.activate("path/to/Project.toml")

プロジェクト(環境)をCクティベート(活性化)する．
Pkg.instantiate()

から Manifest.toml を作る．
から依存関係を解決
と書くこともできる（混乱がなれば Project.toml でよい）．

Project.toml
Manifest.toml
JuliaProject.toml

37 / 87

例: 研究会の⽇程，講演者のアoWトを確認する (1)
$ cd julia_tutorial/playground/code/mdtable

table.jl

を動かすた5の依存関係を解決

$ julia -e 'using Pkg; Pkg.activate("."); Pkg.instantiate()'

または --project=@. で現在のaィレクトリまたは親のaィレクトリ&る Project.toml か:定まるpVェク
トをアクティベートする．Stack Overﬂow での解説も参考
--project,

$ julia --project=@. -e 'using Pkg; Pkg.instantiate()'

38 / 87

－19－

例: 研究会の⽇程，講演者のCoWbを確認する (2)
特定の環境で実⾏する際にも --project または --project=@. が必要.
$ julia --project=@. table.jl 0

--project=@.

を指定するのは⾯倒なので環境変数を設定するnBEル (.bashrc や .zshrc) のな に

export JULIA_PROJECT=@.

を⼊れておくと幸せになれる．

39 / 87

試⾏錯誤の⽅法
を書く
を実⾏する # (´・ω・)起動に時間がかかる ��
を更新する
julia script.jl を実⾏する # (´・ω・)起動に時間がかかる ��
script.jl

julia script.jl
script.jl

JIT

コンパイルが毎回⾛るので（⼈間にとって）効率が悪い．コンパイル結果を使い回す運⽤が必要．
内部に main 関数がるとする．

mylib.jl

julia> include("mylib.jl")
julia> main()
#
julia> include("mylib.jl")
julia> main()
#

別のターミナルで作業して mylib.jl を編集
別のターミナルで作業して mylib.jl を編集

これでも良いが Revise.jl を使うと良い（次のrージへ）

－20－

40 / 87

Revise.jl
mylib.jl

内部に main 関数があるとする．

julia> using Revise
julia> # includet Revise
julia> includet("mylib.jl")
julia> main()

は

から export されている

これで mylib.jl を変更すれば変更の main() が実⾏できる． パッケージの開発で6同様．
using Revise; includet("mylib.jl") の9にYwR ; で繋げて書くと１⾏でYッbアップができるので作業効率が良
い．
VS Code にて Alt-J Alt-O にて REPL を開くと Revise がすでにードされている
Alt-J Alt-W でGJッ]すると? る
41 / 87

Debugger.jl
⽂字通りデバッガ
１⾏１⾏実⾏しその都度変数の状態を監視できる
VS Code で6できる REPL の上で6できる
julia> using Debugger
julia> include("mylib.jl")
julia> @enter main()

VS Code

での使い⽅はこちら

42 / 87

－21－

JET.jl
型不安定なコードや潜在的なIーを検知ができる．
開発のモチベーションは aviatesk/grad-thesis な$ :知ることができる
⽇本語で読むことができる

43 / 87

型安定・型不安定の話 (1)
は JIT コンパイル⽅式で動作をする
Julia は関数の引数に渡された値の型の情報 を元にコンパイルをする
⼊⼒の型から出⼒の型が決定できれ*型安定な実装となり効率の良いコードを⽣成することができる
⼊⼒される値によって出⼒の型が変わると型不安定になる. 速度が求められる箇所では型安定なコードを書く
の必須です.
型安定なコードを書くように意識すると
型システムに親しめるようになる
型システムに親しめるようになる
型安定なコードを書くように意識する
型システムにいては]|ートリアル第⼆部で⾏い2す．（このスライドは第⼀部です）
Julia

44 / 87

－22－

型安定・型不安定の話 (2)
What does "type-stable" mean?
Type annotation make JIT compile faster?

上記の質問に対する Stefan Karpinski さんの回答:
No. You do not generally need type annotations on function arguments (except to control behavior via
dispatch), nor do you need type annotations in local scope. The place that type annotations are essential
for performance is on locations: the ﬁelds of structs and the element types of arrays and other data
structures.

実際 Argument-type declarations にる9に関数の引数に対して型アノテーションが必要な理由は Dispatch,
Correctness, Clarity のでて，実⾏速度の理由で必要とするわけでない．

45 / 87

型安定・型不安定の話 (3)
系のvクロで Julia がどの9&Rードを理解しているか観測がで る
系のvクロの説明は Stack Overﬂow での議論
What is the difference between @code_native, @code_typed and @code_llvm in Julia? の解説が?か;7い．
@code_xxx
@code_xxx

46 / 87

－23－

型不安定な例 (1)
ReLU (rectiﬁed linear unit)

の例が典型的

relu1(x) = x > 0 ? x : 0

上記のRードは下記と等価
function relu2(x)
if x > 0
return x
else
return 0
end
end

これは @code_lowered で確認することができる．playground/code/code_comparision の_ールを⽤いて⽐較ができる．
47 / 87

型不安定な例 (2)
relu1(x) = x > 0 ? x : 0

の型 T が Float64 の場合, x の値によって Int 型である 0 を返すかもしれないし T を型とする x ⾃⾝を返すかもしれ
ない．この曖昧さが型不安定を引き起こす．

x

48 / 87

－24－

型不安定な例 (2)
relu1(x) = x > 0 ? x : 0

の型 T が Float64 の場合, x の値によって Int 型である 0 を返すかもしれないし T を型とする x ⾃⾝を返すかもしれ
ない．この曖昧さが型不安定を引き起こす．

x

処⽅箋
0

の代わりに zero(x) とする. これ) x の型に応じた適切なZ(加法単位元)を与えてくれる.
Get the additive identity element for the type of x

relu3(x) = x > 0 ? x : zero(x)

このよにすると x の⼊⼒の型が T だった場合，常に T を型とする出⼒を与える実装になる(型安定になる)．
49 / 87

型安定な例 (1)
relu3(x) = x > 0 ? x : zero(x)

⼊⼒ x は数学的には実数 R の要素であることに注意する．これは x の型 T が Julia の Real のSoタイプであること
を要請する:
relu4(x::Real) = x > 0 ? x : zero(x)
relu5(x::T) where {T<:Real} = x > 0 ? x : T(0)
relu6(x::T) where {T<:Real} = x > 0 ? x : zero(T)

そして relu6 は等価
Julia 内部では次のよな実装になてる (@less zero(1.0) な$で検証せよ):
relu4, relu5,

zero(x::Number) = oftype(x,0)
oftype(x, y) = convert(typeof(x), y)
zero(::Type{T}) where {T<:Number} = convert(T,0)
convert(::Type{T}, x::Number) where {T<:Number} = T(x)::T

50 / 87

－25－

型安定な例 (2)
@code_llvm relu4(rand())

の結果

; @ string:1 within `relu4`
deﬁne double @julia_relu4_2276(double %0) #0 {
top:
%.inv = fcmp ogt double %0, 0.000000e+00
%1 = select i1 %.inv, double %0, double 0.000000e+00
ret double %1
}

@code_llvm relu5(rand())

の結果

; @ string:2 within `relu5`
deﬁne double @julia_relu5_2278(double %0) #0 {
top:
%.inv = fcmp ogt double %0, 0.000000e+00
%1 = select i1 %.inv, double %0, double 0.000000e+00
ret double %1
}

51 / 87

型不安定・型安定なコードの⽐較 (1)
using Random
"""

型不安定なコード

"""
function main1(N=10)
rng = MersenneTwister(0)
s = 0.0
arr = [relu1(2rand(rng) - 1) for _ in 1:N] # arr
for y in arr
s += y
end
s
end

は Vector{Real} になる

"""

型安定なコード

"""
function main3(N=10)
rng = MersenneTwister(0)
s = 0.0
arr = [relu3(2rand(rng) - 1) for _ in 1:N] # arr
for y in arr
s += y
end

は Vector{Float64} になる
52 / 87

－26－

型不安定・型安定なコードの⽐較 (2)
ほんのちょっとの気遣いで 10 倍コードが速くなる・ちょっとした怠けで 10 倍遅くなる．
julia> using BenchmarkTools
julia> N = 100000;
julia> @benchmark main1($N) #
BenchmarkTools.Trial: 1391 samples with 1 evaluation.
Range (min … max): 3.283 ms … 7.546 ms
GC (min … max): 0.00% … 44.61%

型不安定
┊

┊

Time (median): 3.403 ms
GC (median): 0.00%
Time (mean ± σ): 3.590 ms ± 595.631 μs
GC (mean ± σ): 4.67% ± 10.10%

┊
▃▅█▅▃▂▁
███████▇▅▅▄▄▁▃▃▁▃▁▃▁▁▁▅▃▄▅▆▆▃▆▅▅▅▅▅▃▄▅▅▃▃▄▄▁▅▅▄▅▅▆▄▄▅▅▅▁▆▅

3.28 ms Histogram: log(frequency) by time 6.19 ms <
Memory estimate: 3.83 MiB, allocs estimate: 149983.
julia> @benchmark main3($N) #
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
Range (min … max): 214.528 μs … 9.073 ms
GC (min … max): 0.00% … 96.93%
Time (median): 249.382 μs
GC (median): 0.00%

型安定

┊

┊

┊ GC (mean ± σ): 24.49% ± 12.41%
█▂
▁
██▆▆▄▁▁▃▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃▁▃▃▄▁▁▁▃▅▄▄▄▄▃▃▅▄▅▄▅▅▄▅

Time (mean ± σ): 335.062 μs ± 635.965 μs

215 μs

Histogram: log(frequency) by time

4.45 ms <

53 / 87

Memory estimate: 800.88 KiB, allocs estimate: 14.

型不安定・型安定なコードの⽐較 (3)
code_warntype, JET.report_opt

な$"検出する. 対応するvクロもある．

using JET

⾊々警告が出る．

だと警告は⾚⾊で表⽰される

#
REPL
code_warntype(main1, (Int,))
@code_warntype main1(10) # (InteractiveUtils.code_warntype(main1, (Base.typesof)(10)))
report_opt(main1, (Int,))
@report_opt main1(10) # JET.report_opt(main1, (Base.typesof)(10))
code_warntype(main3, (Int,))
@code_warntype main3(10)
report_opt(main3, (Int,))
@report_opt main3(10)

54 / 87

－27－

「型安定・型不安定の話 (2)」の続き
Type annotation make JIT compile faster?

上記の質問に対する Stefan Karpinski さんの回答:
No. You do not generally need type annotations on function arguments (except to control behavior via
dispatch), nor do you need type annotations in local scope. The place that type annotations are essential
for performance is on locations: the ﬁelds of structs and the element types of arrays and other data
structures.

要するに
構造体のフィールド
パフォーマンスの観点からは (Vector{Real} ではなく Vector{Float64} とするのが良い)
配列をはじめとするデータ構造に関しての要素型
55 / 87

構造体: アフィン変換の例 (1)
アフィン変換 f (x) = W x + b を考える
julia> struct Afﬁne
W
b
end
julia> (aff::Afﬁne)(x) = aff.W * x + aff.b # Python

での `def __call__` 相当

julia> aff = Afﬁne(rand(2,2), rand(2))
Afﬁne([0.05221655288726834 0.4531635136343529; 0.006660217186564732 0.32819491295200176], [0.16834154818251978, 0.6642250801
julia> aff(x)
2-element Vector{Float64}:
0.24022734713529842
0.7082428857823445

動いてそう.
56 / 87

－28－

構造体: アフィン変換の例 (2)
⼈間のからすると W は⾏列, b はqクトルということはわかるが， Julia 側からはわからない:
前のスライドからの続き
フィールドの型情報が⾒えない

julia> #
julia> typeof(aff)
Afﬁne #

julia> @code_warntype aff(x)
MethodInstance for (::Afﬁne)(::Symbolics.Arr{Num, 1})
from (aff::Afﬁne)(x) @ Main REPL[2]:1
Arguments
aff::Afﬁne
x::Symbolics.Arr{Num, 1}
Body::Any
1 %1 = Base.getproperty(aff, :W)::Any
%2 = (%1 * x)::Any
%3 = Base.getproperty(aff, :b)::Any
%4 = (%2 + %3)::Any
return %4

─
│
│
│
└──

フィールcの型情報が⾒えないので Any がい+い．(パフォーマンスの観点から9 ない)

57 / 87

構造体: アフィン変換の例 (3)
フィールドに型をつければいいんでしょ？
⼀旦 REPL を再起動して下記を動かす
julia> struct Afﬁne
W::Matrix{Float64}
b::Vector{Float64}
end
julia> (aff::Afﬁne)(x) = aff.W * x + aff.b
julia> aff = Afﬁne(rand(2,2), rand(2))
Afﬁne([0.6711341109848074 0.24382214292980786; 0.6036132680677607 0.13610672620163822], [0.2609069141726297, 0.4285554314320
julia> aff(rand(2))
2-element Vector{Float64}:
0.9884996898337165
1.04891268376067

動いてそう
58 / 87

－29－

構造体: アフィン変換の例 (4)
良さそうにみえる
julia> @code_warntype aff(rand(2))
MethodInstance for (::Afﬁne)(::Vector{Float64})
from (aff::Afﬁne)(x) @ Main REPL[2]:1
Arguments
aff::Afﬁne
x::Vector{Float64}
Body::Vector{Float64}
1 %1 = Base.getproperty(aff, :W)::Matrix{Float64}
%2 = (%1 * x)::Vector{Float64}
%3 = Base.getproperty(aff, :b)::Vector{Float64}
%4 = (%2 + %3)::Vector{Float64}
return %4

─
│
│
│
└──

59 / 87

構造体: アフィン変換の例 (5)
次の使い⽅をするユーザには適⽤できない ��
julia> struct Afﬁne
W::Matrix{Float64}
b::Vector{Float64}
end
julia> Wf32 = rand(Float32, 2, 2); bf32 = rand(Float32, 2);
julia> Afﬁne(Wf32, bf32).W |> typeof
Matrix{Float64} (alias for Array{Float64, 2}) # Float32
julia> using CUDA
julia> Wcu = CUDA.rand(2,2); bcu = CUDA.rand(2,2)
julia> Afﬁne(Wcu, bcu).W |> typeof
julia> Matrix{Float64} (alias for Array{Float64, 2}) # GPU
julia> W = OffsetArray(rand(2,2), 0:1, 0:1)

で計算したいのに...
のリソース使いたかったのに．．．

アフィン変換は汎⽤性が⾼いので多くの場⾯で"いい感じに"振る舞ってほしい!
60 / 87

－30－

構造体: アフィン変換の例 (6)
次のよにてパラメトリック構造体(Parametric Composite Types) によって定義する:
julia> struct Afﬁne{T1, T2}
W::T1
b::T2
end
julia> (aff::Afﬁne)(x) = aff.W * x + aff.b
julia> Wf32 = rand(Float32, 2, 2); bf32 = rand(Float32, 2);
julia> aff = Afﬁne(Wf32, bf32);
julia> @code_warntype aff(rand(Float32, 2))
MethodInstance for (::Afﬁne{Matrix{Float32}, Vector{Float32}})(::Vector{Float32})
from (aff::Afﬁne)(x) @ Main REPL[7]:1
Arguments
aff::Afﬁne{Matrix{Float32}, Vector{Float32}}
x::Vector{Float32}
Body::Vector{Float32}
1 %1 = Base.getproperty(aff, :W)::Matrix{Float32}
%2 = (%1 * x)::Vector{Float32}
%3 = Base.getproperty(aff, :b)::Vector{Float32}
%4 = (%2 + %3)::Vector{Float32}
return %4

─
│
│
│
└──

61 / 87

構造体: アフィン変換の例 (7)
今回の9& 2x2, 3x3 程度の⼩規模の⾏列では配列のSイズ情報も含めて最適なコードを⽣成する StaticArrays.jl を
⽤いた⽅法を採⽤する#良い．
julia> using StaticArrays
julia> struct Afﬁne{T1, T2}
W::T1
b::T2
end
julia> (aff::Afﬁne)(x) = aff.W * x + aff.b
julia> W = @SMatrix rand(Float32, 2,2)
2×2 SMatrix{2, 2, Float32, 4} with indices SOneTo(2)×SOneTo(2):
0.106584 0.369698
0.182812 0.803011
julia> b = @SVector rand(Float32, 2)
2-element SVector{2, Float32} with indices SOneTo(2):
0.8716512
0.4812594
julia> aff = Afﬁne(W, b)
julia> x = @SVector rand(Float32, 2)
julia> aff(x)

62 / 87

－31－

寄り道 (StaticArrays.jl による計算の様⼦)
⼈間が直接 W[1,1] * x[1] + W[1,2] * x[2] のよに書き下したのとほぼ同じことをしていることがわかる:
前のページの続き

julia> #
julia> @code_typed aff(x)
CodeInfo(
1 %1 = Base.getﬁeld(aff, :W)::SMatrix{2, 2, Float32, 4} # W = aff.W
%2 = StaticArrays.getﬁeld(%1, :data)::NTuple{4, Float32} # W.data <-

─
│
⽣の8ータにアクセス
│ %3 = Base.getﬁeld(%2, 3, false)::Float32 # W.data[3] == W[1,2]
│ %4 = StaticArrays.getﬁeld(x, :data)::Tuple{Float32, Float32} # data = x.data
│ %5 = Base.getﬁeld(%4, 2, false)::Float32 # x.data[2] == x[2]
│ %6 = StaticArrays.getﬁeld(%1, :data)::NTuple{4, Float32} # data = W.data
│ %7 = Base.getﬁeld(%6, 1, false)::Float32 # W.data[1] == W[1,1]
│ %8 = StaticArrays.getﬁeld(x, :data)::Tuple{Float32, Float32} # x.data
│ %9 = Base.getﬁeld(%8, 1, false)::Float32 # x.data[1]
│ %10 = Base.mul_ﬂoat(%7, %9)::Float32 # W[1,1] * x[1] 相当の計算
│ %11 = Base.muladd_ﬂoat(%3, %5, %10)::Float32 # W[1,2] * x[2] + W[1,1] * x[1]
│ %12 = StaticArrays.getﬁeld(%1, :data)::NTuple{4, Float32}
(中略)
└── return %30
) > SVector{2 Float32}

63 / 87

寄り道 (ConcreteStructs.jl でサボる)
julia> struct Afﬁne{T1, T2}
W::T1
b::T2
end

における T1, T2 は形式的にける6のでて実装の部分としては本質的でない．実装時に⼈間が意識したくない．
マクロを提供する

julia> using ConcreteStructs # @concrete
julia> @macroexpand @concrete struct Afﬁne; W; b; end

下記のRードと等価
struct Afﬁne{__T_W, __T_b} <: Any
W::__T_W
b::__T_b
function Afﬁne(W::__T_W, b::__T_b) where {__T_W, __T_b}
return new{__T_W, __T_b}(W, b)
end
end)

64 / 87

－32－

寄り道 (サボり⽅)
型を意識しなくて済む.
julia> using ConcreteStructs; @concrete struct Afﬁne; W; b; end
julia> using StaticArrays
julia> W = @SMatrix rand(2,2);
julia> b = @SVector rand(2)
2-element SVector{2, Float64} with indices SOneTo(2):
0.620866754177038
0.5075755639223726
julia> aff = Afﬁne(W, b)
Afﬁne{SMatrix{2, 2, Float64, 4}, SVector{2, Float64}}([0.21089312982809838 0.4062630292552363; 0.039518291001160444 0.68047617533

今までの議論から分か<9&，サボ<&)=なりの教養が必要．
65 / 87

JET.jl

を⽤いた潜在的なエラーの発⾒ (1)

"""
sumevens(N::Integer)

から までの範囲の偶数を収集する 収集した値を全て⾜し，その値を返却する.

1
N
.
"""
function sumevens(N::Integer)
N ≥ 1 || throw(DomainError(N, "`N` cannot be less than 1."))
arr = []
for n in 1:N
if iseven(n)
push!(arr, n) # Python
`arr.append(n)`
end
end
return sum(arr)
end

での

相当

julia> @assert sumevens(8) == 2 + 4 + 6 + 8 == 20

66 / 87

－33－

JET.jl

を⽤いた潜在的なエラーの発⾒ (2)
おっとっと？

julia> sumevens(1) #
ERROR: MethodError: no method matching zero(::Type{Any})

の時は arr = [] の次にあるループは実質何もしない. ⼀番最後で sum([]) を実⾏することになる．
は Vector{Any} を型とする要素数が 0 の配列
Vector{Any} は任意の値が格納で る. 任意の型に対するZ元を定義するのがで ないので Julia ではエラ
ーを返す
no method matching zero(::Type{Any}) が出るのはこのため

N=1
[]

は⼤雑把にいえば下記のようなことをする
初期値 これをどう定義するか？

# sum(arr)

s=<
>#
for a in arr
s += a
end

67 / 87

JET.jl
JET.jl

を⽤いた潜在的なエラーの発⾒ (3)

は sumevens(1) を実⾏しなくてもまずい部分を検出してくれる．

julia> using JET; report_call(sumevens, (Int,))
1 possible error found
sumevens(N::Int64) @ Main ./REPL[1]:13
sum(a::Vector{Any}) @ Base ./reducedim.jl:994

═════
═════
┌
│┌
││┌ sum(a::Vector{Any}; dims::Colon, kw::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}) @ Base ./reducedim.jl:994
│││┌ _sum(a::Vector{Any}, ::Colon) @ Base ./reducedim.jl:998
││││┌ _sum(a::Vector{Any}, ::Colon; kw::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}) @ Base ./reducedim.jl:998
│││││┌ _sum(f::typeof(identity), a::Vector{Any}, ::Colon) @ Base ./reducedim.jl:999
││││││┌ _sum(f::typeof(identity), a::Vector{Any}, ::Colon; kw::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}) @
│││││││┌ mapreduce(f::typeof(identity), op::typeof(Base.add_sum), A::Vector{Any}) @ Base ./reducedim.jl:357
││││││││┌ mapreduce(f::typeof(identity), op::typeof(Base.add_sum), A::Vector{Any}; dims::Colon, init::Base._InitialValue) @ Base
│││││││││┌ _mapreduce_dim(f::typeof(identity), op::typeof(Base.add_sum), ::Base._InitialValue, A::Vector{Any}, ::Colon) @ Bas
││││││││││┌ _mapreduce(f::typeof(identity), op::typeof(Base.add_sum), ::IndexLinear, A::Vector{Any}) @ Base ./reduce.jl:427
│││││││││││┌ mapreduce_empty_iter(f::typeof(identity), op::typeof(Base.add_sum), itr::Vector{Any}, ItrEltype::Base.HasEltype
││││││││││││┌ reduce_empty_iter(op::Base.MappingRF{typeof(identity), typeof(Base.add_sum)}, itr::Vector{Any}, ::Base.Ha
│││││││││││││┌ reduce_empty(op::Base.MappingRF{typeof(identity), typeof(Base.add_sum)}, ::Type{Any}) @ Base ./reduc68 / 87
││││││││││││││┌ mapreduce empty(::typeof(identity) op::typeof(Base add sum) T::Type{Any}) @ Base /reduce jl:367

－34－

を⽤いた潜在的なエラーの発⾒ (4)

JET.jl

処⽅箋は次のと;．arr = Int[] とすればよい:
function sumevens(N::Integer)
N ≥ 1 || throw(DomainError(N, "`N` cannot be less than 1."))
arr = Int[] #
for n in 1:N
if iseven(n)
push!(arr, n)
end
end
return sum(arr)
end

この⾏を修正した

# sum(arr)

は⼤雑把にいえば下記のようなことを$
型の3H元

s = 0 # Int
for a in arr
s += a
end

Int

型の値を格納する配列に対し 空の配列の sum は 0 とするのが⾃然だから sum(arr) が 0 として計算できる

JET.jl

69 / 87

周り学習リソース(JuliaCon)
JET.jl

の紹介

レクチャー動画

JET.jl: The next generation of code checker for Julia | Shuh…
Shuh…

Package development: improving engineering quality & late…
late…

70 / 87

－35－

⾼速化周り学習リソース(JuliaCon)
佐藤さんによる Julia⾼速化の常識･⾮常識 @ Bio”Pack”athon2022#12
⾼速化の常識･⾮常識 @ Bio”Pack”athon2022#12

Julia

71 / 87

Cthulhu.jl
クトゥルフ と発⾳するらしい
ローカル変数がどのような型になっているかを観察できる
誤解を恐れずに⾔えば Debugger.jl の型バージョン
Making Abstract Interpretation Less Abstract in Cthulhu.jl | …

72 / 87

－36－

BenchmarkTools.jl
ベンチマークをとることができる
型不安定な実装
┊
┊
Time (mean ± σ): 3.540 ms ± 588.671 μs ┊ GC (mean ± σ): 4.80% ± 10.34%
▆▇█▅▄▂
███████▇▆▅▅▁▄▄▅▁▄▅▄▁▄▁▄▅▅▅▅▆▅▅▅▅▅▅▄▅▅▄▄▅▅▆▅▅▆▅▆▅▇▅▅▅▆▆▅▆▆▅

julia> @benchmark main1($N) #
BenchmarkTools.Trial: 1410 samples with 1 evaluation.
Range (min … max): 3.257 ms … 6.721 ms
GC (min … max): 0.00% … 39.38%
Time (median): 3.367 ms
GC (median): 0.00%

3.26 ms Histogram: log(frequency) by time 5.97 ms <
Memory estimate: 3.83 MiB, allocs estimate: 149983.
julia> @benchmark main3($N) #
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
Range (min … max): 215.560 μs … 10.235 ms
GC (min … max): 0.00% … 97.30%
Time (median): 249.598 μs
GC (median): 0.00%
Time (mean ± σ): 338.071 μs ± 647.998 μs
GC (mean ± σ): 24.88% ± 12.44%

型不安定な実装
┊
┊
┊
█▂
▁
██▆▅▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄▄▄▄▄▄▄▄▆▄▄▄▆▄▄

216 μs
Histogram: log(frequency) by time
4.59 ms <
Memory estimate: 800.88 KiB, allocs estimate: 14.

73 / 87

Proﬁle.jl/ProﬁleView.jl/ProﬁleSVG.jl
どの部分が時間がかかているかを調べることができる．
先ほどの Cthulhu.jl と連携することもできる
Tools for xing performance problems in Julia: Pro leView …

74 / 87

－37－

Replay.jl
での作業を⾃動化する．
ドキュメントはこちら
hyrodium さんありがとうございます
REPL

75 / 87

Pluto.jl

⼊⾨
76 / 87

－38－

Pluto.jl

について

が動作するノートブック
ソースコード: https://github.com/fonsp/Pluto.jl
使い⽅: https://github.com/fonsp/Pluto.jl/wiki
ユーザが書いコードが更新されるとセル間の依存関係を⾃動で解決
パッケージの依存関係が⼀つのファイル内に記録される
動作の再現がしやすい(=他者と共有しやすい)
PlutoUI.jl と連携し簡易 UI を構築するとがで る
頑張ると Introduction to Computational Thinking の9%ッ]%教育資源を構築するとがで る
Julia

77 / 87

Pluto.jl

の使い⽅

インストール
$ julia -e 'using Pkg; Pkg.add("Pluto")'

起動
デフォルトのポート番号は 1234
$ julia -e 'using Pluto; Pluto.run()'

Docker

の Julia 公式イメージから出発すると次の9&起動すると " る:

$ docker run --rm -it -v $PWD:/work -w /work \
-p 1234:1234 \
julia:1.9.2 \
julia -e 'using Pkg; Pkg.add("Pluto"); using Pluto; Pluto.run(host="0.0.0.0")'

78 / 87

－39－

Demo
ソースコードはこちら
Lévy C curve
反復関数系
RandomLogos.jl

79 / 87

⾖知識
ローカルにある⾃作パッケージを Pluto で動かす場合
julia --project=@. -e 'using Pluto; Pluto.run()' で Pluto を起動
下記のセルを追加
begin
Pkg.activate(Base.active_project())
end

terasakisatoshi/Sacabambaspis.jl
Pluto's built-in package management

詳しいことは

を読むと良い

80 / 87

－40－

Pluto.jl

周り学習リソース(JuliaCon)
Pluto.jl

の紹介

開発者らの最近のー動画

JuliaCon 2020 | Interactive notebooks ~ Pluto.jl | Fons van …

Julia for HPC webinar - Fons van der Plas - Pluto.jl - Reprod…
Reprod…

81 / 87

VS Code (1)
julia-vscode/julia-vscode
Julia in VS Code - What's New | David Anthoff, Sebastian P…
P…

82 / 87

－41－

VS Code (2)
VS Code 1.80.0

から VS Code の統合ターミナル上で画像を表⽰で <9&%:

using ImageInTerminal
using TestImages
m = testimage("mandril_color")
using FileIO, Sixel, Plots
gr()
buf = IOBuffer()
show(buf, MIME("image/png"), plot(sin, size=(500, 300)))
buf |> load |> sixel_encode
using OpenCV
using ImageCore: normedview, colorview, RGB
using ImageInTerminal
img_bgr = OpenCV.imread("sin.png")
img_rgb = OpenCV.cvtColor(img_bgr, OpenCV.COLOR_BGR2RGB)
# Julia
jlimg = colorview(RGB, normedview(img_rgb))
#
jli '

のエコシステムが理解できるデータ構造に変換する
メモリレイアウトの関係上， 空間⽅向の座標が反転していることに注意. 転置する.

83 / 87

VS Code (3)
によって部分的にテストを実⾏可能
機能 A はテスト 通った
機能 B のデバッグとテストに集中したい

@testitem

JuliaTesting/ReTestItems.jl
terasakisatoshi/Sacabambaspis.jl

84 / 87

－42－

まとめ
85 / 87

まとめ
⾔語の導⼊, REPL 周りのエコシステムを紹介した
ここで紹介し =なかったこと)いっ+いある
思ったよりも便利な機能があると思います
皆さんの感想を教えてください
エコシステムが個⼈や集団からなるコミュニティベースで作ら=ていることが?かる
Julia⾔語を多くの⼈に試してもらい作る側の⽴場に⽴てる⽅が増え=*幸いです
Julia

86 / 87

－43－

以上
87 / 87

－44－

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

型と多重ディスパッチ
数学と物理におけるJuliaの活用 チュートリアル講演(2)
2023/07/10 後藤俊介

本講義の予定
1. はじめに（5分）
2. 型システムの基本 (10分)
3. 複合型と抽象型 (10分)
4. 多重ディスパッチの基本（10分）
5. 実習：Julia の多重ディスパッチと型定義（10分）
6. まとめ（5分）

－45－

はじめに

講義の目的と目標

● 型（の基本）の理解
○
○
○
○

型とは何か？
Julia の型システム
型の定義
型パラメータ、型制約

● 多重ディスパッチ（の基本）の理解
○ 多重ディスパッチとは？
○ メソッドの選択の仕組み

－46－

実験コード
● 本講演資料に記載のサンプルコードを 以下のURLで公開しています。
https://antimon2.github.io/julia_imi_workshop2023_tutorial/1st_examp
le_type_md.html

● チュートリアル(1) で Pluto.jl を導入済の方は、上記URL→[Edit or Run]
ボタン→「On your computer」の手順に従って手元の環境で実行できます。
● 実習課題も編集可能な状態で用意してありますので
ご利用ください。

型システムの基本

－47－

型とは何か（1）

julia> typeof(1)
Int64

プログラミング全般において…

julia> typeof("文字列")
String

● その値（オブジェクト）
が『何』であるかを表
すもの。
● （Julia の場合）
typeof() 関数でそ
の値の型が分かる

julia> typeof([1.0, 2.0, 3.0])
Vector{Float64} (alias for Array{Float64, 1})

型とは何か（2）
少し概念的な例えとして…

● 物理で言えば『単位』
●
●

物理：『単位』を意識すること
で現象や状況を理解
プログラミング：変数や式の
値の『型』を意識することで
その式の意味や状況を理解

－48－

型とは何か（3）
さらに抽象的な例えとして…

● 数学（圏論）で言う
『対象』

画像は Wikipediaより引用
A, B, C が『対象』、f, g, g∘f は『射』（Julia における関数に対応）

型とは何か（4）
要するに…

● 値（オブジェクト）を分類（もしくはラベル付け）した
もの
● その値（オブジェクト）が『どのような性質を持つの
か』を表す
● その値に対して『どのような操作ができるのか』を
規定する

－49－

Julia の型システム（1）
Julia の型システムの基本…

●

●

●

公称型システム
(Nominative Type
System)
○ 名前で型が決まる！
公称的サブタイピング
(Nominal Subtyping)
○ 明示的な宣言で型の
基本型-派生型の関
係が決まる！
typeof() 関数、 <: 演算
子、isa 演算子等

julia> typeof(1) === Int
true

# Int という型は唯一

julia> Int <: Signed <: Integer <: Real <: Number <: Any
true

# ↑ `<:` は派生型演算子（左辺が右辺の派生型なら `true`）

julia> 1 isa Number
true

# `a isa T` は `typeof(a) <: T` と（ほぼ）等価

julia> supertype(Int)
Signed

Julia の型システム（2）
サブタイピングについて…

●

●

●

（直接の）基本型はただ1
つだけ
○ supertype() 関数
派生型はいくらでもあり
うる
○ subtypes() 関数
Any型
○ 全ての型の基本型

julia> supertype(String)
AbstractString
julia> subtypes(AbstractString)

# 環境によって結果は変る

6-element Vector{Any}:
Core.Compiler.LazyString
LazyString
String
SubString
SubstitutionString
Test.GenericString
julia> length(subtypes(Any))
587

－50－

# 環境によって結果は変る

Julia の型システム（3）
型階層…

●

●

型どうしのサブタイピング
関係を図にするときれい
な階層構造が出現（型階
層）
（Anyを根とした）型階層
の図を『型ツリー』とも呼
ぶ

画像は拙著 実践Julia入門 より引用

代表的な型の紹介（1）

julia> typeof(4294967296)
Int64

整数型

●

●

●
●

符号付き整数
○
Int8, Int16, Int32,
Int64, Int128
○
Int は Int64 のエイリアス
符号なし整数
○
UInt8, UInt16, UInt32,
UInt64, UInt128
○
UInt は UInt64 のエイリア
ス
多倍長整数（BigInt）
リテラル
○
10進表記なら符号付き、16
進表記（など）なら符号なし
○
桁数や大きさで型が自動的に
決まる

julia> typeof(9223372036854775808)
Int128
julia> typeof(170141183460469231731687303715884105728)
BigInt
julia> typeof(0x01)
UInt8
julia> typeof(0x0123456789abcdef)
UInt64

－51－

julia> typeof(1.0) === typeof(0.) === typeof(-.0) === Float64

代表的な型の紹介（2）

true

# よくある表記は大抵 Float64

julia> typeof(1e308) === typeof(0x193p-2) === Float64

浮動小数点数型

●
●
●

Float16, Float32,
Float64
BigFloat
リテラルは様々な書式が
ある

true

# `0x《16進表記》p《2の冪数》`という書式も OK

julia> typeof(3.14f0)
Float32

# `xxxfyy` と言う書式にすると Float32 になる

julia> Float16(3.1416)

# Float16 は明示的に指定

Float16(3.14)
julia> big(π)

# BigFloat の例

3.1415926535897932384626433832795028841971693993751058209749445
92307816406286198

julia> typeof(1//2)

代表的な型の紹介（3）

Rational{Int64}
julia> -12//20

その他の数値型

●
●

有理数
○ ○//△
複素数
○ ○ + △im

-3//5

# 結果は適宜約分される

julia> 3//0

# これもOK（`0//0` のみエラー）

1//0
julia> im

# 虚数単位（定数として定義済）

im
julia> typeof(1 + 2im)
Complex{Int64}
julia> typeof(0.0 + 1.0im)
ComplexF64 (alias for Complex{Float64})

－52－

代表的な型の紹介（4）

julia> typeof('a')
Char

文字型と文字列型

●
●

●

Julia では 文字型 と
文字列型 は別物
シングルクォーテーション
で括ったものは文字型
(Char)
ダブルクォーテーションで
括ったものは文字列型
(String)

julia> typeof("a")
String
julia> typeof("""
複数行にわたる文字列
これも String 型
""")
String

代表的な型の紹介（5）

julia> typeof([1.0, 2.0, 3.0])

配列型

# ↑リストではなくベクトル（＝1次元配列）

●

●

●

Julia では（同じ型の）値
の列はリストではなく（1
次元）配列
Julia の配列は1次元だ
けでなく 多次元配列 も
OK
要素の型と次元数は型パ
ラメータ（後述）で表現

Vector{Float64} (alias for Array{Float64, 1})

julia> typeof([
1 2
3 4
])

# 行列（＝2次元配列）

Matrix{Int64} (alias for Array{Int64, 2})
julia> typeof([1;2;;3;4;;;5;6;;7;8])
Array{Int64, 3}

－53－

# 2×2×2 の3次元配列

julia> Dict("Alice"=>1, "Bob"=>2, "Carol"=>3)

代表的な型の紹介（6）

Dict{String, Int64} with 3 entries:
"Carol" => 3

辞書型・集合型

●

●
●

辞書型(Dict)：キーと値
のマッピングで管理するコ
レクション
集合型(Set)：重複を許さ
ない値のコレクション
配列型、辞書型、集合型と
後述のタプル・名前付きタ
プルを総称して コレク
ション型 と言う

"Alice" => 1
"Bob"

=> 2

# 結果は順不同

julia> Set([3, 1, 4, 1, 5, 9, 2, 6, 5, 3])
Set{Int64} with 7 elements:
5

# 重複は排除される

4
6
2
9
3
1

代表的な型の紹介（7）

# こちらも順不同

julia> typeof((1, 'b', "三"))
Tuple{Int64, Char, String}

タプル型・名前付きタプル型

●
●

●

タプル型(Tuple)：（同じ
型とは限らない）値の列
名前付きタプル型
(NamedTuple)：タプル
の各値に名前（キー）が紐
付いているもの（辞書に類
似）
これらは、各値の型も情報
として保持するのが特徴

julia> typeof((a=1, b='b', c="三"))
NamedTuple{(:a, :b, :c), Tuple{Int64, Char, String}}

－54－

julia> let

型アノテーション

x::Int = sin(π)

# == 0.0 なので変換されて 0 が代入される

x

●
●

●

●

ｘ：：T という書式で 型ア
ノテーション を付けられ
る
宣言時や代入時の左辺の
場合は（変数の）型の指定
（変換に失敗するとエ
ラー）
代入時の右辺などの場合
は 型検査 （互換性のない
型だとエラー）
※型ヒントではない！
（重要）

end
0
julia> let
x = (sin(π)::Int)

# こちらはエラー（!(sin(π) <: Int) なので）

x
end
ERROR: TypeError: in typeassert, expected Int64, got a value of type Float64
Stacktrace:
[1] top-level scope
@ REPL[XX]:2

複合型と抽象型

－55－

複合型(1)

julia> struct SSample
x

（構造体）
●
●
●

struct 《型名》 ～ end
で 構造体 を定義可能
構造体として定義された
型のことを 複合型 と呼ぶ
※基本的な型以外のほと
んどの型は複合型
（有理数型・複素数型、辞
書型・集合型・名前付きタ
プル型も実は複合型）

複合型(2)

y
end
julia> ssample = SSample(1, 2)
SSample(1, 2)
julia> typeof(ssample)
SSample

julia> struct SISample
x::Int

（フィールドの型指定）
●

●

●

構造体のフィールドには
型が指定できる（型アノ
テーション）
※型アノテーションが付い
ていないフィールドは実
は ::Any と同じ意味（詳
細略）
補足：目的がはっきりして
いるときはそのフィールド
には型アノテーションを付
けましょう（詳細後述）

y::Int
end
julia> sisample = SISample(1, 2)
SISample(1, 2)
julia> typeof(sisample.x)
Int64
julia> SISample(2.0, 3.14)

# `3.14` の方でエラーが発生

ERROR: InexactError: Int64(3.14)
# :《以下略》

－56－

複合型(3)

julia> struct MyDecimal <: Real
value::BigInt

（基本型の指定）
●

●

型定義時に基本型を指定
（その型の派生型という宣
言）ができる
※基本型の指定がない型
定義は、<: Any（Any型
の派生型の定義）と同じ意
味（詳細略）

point::Int
end

# 固定小数点数を意図した型定義

julia> MyDecimal(1, 0) isa Number
true
julia> MyDecimal(1, 0) isa Integer
false

julia> abstract type AbstractFPoint end

抽象型(1)
●

●

abstract type 《型名》
end で 抽象型 を定義で
きる
抽象型は他の抽象型や複
合型等の基本型にできる
（＝型階層が作れる）

julia> struct FPoint2D <: AbstractFPoint
x::Float64
y::Float64
end
julia> struct FPoint3D <: AbstractFPoint
x::Float64
y::Float64
z::Float64
end
julia> FPoint2D(1.0, 2.0) isa AbstractFPoint
true
julia> FPoint3D(3, π, 99.9) isa AbstractFPoint
true

－57－

julia> struct TWrapperSample{T}

型パラメータ(1)
●

●

型定義時に {T} のように
型パラメータ を指定でき
る
型パラメータの使い途：
○ フィールドの型指定
○ 性質を表すパラメー
タ（フラグ）（詳細後
述）

value::T
end
julia> intwrapper = TWrapperSample(1)
TWrapperSample{Int64}(1)
julia> typeof(intwrapper)
TWrapperSample{Int64}
julia> typeof(TWrapperSample("文字列"))
TWrapperSample{String}

julia> abstract type AbstractPoint{T <: Real} end

型制約(1)
julia> struct Point2D{T} <: AbstractPoint{T}

●

●

型パラメータ指定時に {T
<: Real} のようにその
型に制約を設けることが
できる（型制約）
制約に合わない型が来る
とエラーになる

x::T
y::T
end
julia> Point2D(1.0, 3.2) isa AbstractPoint{Float64}
true
julia> Point2D(1 + 0im, 0 + im)
ERROR: TypeError: in AbstractPoint, in T, expected
T<:Real, got Type{Complex{Int64}}
# ：《以下略》

－58－

我々は何を見せられているのか…

型パラメータ(2)

julia> abstract type AbstractPoint{T <: Real} end

＋型制約(2)
julia> struct Point2D{T} <: AbstractPoint{T}

この形
→
は非常によく見るので、
ぜひ覚えておきましょう！

x::T
y::T
end

利点：
●
●
●

フィールドの型を特定しな
いので柔軟な設計が可能
実行時には確定するので
型安定性 に繋がる
型制約と組み合わせると
さらに賢い設計に！

－59－

多重ディスパッチの基本

julia> add(x, y) = x + y

多重ディスパッチとは(1)
関数の多重定義(1)

●

●
●

Julia では同名の関数を
引数の違いで 多重定義
できる
個々の実体のことを メ
ソッド と呼ぶ
引数の違いとは
○ 引数の個数の違い
○ 引数の型の違い
○ その組み合わせ

add (generic function with 1 method)
julia> add(x, y, z) = x + y + z
add (generic function with 2 methods)
julia> add(x, y, z...) = add(x + y, z...)
add (generic function with 3 methods)
julia> methods(add)
# 3 methods for generic function "add" from Main:
[1] add(x, y)
[2] add(x, y, z)
[3] add(x, y, z...)

－60－

julia> double(x) = 2x

多重ディスパッチとは(2)
関数の多重定義(2)

●

●
●

Julia では同名の関数を
引数の違いで多重定義で
きる
個々の実体のことをメ
ソッドと呼ぶ
引数の違いとは
○ 引数の個数の違い
○ 引数の型の違い
○ その組み合わせ

double (generic function with 1 method)
julia> double(s::AbstractString) = s ^ 2
double (generic function with 2 methods)
julia> double(x, y) = string(double(x), double(y))
double (generic function with 3 methods)
julia> double(x::Number, y::Number) = double(x) + double(y)
double (generic function with 4 methods)
julia> methods(double)
# 4 methods for generic function "double" from Main:
# ：《略》

julia> double("文字列")
# "文字列" isa AbstractString なので s ^ 2 の実装（メソッド）が選択される

多重ディスパッチとは(3)

"文字列文字列"

ディスパッチ

julia> double(1, 2)
# 1 isa Number かつ 2 isa Number なので double(x) + double(y)

●

多重定義された関数は、
呼び出し時の実引数の組
み合わせで適切なメソッド
が選択されて実行（＝ディ
スパッチ）される

6
julia> double(160, "円")
# どうしてこうなるのか考えてみよう！
"320円円"

－61－

そもそも関数とは？（1）
抽象的な例えとして…

● 数学（圏論）で言う
『射（射影）』

画像は Wikipediaより引用 （再掲）
A, B, C が『対象』、f, g, g∘f は『射』（Julia における関数に対応）

そもそも関数とは？（2）
少し概念的な例えとして…

● 物理で言えば『式』
●

●

物理：似たような『式』でも
『単位』を意識することで表
される現象や得られる結果
が変わる
Julia：『関数』に渡される
『型』によってその挙動の意
味や得られる結果が変わる

－62－

多重ディスパッチとは(4)
重要な概念として…

● 関数の名前を固定することで、その関数に『役割』を与える。
○ 物理学において「だいたい同じ形の式で書けるものは似たような現
象を表す」ことにだいたい相当

● その関数を型（の組み合わせや個数）の違いで多重定義（メ
ソッドを定義）することで、特定の状況に置ける意味とその
結果を与える（実装する）。
○ 同じ式でも『単位』の違いで状況・意味・結果が変ることに類似

Juliaにおける多重ディス
パッチの利用例(1)

julia> log(20)

引数の個数による意味分け

julia> log(5, 20)

●
●
●

同じ関数で、引数の個数
によって意味を分ける
例：log() 関数
例：atan() 関数

# 自然対数

2.995732273553991
# 底を5とする対数

1.8613531161467862
julia> atan(0.3)

# tan(θ) == 0.3 となる θ

0.2914567944778671
julia> atan(1, 3)
0.3217505543966422

－63－

# tan(θ) == 1/3 となる θ

Juliaにおける多重ディス
パッチの利用例(2)

julia> 1 * 2

演算子オーバーロード

julia> 3.0 * π

●

●

同じ演算子でも両辺（被演
算子）の型で挙動（結果の
値や型）に違いが出る
わかりやすい例：* 演算子
○ 数値同士なら乗算
○ 文字列同士なら結合
○ 数値同士でも、型の
違いによって適切な
処理

# 整数どうしの `*` 演算は結果も整数

2
# 浮動小数点数と無理数の `*` 演算は浮動小数点数

9.42477796076938
julia> "Hello, " * "Julia!"

# 文字列同士の `*` 演算は結合

"Hello, Julia!"

Juliaにおける多重ディスパッチの利用例(3)
その他の例…

● Holyトレイト
○ 特殊な型による定義分けを利用して、関係ない型での共通の挙動（振る舞
い）の定義（実装）を共有化するテクニック

● アルゴリズムの分離・モジュール化
○ ソートアルゴリズムの定義など

● 擬似値ディスパッチ
● ポリモーフィズム

－64－

実習：Julia の多重ディスパッチと
型定義

お題：double()関数(1)
Step1: 型定義

●

●

MyType という型（複合
型）を定義する
○ フィールドは value
1つだけ
○ フィールドの型は型
パラメータ T
※これに関しては答えを
示しておきます→

julia> struct MyType{T}
value::T
end
julia> MyType(1)
MyType{Int64}(1)

－65－

julia> struct MyType{T}

お題：double()関数(2)

value::T
end

Step2: 関数の多重定義(1)

●

●

先ほど例示した
double() 関数を、
MyType 型が扱えるよう
に多重定義する。
まずは引数に
MyType{T} 1つだけを
受け取るメソッドの定義：
○ フィールド value
に double() 関数
を適用した結果を格
納する Mytype を
返す

julia> MyType(1)
MyType{Int64}(1)
julia> double(mytype::MyType) = ...

# ここを適切に実装

double (generic function with 5 methods)
julia> double(MyType(1))
MyType{Int64}(2)
julia> double(MyType("ABC"))
MyType{String}("ABCABC")

julia> double(x::MyType, y::MyType) = ...

# ここを適切に実装

double (generic function with 6 methods)

お題：double()関数(3)
julia> double(x::MyType, y) = ...

Step3: 関数の多重定義(2)

●
●

double() 関数を、
MyType 型が扱えるように
多重定義する。
引数を2つ（一方が
MyType） を受け取るメ
ソッドの定義（難易度高め）：
○ 結果は MyType、た
だし
MyType{MyType{
...}} のような入れ子
にならないように注
意！

# ここを適切に実装

double (generic function with 7 methods)
julia> double(x, y::MyType) = ...

# ここを適切に実装

double (generic function with 8 methods)
julia> double(2, MyType(π))
MyType{Float64}(10.283185307179586)
julia> double(MyType(35), "v")
MyType{String}("70vv")
julia> double(MyType("ABC"), MyType("XYZ"))
MyType{String}("ABCABCXYZXYZ")

－66－

julia> double(x::Char) = ...

# ここを適切に実装

double (generic function with ? methods)

お題：double()関数(4)
julia> double(x::Tuple) = ...

Extra Step（時間が余った人用）

●

●

double() 関数を、他の型
にも適用しよう
○ 例：文字型⇒文字列と
して扱う
○ 例：コレクション型⇒各
要素にdouble()を適
用したコレクションに
する
○ 他⇒自分で仕様を考
えて実装してみよう
型パラメータの利用
○ 型パラメータの違いで
処理を分ける、など

# ここを適切に実装

double (generic function with ? methods)
julia> double(x::Array) = ...

# ここを適切に実装

double (generic function with ? methods)
julia> double((1, `b`, [1.0, 2.0, 3.0]))
(2, “bb”, [2.0, 4.0, 6.0])
julia> double(x::MuType{T}) where {T <: Tuple} = ...
double (generic function with ? methods)
julia> double(Mytype((1, `b`, “三”)))
#> MyType((1, `b`, “三”, 1, `b`, “三”))

まとめ

－67－

# ←一例

# ここを適切に実装

講義の目的と目標（再掲）

✓ 型（の基本）の理解
☑
☑
☑
☑

型とは何か？
Julia の型システム
型の定義
型パラメータ、型制約

✓ 多重ディスパッチ（の基本）の理解
☑ 多重ディスパッチとは？
☑ メソッドの選択の仕組み

お疲れさまでした

－68－

̢͕̬̲༊Ⴓˡ‫ܤ‬஽
ଲ‫˛ڮ‬໲ၟ˞ʽ˄˾Ɋɽɴɮɥˡ‫ဳھ‬ȝ̪͎ʏ̱͓͔̋࣠‫؃‬ȧȲȩ
ȱȯȱȲȮȯȶȮȰȯȝ࢑൳੕٨

ྂ࣠‫ݓ‬ˡဟ೾
Ȱȭ ˢˋ˴˞Șȴ໴ș
ȱȭ ̢͕̬̲༊Ⴓ˛ˢ‫ط‬ʾȧȰȯ໴ȩ
Ȳȭ Ɋɽɴɮɥ˞ʽ˄˾̢͕̬̲༊Ⴓˡ‫ྂܤ‬ȧȰȯ໴ȩ
ȳȭ ০ႣȘȰȯ໴ș
ȴȭ ০ਮȦ༊Ⴓ‫̺̒̾˾˻˞ص‬ʏ̢͇͜‫ٶ‬஫ȘȰȯ໴ș
ȵȭ ಣ֑ല˛˱˛˴Șȴ໴ș

－69－

ˢˋ˴˞

࣠‫ݓ‬ˡ࿠ജ˛࿠ຢ

● ̢͕̬̲༊ႳȘˡ‫ྂܤ‬șˡၟ٪
○
○
○
○

༊Ⴓ੪ၟ˛ˢ‫ط‬ʾ
̢͕̬̲ˡ‫ྂܤ‬ജ˝‫ڋ‬ෞ
̢͕̬̲ˡၚല˛୅࡮
Ș̢͕̬̲ˡ̞̘͔͒̍̾̍ș

● Ɋɽɴɮɥ˞˻˾̢͕̬̲༊Ⴓˡၟ٪
○ य೷࡮ˡ়ப੢຃
○ ‫ڜ‬ਓɀɓɉ
○ ̢͕̬̲̤ʏ̾Ș˛̰̬̲͖̬̘ș

－70－

০࡟̜ʏ̲
● ྂ࣠‫؃‬ষႁ˞ܿॄˡ̞͔̜̀͜ʏ̲̅ȝֆ‫ش‬ˡɘɕɎ˚ࢤ‫˙ˊځ‬ʷ˱ˌʇ
ɬɼɼɸɻȹȮȮɥɶɼɮɵɷɶȱȭɫɮɼɬɽɦȭɮɷȮɯɽɴɮɥɢɮɵɮɢɿɷɺɱɻɬɷɸȱȯȱȲɢɼɽɼɷɺɮɥɴȮȰɻɼɢɩʀɥɵɸ
ɴɩɢɶɼɬɺɩɥɨɻȭɬɼɵɴ

● ̪͎ʏ̱͓͔̋ȧȰȩȝ˚ȝɓɴɽɼɷȭɯɴȝ̅ඃ෎हˡཅˢʆૅܿɘɕɎ@ɞɄɨɮɼȝɷɺȝɕɽɶɠ
̨͜ͅ@ʚɒɶȝʁɷɽɺȝɧɷɵɸɽɼɩɺʛˡ਎੩˞ੀ˖˙਎ࡡˡ۸‫˚ޝ‬০ࣞ˚ˀ˱ˌʇ
● ০ਮٕూ˵༜਺‫ؼ‬෫˝૕మ˚ဳ֑ˊ˙ʵ˽˱ˌˡ˚
ˇၚဳ˂˓ˈʷʇ

̢͕̬̲༊Ⴓ˛ˢ‫ط‬ʾ

－71－

༊Ⴓ੪ၟ˛ˢ‫ط‬ʾȘȰș
͖̙͈̙̀͒͜ம๘˞ʽʷ˙

● 並行計算(Concurrent Computing)Ȧ
໫ଲˡࠥঁ̅ඁৄ˞০ࣞˌ˾ˆ˛
● 並列計算(Parallel Computing)Ȧ
細分化された独立な処理ȝ̅༊ࣞࠥঁˌ˾ˆ˛
● Șྂ࣠˚ˢș༊Ⴓࠥঁ੪ၟˡˆ˛̅ঞˊ˙
並列処理ȝ˛ࡲ˩ˆ˛˞ˌ˾

༊Ⴓ੪ၟ˛ˢ‫ط‬ʾȘȱș
˝ˏຖ္ʾ˛ʷʹʽს

● ిܾ࿓˝ࠥঁ̨̢̘̅ࢦၨജ˞੪ၟˌ˾˒˴ˡ਎ಀ
Ș‫ځ‬฻ʈࡒ‫̢̽ݸ‬ʏ̲̅ૅ˅˾˒˴ș
● ͓̦ʏ̢ˡࢦၨജ˝ၚဳ
Șɂɓɘ̜̋̅မˣˎˍ͔̾‫˴˒˾ˌဳھ‬ș
● ̨͓͔͉͓̋̍̋̀̚ʏ̠͐͜ʾ˼ˡਟ္
Ș໫ଲˡ̨̢̘̅ඁৄ˞କఆ˞੪ၟˌ˾˒˴ș

－72－

julia> Threads.nthreads()

༊Ⴓ੪ၟ˛ˢ‫ط‬ʾȘȲș
̋ʏ̖̯̘̪͌ജ˝ს

● ֨Ⴠˡ੪ၟȘ̨̢̘ș
̅ඁৄ˞௰˼ˎ˙੪
ၟ̅ˌ˾ˆ˛
● ࿠ജˢȝスループット
向上Ș̺̾̒ʏ̢͇͜ࢪૅș
● ʷ˂˗ʾਓႠʿʵ˾Ȧ
○ スレッド並列
○ プロセス並列

6
julia> fib(n) = n ≤ 1 ? n : fib(n - 2) + fib(n - 1)
fib (generic function with 1 method)
julia> @time fib(40)
0.510351 seconds
102334155
julia> @time for _=1:4 fib(40) end
2.073198 seconds (8.14 k allocations: 600.650 KiB, 1.68% compilation time)
julia> @time Threads.@threads for _=1:4 fib(40) end
0.783171 seconds (34.33 k allocations: 2.571 MiB, 19.96% compilation time)

̢͕̬̲༊Ⴓ˛ˢȘȰș
˱ˍˢဳ࢘ୌၟ

● ̨̢̘Ȧ
֨Ⴠˡ੪ၟˡ˱˛˱
˽Șྂ࣠˚ˢș
● ̢͕̬̲Ȧ
ɂɓɘˡय઎੪ၟ౧ֈ
● ̢͖̤̀Ȧ
͖̙͉̀͒ˡ০ࣞ౧
ֈ

－73－

̢͕̬̲༊Ⴓ˛ˢȘȱș

● ೠ૏ˡ͖̙͉̀͒˚ˢʆ౧̢͖̤֨̀ʈ౧̢͕̬̲֨˚
Ȱ˗ˡ̨̢̘ʿ඀˂˓˄Ș̢̠̙͔͖̤̀͜ʈ̢̠̙͔͜
͕̬̲ș
● ౧̢͖̤֨̀ಜ˚໫ଲˡ̢͕̬̲̅ඁৄ˞‫۽‬ၟˊʆ໫ଲ
ˡ̨̢̘̅ඁৄ০ࣞȘr༊Ⴓ੪ၟș˚ˀ˾Ș̠̙͔͖̀͜
̢̤ʈマルチスレッドș
● ໫ଲ̢͖̤̀̅ඁৄ˞ၩ˔ૅ˅˙໫ଲˡ̨̢̘̅ඁৄ
০ࣞȘr༊Ⴓ੪ၟșˌ˾ˆ˛˵˚ˀ˾Șマルチプロセスș

̢͕̬̲༊Ⴓ˛ˢȘȲș

● ೠ૏ˡ͖̙͉̀͒˚ˢʆ౧̢͖̤֨̀ʈ౧̢͕̬̲֨˚
Ȱ˗ˡ̨̢̘ʿ඀˂˓˄Ș̢̠̙͔͖̤̀͜ʈ̢̠̙͔͜
͕̬̲ș
● 単一プロセス中で複数のスレッドを同時に管理し、複数
のタスクを同時実行（≒並列処理）できる（シングルプ
ロセス・マルチスレッド）⇒スレッド並列
● ໫ଲ̢͖̤̀̅ඁৄ˞ၩ˔ૅ˅˙໫ଲˡ̨̢̘̅ඁৄ
০ࣞȘr༊Ⴓ੪ၟșˌ˾ˆ˛˵˚ˀ˾Șマルチプロセスș

－74－

̢͕̬̲༊Ⴓˡၚല˛୅࡮

● 5໫ଲ̨̢̘ඁৄ০ࣞ˞˻˾̢͔ʏ̬̱̀ࢪૅȘ̺
̾̒ʏ̢͇͜ࢪૅșȘr࣬ఆ‫ص‬ș
● 3౧̢͖̤֨̀˝ˡ˚‫۽‬ၟʿ๲‫ک‬ജ‫ڰ‬Șɾȭɻȭ̢͇͔̪͖̤̀ș
● :౧̢͖̤֨̀˝ˡ˚͓͊͋ˢ‫ޖ‬ဏȘ̢͖̤̀‫܊‬ೠ૬˞
˻˾ిˀ˝̓ʏ̹ʏ̬̲́ˢ˝ʷș
● ;Șːˡఽ̂˽ș͓͊͋‫˞ࣱޕ‬ˢಣ֑Ș̢͕̬̲̤ʏ֑̾
৕ˎ˻ʆ઼ु࢑੔ș
● ;̰̬̲͖̬̘ಣ֑Ș‫࢑္ڋ‬੔ș

༬అȦʜ૽ˡ̢͕̬̲༊Ⴓʝ
ৄ‫܊‬ʿʵ˿ˣȘ˝˄˿ˣ຀ˣˌˡ˚඙̆˚ʽʷ˙˂˓ˈʷș

● ఖ࡬࢘˞˵̢͕̬̲Ș̢͇͔̪͕̬̲ș˛ʷʹܳ෫ˢʵ˾
● ֨ໞ࡬࢘ˡ̢͕̬̲ˢʆɆɉɎȘɆɴɷɦɥɴȝɉɶɼɩɺɸɺɩɼɩɺȝɎɷɧɱș˛ʷʹܳࣇ˞˻˽
༊Ⴓ඀ॖʿ୅࡮ˈ˿˙ʷ˾
● Ɋɽɴɮɥȝ˞ˢːˡ୅࡮ʿ˝ʷȘʜ૽ˡ̢͕̬̲༊Ⴓʝș
● ఖ࡬࢘˚ʚ༊Ⴓ੪ၟˌ˾˝˼̢͇͔̪͖̤̀ʛ˛ˊ˙ഇւʈഇઇˈ˿˙ʷ˾়ւ
˚˵ʆɊɽɴɮɥȝ˝˼̢͕̬̲༊Ⴓ˚‫ʹࣱ˞܊‬Ș˳ˊ̀ː˔˼ˡཅʿ਎˖਌˽௛ʷș˵
ˡ˵గʷ
● Ɋɽɴɮɥȝ˚ˢʚ̓ʏ̹ʏ̬̲́˞˻˾̰͓̬̱͊˻˽˵ʜ̢͕̬̲̤̋͜ʏ̾˷̰̬̲
͖̬̘ˡ૲ฉ̅ˊ˝˂˙ह˳ʝ˛ʷʹ͓̬̱͊ˡཅʿઁ˾ʛ˻ʹ˝ૌࣱ˞ˢȝ͖̤̀
̢༊Ⴓȝ˵தొ੒၊˾Șதొমʿగʷș˛ࡈ୻ജ˞ࣗʻ˾ˆ˛˵˚ˀ˾

－75－

༬అȦ̢͕̬̲ˡ̞̘͔͒̍̾̍
ˈ˼˞‫ڋ‬ෞജʾ˗඲ໞঋုജ˝სȘৄ‫܊‬ʿ˝˄˿ˣ຀ˣˌˡ˚඙̆˚ʽʷ˙˂˓ˈʷș

● ֨๘ജ˝̢͕̬̲ˡ૕మȘ̢̯ʏ̡șȦʚȰȭȝ૵ܾȧɑɩɿȩʛʚȱȭȝ০ࣞ‫ؼ‬෫
ȧɕɽɶɶɥɦɴɩȩʛʚȲȭȝ০ࣞಜȧɕɽɶɶɮɶɫȩʛʚȳȭȝ͖̬̘̿ȧɁɴɷɧɱɩɨȩʛʚȴȭȝਬၻ
ȧɗɩɺɵɮɶɥɼɩɨȩʛ
● ೷͕͔͂˚̢͕̬̲̅ճʹຖ္ˡʵ˾۸‫˚ޝ‬ˢʆˆ˿˼ˡ૕మȘ̢̯ʏ̡ș̅ഞ୶
˞ෲժʈ‫۽‬ၟˌ˾ຖ္ʿʵ˾
● Ɋɽɴɮɥȝ˚ˢʵ˱˽ૺ˂ࣗʻ˾ຖ္ˢ˝ʷ
● ຖ္˞˝˖˒˛ˀ˞ȝɊɽɴɮɥȝʿ඲ໞ˚̢͕̬̲̅ၩ˔ૅ˅ȘʚȰȭȝɑɩɿʛș̨̢̘̅
௰˼ˎ˙ȘʚȱȭȝɕɽɶɶɥɦɴɩʛʚȲȭȝɕɽɶɶɮɶɫʛșഞ‫ޘ݋‬ೈ඀ॖȘʚȳȭȝɁɴɷɧɱɩɨʛșˊ
˝ʿ˼̨̢̘ʿਬၻˊ˒˼ਬၻȘʚȴȭȝɗɩɺɵɮɶɥɼɩʛșˊ˙˂˿˾
● ˡ˚ˆˡ૕మȘ̢̯ʏ̡șˢʆ‫ڋ‬ෞ˛ˊ˙ॸࣗ˱˚˞

Ɋɽɴɮɥ˞ʽ˄˾̢͕̬̲༊Ⴓˡ
‫ྂܤ‬

－76－

य೷࡮ˡ়ப੢຃

> julia -t 4

# 最大スレッド数は4

> julia -t auto

Ɋɽɴɮɥȝ݁඀ৄ˞यి̢͕̬
̲ଲ̅ঞ೾ˌ˾
●
●
●
●

ɯɽɴɮɥȝȬɼȝȳ
Ș஀నଲঞ೾ș
ɯɽɴɮɥȝȬɼȝɥɽɼɷ
Șၚဳ˚ˀ˾ɂɓɘ̜̋य
ి‫ڤ‬༥ș
۸‫༙ޝ‬ଲ˞˻˾ঞ೾˵‫ؼ‬
෫
ɓɴɽɼɷȝˡૌࣱ̰͔̱̾̒
˚ʚၚဳ˚ˀ˾यిェȱʛ
ʿ୻೾ˈ˿˾

# 確保できる分だけスレッド数を確保

> env JULIA_NUM_THREADS=XXX julia

# Bash/Zsh などで環境変数指定と同時にJulia起動

# Julia 内では ↓ で利用可能なスレッド数が取得できる
julia> Threads.nthreads()
6

‫ڜ‬ਓɀɓɉȘȰș
ȿɼɬɺɩɥɨɻ

● ̢͕̬̲༊Ⴓˡ‫ྂܤ‬
● ɪɷɺȝⳒȝɩɶɨȝ̅‫۾‬౧
˞̢͇͔̪͕̬̲‫˚ص‬
ˀ˾
● ྂ࣠˚ˢ̡͎͋ʏ͔
̢͕̬̘̀̾̌຿ˀ˚ȝ
Threads.@threadsȝ
˛ʷʹວܿ̅ၚဳȘֆ
‫ش‬ඁုș

julia> Threads.nthreads()
6
julia> fib(n) = n ≤ 1 ? n : fib(n - 2) + fib(n - 1)
fib (generic function with 1 method)
julia> @time fib(40)
0.510351 seconds
102334155
julia> @time for _=1:4
fib(40)
end
2.073198 seconds (8.14 k allocations: 600.650 KiB, 1.68% compilation time)
julia> @time Threads.@threads for _=1:4
fib(40)
end
0.783171 seconds (34.33 k allocations: 2.571 MiB, 19.96% compilation time)

－77－

‫ڜ‬ਓɀɓɉȘȱș
ȿɻɸɥɿɶ

● ̢͕̬̲ˡ݁඀
Șːˡെോྂ൨˞૵ˊʷ̢͕̬̲
ʿ݁඀ˌ˾̂˄˚ˢ˝ʷʆ઼ु
‫ں‬բș

● ఋ˂৔̅Ȱ˗ˡ̨̢̘
˛ˊ˙˱˛˴ʆȘߥʷ
˙ʷ˾ș̢͕̬̲˚੪
ၟˌ˾
● ȿɻʁɶɧȝ͇̘͖˚ඁ
ܰ̅਌˾ຖ္ʵ˽

julia> Threads.nthreads()
6
julia> @time Threads.@threads for _=1:4
fib(40)
end
0.783171 seconds (34.33 k allocations: 2.571 MiB, 19.96% compilation time)
julia> @time @sync for _=1:4
Threads.@spawn fib(40)
end
0.715108 seconds (2.01 k allocations: 135.894 KiB, 8.50% compilation time)

julia> let
hist = zeros(Int, Threads.nthreads())
@time Threads.@threads for _=1:(10000*Threads.nthreads())
fib(15)

‫ڜ‬ਓɀɓɉȘȲș

# 毎回ほぼ同じ負荷の処理が実行される例

hist[Threads.threadid()] += 1
end
hist

ȿɼɬɺɩɥɨɻȝɾȭɻȭȝȿɻɸɥɿɶȝȧȰȩ

● ిႌˡुʾʷȘၳോʿ
ˮ˯ఏ˖˙ʷ˾ș̨̢
̘̅༊Ⴓ੪ၟˊ˒ʷ
ૌࣱˢȝȿɼɬɺɩɥɨɻȝ
ˡཅʿဇਪ

end
0.134166 seconds (22.88 k allocations: 1.588 MiB, 169.92% compilation time)
#> [10000, 10000, 10000, 10000, 10000, 10000]
julia> let
hist = zeros(Int, Threads.nthreads())
@time @sync for _=1:(10000*Threads.nthreads())
Threads.@spawn begin
fib(15)

# 毎回ほぼ同じ負荷の処理が実行される例

hist[Threads.threadid()] += 1
end
end
hist
end
0.301360 seconds (569.85 k allocations: 35.267 MiB, 102.04% compilation time)
#> [0, 12151, 12635, 11490, 11973, 11751]

－78－

julia> let
hist = zeros(Int, Threads.nthreads())
@time Threads.@threads for _=1:(250*Threads.nthreads())
fib(rand(25:35))

‫ڜ‬ਓɀɓɉȘȳș

# 毎回わずかに異なる負荷の処理が実行

hist[Threads.threadid()] += 1
end
hist

ȿɼɬɺɩɥɨɻȝɾȭɻȭȝȿɻɸɥɿɶȝȧȱȩ

● ̨̢̘ˡଲʿːˆ˱˚
ిႌ˚˝˂ʆၳോʿ̹
̹͒͒ˡૌࣱˢʆ
ȿɻɸɥɿɶˡཅʿ̺࣬
̾̒ʏ̢͇͜˛˝˾ˆ
˛˵ʵ˾
Ș໓ّ໴ॼș

end
4.370164 seconds (34.68 k allocations: 2.347 MiB, 6.02% compilation time)
#> [250, 250, 250, 250, 250, 250]
julia> let
hist = zeros(Int, Threads.nthreads())
@time @sync for _=1:(250*Threads.nthreads())
Threads.@spawn begin
fib(rand(25:35))

# 毎回わずかに異なる負荷の処理が実行

hist[Threads.threadid()] += 1
end
end
hist
end
4.310160 seconds (19.18 k allocations: 1.329 MiB, 1.33% compilation time)
#> [245, 244, 224, 283, 271, 233]

‫ڜ‬ਓɀɓɉȘȴș
ːˡఖ

● ɂɬɥɶɶɩɴȧɻɸɥɿɶȼɼɺɽɩȩȝɨɷȝⳒȝɩɶɨ
○
○

ඓ೾ˡ੪ၟ˓˄̬̅͒̀ˊ˙༔̢͕̬̲˚੪ၟ˚ˀ˾
ഞ‫݋‬ಃˡਙ˄ുˊʿ˚ˀ˾
Șࡡ ɂɬɥɶɶɩɴˢ̨̢̘‫˚܊‬ಃ̅ਙ˄ുˊˊ˝ʿ˼‫ޘ‬ೈ඀ॖˌ˾˒˴ˡ˵ˡș

● ɗɬɺɩɥɨɻȭɪɷɺɩɥɧɬȧȩ
○
○
○

Ɂɥɻɩȭɪɷɺɩɥɧɬȧȩȝˡ̢͇͔̪͕̬̲๓
੪ၟȘ‫܌‬ଲș˛̯͕̍ʏ̨Șɂɬɥɶɶɩɴ˚̬͒̀ˊ˒˵ˡș̅ਙ˄਌˽̢͇͔̪͕̬̲˚੪ၟ
໓ّ໴ॼ˞ဏဳ

● >ˆ˿˼ˡ઼ुˢ‫ں‬բ

－79－

০Ⴃ

julia> function my_matmul(A::AbstractMatrix, B::AbstractMatrix)
T = promote_type(eltype(A), eltype(B))

০ႣȘȰș

C = Matrix{T}(undef, (size(A, 1), size(B, 2)))
Threads.@threads for x = axes(B, 2)
Threads.@threads for y = axes(A, 1)

ȿɼɬɺɩɥɨɻȝˡႣȧȰȩ

● ฉႳˡ֗˝˾̰̬̍͜
̢̘˞ඁုˡ‫̅ঁ؃‬
০ঢˊ˙ಃ̅‫ڠ‬෪ˌ
˾̨̺ʏ͜Șˆˡ̺
̨ʏ͜ʿగʷș
● ȿɼɬɺɩɥɨɻȝˢ෎˿
ঘ˞˚ˀ˾ȘɊɽɴɮɥȝ
ɾȰȭȷȝֆࣩș

C[y, x] = @view(A[y, :])' * @view(B[:, x])
end
end
C
end
my_matmul (generic function with 1 method)
julia> let A=[1 2; 3 4; 5 6; 7 8], B=[1 2 3; 4 5 6]
C = my_matmul(A, B)
@assert C == A * B
@show C;
end
C = [9 12 15; 19 26 33; 29 40 51; 39 54 69]

－80－

julia> function my_matmul_st(A::AbstractMatrix, B::AbstractMatrix)
T = promote_type(eltype(A), eltype(B))
C = Matrix{T}(undef, (size(A, 1), size(B, 2)))
for x = axes(B, 2)

০ႣȘȱș

for y = axes(A, 1)
C[y, x] = @view(A[y, :])' * @view(B[:, x])
end

ȿɼɬɺɩɥɨɻȝˡႣȧȱȩ
end

● ப̢̲͒̍˚ৌˊ˒̞
͔̀͜ˡ̪͇͂͜ʏ̘
● Ș>˝ʽ໋ೠ˞ȝɀȝȝɁȝ
̅ࠥঁˊ˒ཅʿ͓͊͋
঎ဳႌ˵એ˝ʷˊ౽
஭ఆʷȘʵ˂˱˚ॸ
ࣗșș

C
end

# シングルスレッド版

my_matmul_st (generic function with 1 method)
julia> using BenchmarkTools
julia> @btime my_matmul_st(A, B) setup=(A=rand(100, 100); B=rand(100, 100));
374.714 μs (2 allocations: 78.17 KiB)
julia> @btime my_matmul(A, B) setup=(A=rand(100, 100); B=rand(100, 100));
329.578 μs (3239 allocations: 419.11 KiB)
julia> # 参考
@btime (A * B) setup=(A=rand(100, 100); B=rand(100, 100));
48.012 μs (2 allocations: 78.17 KiB)

০ႣȘȲș

julia> function threaded_map(fn, array::AbstractArray)
tasks = [Threads.@spawn(fn(v)) for v in array]
[fetch(task) for task in tasks]
end

ȿɻɸɥɿɶȝˡႣ

● ຢ੢ˡȝɵɥɸȝ‫܌‬ଲ̅
̢͇͔̪͕̬̲‫˙ˊص‬
˲˒
● ɪɩɼɧɬȧȩȝˢʆඁܰ̅
਌˖˙࠹‫ل‬ˡಃ̅਌
ඐˌ˾‫܌‬ଲ

threaded_map (generic function with 1 method)
julia> threaded_map(fib, 15:40) == map(fib, 15:40)
true
julia> using BenchmarkTools
julia> @btime map(fib, 15:40);
1.419 s (1 allocation: 272 bytes)
julia> @btime threaded_map(fib, 15:40);
652.908 ms (167 allocations: 13.77 KiB)

－81－

julia> mutable struct UnsafeCounter
count::Int
UnsafeCounter() = new(0)

ಣ֑ȘȰș
̢͕̬̲̤ʏ̾˞˗ʷ˙ȧȰȩ

● ˝̆˚˵̢͇͔̪͕̬
̲‫ܰ˙ˊص‬బೠ˽˞඀
ॖˌ˾̂˄˚ˢ˝ʷ
● ႣȦඁˋ༙ଲȘ̓̿
̡̘̱̐șˡ඲ဨ̅ࣃ
૵ˌ˾ૌࣱ˝˜
● _̢͕̬̲̤̋͜ʏ̾

end
julia> begin
counter1 = UnsafeCounter()
for n=1:1000
counter1.count += 1
end
counter1.count
end
1000
julia> begin
counter2 = UnsafeCounter()
Threads.@threads for n=1:1000
counter2.count += 1
end
counter2.count
end
650

# 想定より少ない値になってしまう（ことがある）

julia> mutable struct AtomicCounter
Threads.@atomic count::Int
AtomicCounter() = new(0)

ಣ֑Șȱș
̢͕̬̲̤ʏ̾˞˗ʷ˙ȧȱȩ

● ࣃ૵ˌ˾̾̌ʏ͔̲˞
ȿɥɼɷɵɮɧȝ̅຿˄˾
● ːˡ̾̌ʏ͔̲̅ࣃ૵
ˌ˾˛ˀ˵ȝ
ȿɥɼɷɵɮɧȝ̅຿˄˾
● _̢͕̬̲̤ʏ̾˞˝
˾

end
julia> begin
counter3 = AtomicCounter()
Threads.@threads for n=1:1000
Threads.@atomic counter3.count += 1
end
counter3.count
end
1000

# 期待通りの結果になる

julia> begin
counter4 = AtomicCounter()
for n=1:1000

# ↓シングルスレッドでも `Threads.@atomic` 必須

Threads.@atomic counter4.count += 1
end
counter4.count
end
1000

－82－

ಣ֑ȘȲș
ːˡఖ

● ˵˖˛ࢷʷ๛֋̢͕̬̲̤̅ʏ̾˞ˊ˒ʷ
_͖̬̘ܳࣇȘȿɴɷɧɱȝʘɴɷɧɱ̡̘̱̓̿̐ʙȝɦɩɫɮɶȝⳒȝɩɶɨș˞˻˾ȝ排他制御ȝ
̅ඃ෎ˌ˾ຖ္ʵ˽
● デッドロックȘᡦ˲ˡ૕మ˞˝˖˙˜ˡ̢͕̬̲˵͖̬̘̿૕మȘɁɴɷɧɱɩɨș˞
˝˖˙ˊ˱ʹș˞˵ಣ֑
● >ˆ˿˼ˡ઼ुˢ‫ں‬բ

০ਮȦ༊Ⴓ‫̺̒̾˾˻˞ص‬ʏ̢͇͜
‫ٶ‬஫

－83－

ʽూȦɑȬɔɽɩɩɶ࿨ూȧȰȩ
ʚɑȬɔɽɩɩɶ࿨ూʛ˛ˢ

●

●

ȷȬɔɽɩɩɶȝ࿨ూ˛ˢʆ
ȷウȷˡ̢͇࿠Ș̢̪̐
๣ș˞ȷࡰˡ̘̍ʏ̅͜ʆࢋ
ʷ˞ࣀ࠰ˊࣱ̂˝ʷ˻ʹ
˞ฉಌˌ˾̨̺ʏ͜ʿ‫ط‬
ೠ˽ʵ˾ʾʆ˛ʷʹ̺̣͔
࿨ూ
ˆ˿̅ȝɑウɑȝˡ̢͇࿠˞
ɑࡰฉಌˌ˾ʆ˛֨๘‫ˊص‬
˒ˡʿȝɑȬɔɽɩɩɶȝ࿨ూ
ゴȬɔɽɩɩɶȝ࿨ూˡ٪Ș֨Ⴃș

ʽూȦɑȬɔɽɩɩɶ࿨ూȧȱȩ

function nQueen(n::Int)
counter = AtomicCounter()

# 先ほど作ったカウンタを利用

for y=1:n
Threads.@atomic counter.count += nQueen_sub(n, 1, [y])

ɑȬɔɽɩɩɶ࿨ూ̅٪˂ȧȰȩ

●
●

Ɋɽɴɮɥȝˡ̜ʏ̲Ⴃ̅ৌˌ
ːˡȰȘப๊șȦ
͊̍͜ໞʽ˻˦मܴໞ

end
counter.count
end
function nQueen_sub(n::Int, k::Int, params::Vector{Int})
k == n && return 1
board = create_board(n, params)
counter = AtomicCounter()
for y=1:n
if issafe(board, y, k + 1)
Threads.@atomic counter.count += nQueen_sub(n, k+1, [params; y])
end
end
counter.count
end

－84－

ʽూȦɑȬɔɽɩɩɶ࿨ూȧȲȩ

function create_board(n, params)
board = falses(n, n)
for (x, y) in pairs(params)
board[y, x] = true

ɑȬɔɽɩɩɶ࿨ూ̅٪˂ȧȱȩ

●
●

Ɋɽɴɮɥȝˡ̜ʏ̲Ⴃ̅ৌˌ
ːˡȱȘ࢑๊șȦ
๣࿐̅ວˌȝɁɮɼɏɥɼɺɮʀȝ
̅୔୊ˌ˾‫܌‬ଲ˛ʆːˡ๣
࿐˞ঞ೾ֈಌ˞̜͇̅ಌ
ʷ˙վமʾ้̅೾ˌ˾‫܌‬
ଲ

end
board
end
function issafe(board::BitMatrix, y, x)
h, w = size(board)
any(board[y, :]) && return false
any(board[:, x]) && return false
any(board[y-i, x-i] for i=1:min(y, x)-1) && return false
any(board[y+i, x+i] for i=1:min(h-y, w-x)) && return false
any(board[y-i, x+i] for i=1:min(y-1, w-x)) && return false
any(board[y+i, x-i] for i=1:min(h-y, x-1)) && return false
true
end

julia> using BenchmarkTools

ʽూȦɑȬɔɽɩɩɶ࿨ూȧȳȩ

julia> @btime nQueen(8)
3.699 ms (75032 allocations: 2.95 MiB)
92

ɑȬɔɽɩɩɶ࿨ూ̅٪˂ȧȲȩ

●
●

Ɋɽɴɮɥȝˡ̜ʏ̲Ⴃ̅ৌˌ
ːˡȲȦ
০ࣞႣȘɑȼȷⳒȰȱ˱˚ș

●

>ঁ੒ˊ˒٪ˢ੬֞ȝバラエ
ティ解Ș٫രʈ๋ര˚ඁ֨঱
˚ˀ˾໲˵ࡰ ˞Ⴓ‫ޅ‬ș˚
ʵ˽ʆ基本解Ș٫രʈ๋ര˚
ඁ֨঱˚ˀ˾໲̅Ȱ˗˛ଲ
ʻ˾ཅཇș˚ˢ˝ʷˆ˛˞˵
ಣ֑

julia> @btime nQueen(9)
16.535 ms (326209 allocations: 13.11 MiB)
352
julia> @btime nQueen(10)
82.212 ms (1483298 allocations: 60.40 MiB)
724
julia> @btime nQueen(11)
436.700 ms (7371379 allocations: 303.26 MiB)
2680
julia> @btime nQueen(12)
2.611 s (39709620 allocations: 1.61 GiB)
14200

－85－

# ↓マルチスレッド化してみよう！
function nQueenMT(n::Int)

ʽూȦɑȬɔɽɩɩɶ࿨ూȧȴȩ

counter = AtomicCounter()

# 先ほど作ったカウンタを利用

for y=1:n
Threads.@atomic counter.count += nQueen_subMT(n, 1, [y])

ɑȬɔɽɩɩɶ٪ཇˡ༊Ⴓ‫ص‬ȧȰȩ

end
counter.count

●

●

ப๊ˡ̜ʏ̲
Șɶɔɽɩɩɶȧȩȝ‫܌‬ଲʆ
ɶɔɽɩɩɶɢɻɽɦȧȩȝ‫܌‬ଲș
̢͇͔̪͕̬̲̅‫ص‬Ș̢
͕̬̲༊Ⴓ‫ص‬șˊ˙˲˻
ʹ
̻̱͜Ȧऋ෌‫ڦ‬ʻ˒ɀɓɉ
Șȿɼɬɺɩɥɨɻʆȿɻɸɥɿɶȝ
˝˜ș̅঎ʻˣɒɌ

end
# ↓マルチスレッド化してみよう！
# ↓適宜編集して最適な挙動にしてみよう！
function nQueen_subMT(n::Int, k::Int, params::Vector{Int})
k == n && return 1
board = create_board(n, params)
counter = AtomicCounter()
for y=1:n
if issafe(board, y, k + 1)
Threads.@atomic counter.count += nQueen_sub(n, k+1, [params; y])
# Threads.@atomic counter.count += nQueen_subMT(n, k+1, [params; y])
end
end
counter.count
end

julia> using BenchmarkTools

ʽూȦɑȬɔɽɩɩɶ࿨ూȧȵȩ

julia> @btime nQueenMT(8)
?.??? ms (????? allocations: ?.?? MiB)
92

ɑȬɔɽɩɩɶ٪ཇˡ༊Ⴓ‫ص‬ȧȱȩ

●

●

julia> @btime nQueenMT(9)

০ࣞȘɑȼȷⳒȰȱ˱˚șˊ
˙˲˙ʆ̺̾̒ʏ̢͇͜ࢪ
ૅˊ˙ʷ˾ʾ‫˻ˊ෕ڤ‬
ʹ
Ș˓ʷ˒ʷȰȮȱⳒȰȮȲ˂˼
ʷˡ੪ၟৄ‫˾˝˞܊‬ˢˍș
>ʹ˱˂༊Ⴓ‫˚ص‬ˀ˙ʷ˝
ʷ˛‫̺̒̾˞ݧ‬ʏ̢͇͜೷
‫ش‬ȘȱⳒȰȯชˡ੪ၟৄ‫܊‬ș˞
˝˾ˆ˛˵ʵ˽˱ˌʆˇಣ
֑

?.??? ms (?????? allocations: ??.?? MiB)
352
julia> @btime nQueenMT(10)
??.??? ms (??????? allocations: ??.?? MiB)
724
julia> @btime nQueenMT(11)
???.??? ms (??????? allocations: ???.?? MiB)
2680
julia> @btime nQueenMT(12)
?.??? s (???????? allocations: ?.?? GiB)
14200

－86－

༬అȦ͉̩̋ʏ͔ˡཇ௿
● ̢͕̬̲ଲȝɑȝ˚༊Ⴓ‫˵˙ˊص‬ʆ੪ၟৄ‫܊‬ˢȝȰȮɑȝ˛ˢ˝˼˝
ʷ
● மదˡȝȸȯȣȝʿ༊Ⴓ‫˚ص‬ˀ˾˛ˊ˒ૌࣱʆ̺̾̒ʏ̢͇͜
ˢयి˚˵ȝȰȮȧȧȰȬȯȭȸȩȪȯȭȸȮɑȩȝชȝȘɑ@ ȝ˚ȝȰȯชș˞
ˊʾ˝˼˝ʷ
● ʚʷʾ˞༊Ⴓ‫ص‬຾‫ؼ‬෫˝ໞ໴̅એ˝˂ˌ˾ʾʛʿٕూ
● ॸࣗȦȝɚɮɱɮɸɩɨɮɥɞɯɥɠȝ͉̩̋ʏ͔ˡཇ௿ȝ൫

˱˛˴

－87－

࣠‫ݓ‬ˡ࿠ജ˛࿠ຢ

İ ̢͕̬̲༊ႳȘˡ‫ྂܤ‬șˡၟ٪
☑
☑
☑
☑

༊Ⴓ੪ၟ˛ˢ‫ط‬ʾ
̢͕̬̲ˡ‫ྂܤ‬ജ˝‫ڋ‬ෞ
̢͕̬̲ˡၚല˛୅࡮
Ș̢͕̬̲ˡ̞̘͔͒̍̾̍ș

İ Ɋɽɴɮɥ˞˻˾̢͕̬̲༊Ⴓˡၟ٪
☑ य೷࡮ˡ়ப੢຃
☑ ‫ڜ‬ਓɀɓɉ
☑ ̢͕̬̲̤ʏ̾Ș˛̰̬̲͖̬̘ș

ʽ๴˿ˈ˱˚ˊ˒

－88－

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

で始めるB-spline

堀川由人(@Hyrodium)
2023年7月11日

－89－

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

－90－

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

－91－

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

－92－

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

－93－

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

－94－

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

－95－

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

－96－

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

－97－

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

－98－

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

－99－

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

－100－

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

－101－

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

－102－

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

－103－

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

－104－

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

－105－

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

－106－

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

－107－

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

－108－

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

－109－

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

－110－

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

－111－

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

－112－

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

－113－

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

－114－

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

－115－

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

－116－

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

－117－

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

－118－

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

－119－

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

－120－

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

－121－

Ԡ༻਺ֶͷ֓೦Λ༻͍ͨ෺ཧϞσϧΛ͍͔ʹ Julia ͰϓϩάϥϛϯάΛߦ͏͔
チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

߱ដ େհ (େࡕେֶαΠόʔϝσΟΞηϯλʔ)
daisuke.furihata.cmc@osaka-u.ac.jp

਺ֶͱ෺ཧʹ͓͚Δ Julia ͷ‫۝@༻׆‬भେֶ IMI, 2023.07.11

1/19

Tips: ਺஋‫ࢉܭ‬ϓϩάϥϜΛ૊Ή্Ͱ Julia ͩͱ‫ͱ͍͜͠خ‬
ૈ͍ۙࣅղΛ͋Δఔ౓ྑ͍ۙࣅղ΁͚ۙͮΔ (‫܁‬Γฦ͠) ख๏ɾΞϧΰϦζϜ͸ͨ͘͞Μ͋Δ͠ɼ‫ݪ‬
ཧ΋͔ͬ͠Γͨ͠ɼ໌շͳ΋ͷ͕ଟ͍ɽ͔ͦ͠͠ͷ‫ޙ‬ɼ
͋Δఔ౓ྑ͍ۙࣅղ͔Βɼ૬౰ʹྑ͍ۙࣅղ ∼
= ΄΅ਖ਼͍͠ղΛಘΔ
⇐⇒ ඇઢ‫ํܗ‬ఔࣜͷࠜΛ‫ٻ‬ΊΔ
ํ๏͸݁ߏ࣮૷͕໘౗ɽ͜Εʹର͠ɼJulia ʹ͸ “NLsolve” package ͕͋Γɼͱͯ΋ॿ͔Δʜ
Պֶٕज़‫ࢉܭ‬ϓϩάϥϜͰ͸࿈ཱҰ࣍ํఔࣜͷղΛ‫ٻ‬ΊΔγʔϯ͕ଟ͍ͷʹɼ௨ৗͷίϯϐϡʔλ
‫ʹޠݴ‬͸ͦ͏͍ͬͨϥΠϒϥϦ͸ඪ४૷උͰ͸ͳ͍ɽϥΠϒϥϦΛ୳ͯ͠ɼؔ਺ΛϚχϡΞϧͰௐ
΂ͯɼ࢖͍ํΛ֬ೝͯ͠ʜ ໘౗ͩʂ
=⇒ Julia ͩͱ “LinearAlgebra” package ͕ඪ४૷උͳͷͰɼͱͯ΋ָɽ
ྫ: ߦྻ A ͱϕΫτϧ x, b ʹରͯ͠ Ax = b Λղ͍ͯ x Λ‫ٻ‬Ί͍ͨ৔߹ɼx = A \ b ͱॻͩ͘
͚Ͱྑ͍ʂ
‫ʹط‬ଟ͘ͷՊֶٕज़‫ܥ‬ͷϥΠϒϥϦ͕ଘࡏ͢Δʂ ͦͷ͏͑ɼdump ίϚϯυΛ͏·͘࢖͏ͱɼϥΠ
ϒϥϦͷ‫݁ࢉܭ‬ՌΛࣗ༝ʹநग़Ͱ͖ͨΓ͢Δɽ
=⇒ ϥΠϒϥϦͷιʔεΛಡ·ͳͯ͘΋ɼ಺෦݁ՌͳͲΛ஌Δ͜ͱ͕Ͱ͖Δɽ
2/19

－123－

ͯ͞ʜ ‫ڀݚ‬໨తɾഎ‫ܠ‬
Cahn–Hilliard ํఔࣜ‫ܥ‬౷ͷ໰୊ (૬෼཭‫ݱ‬৅ͳͲͷϞσϧํఔࣜ) ͷ਺஋ղੳΛ͍ͨ͠
movie(method of line ʹΑΔ਺஋ղ)
movie(DVDM ʹΑΔ਺஋ղ)

͜ΕΒͷ໰୊ʹ͸࣭ྔอଘɼΤωϧΪʔ‫ݮ‬গੑͳͲͷେҬతͳੑ࣭͕͋Δɽͦ΋ͦ΋
Cahn–Hilliard ʹ͔͗Βͣ͜͏ͨ͠໰୊͸ଟ͍ɽ
͜ΕΒΛอଘ͢Δ਺஋ղੳ๏ ߏ଄อଘ਺஋ղ๏ structure-preserving method ͷ‫ڀݚ‬͸ͦΕͳΓʹ
಺֎Ͱൃల͍ͯ͠Δ (ྫ: SciCADE ͱ͍͏ࠃࡍ‫ूڀݚ‬ձͰ͸͜ͷ಺༰ͷηογϣϯ͕͍ͭ΋͋Δ)ɽ

structure-preserving method ͸਺஋ղ͕༏Ε͍ͯΔ͜ͱ͕ଟ͍ɽ
ODE ͩͱϋϛϧτϯ‫ ܥ‬+α Ͱߟ͑Δ͜ͱ͕ଟ͍͕ɼPDE ͩͱม෼ߏ଄Λհͯ͠ߟ͑Δ͜ͱ͕ଟ
͍ (PDE ͷͦΕ͸཭ࢄม෼ಋؔ਺๏ͳͲͱশ͍ͯͯ͠೔ຊൃ)ɽ
͜ΕΒͷख๏͸‫ج‬ຊతʹඍੵΛ࢝Ίͱͨ͠ʮ࿈ଓ‫Ͱݶۃ‬ఆٛ͞ΕΔ operatorʯΛ཭ࢄతʹఆٛͯ͠
ͦͷҰ؏ੑΛ୲อ͢Δํ๏ɽ͋Δҙຯ਺ֶతʹ‫ີݫ‬ɽͨͩ͠‫ྔࢉܭ‬͸େ͖ΊͰɼ͕ͦ͜೉఺͕ͩɼ
ۙ೥Ͱ͸ߴ଎Խख๏΋ॆ࣮͖ͯͨ͠ɽ
্ͷ࿮૊Έ͸ͦΕͳΓʹॆ࣮͖͍ͯͯ͠ΔͷͰɼͦΖͦΖࠜຊతʹҧ͏੾Γ‫ޱ‬΋ߟ͍͑ͨ
3/19

Cahn–Hilliard ํఔࣜ ֓ཁ (1)
૬෼཭‫ݱ‬৅ͷϞσϧํఔࣜͷҰͭͰɼ఺ (x, t) ʹ͓͍ͯঢ়ଶ͕ A ૬,B ૬ͷ͍ͣΕʹ͍͔ۙΛࣔ͢෺ཧ
ྔ u(x, t) ʹର͠ҎԼͷΑ͏ͳ PDE Ͱ‫ه‬ड़͞ΕΔ΋ͷ (ଟ͘ͷ৔߹ɼA ૬ঢ়ଶΛ u = 0, B ૬ঢ়ଶΛ
u = 1 ͳͲͱ͢Δ)ɽ


d
∂u
for t ∈ [0, ∞), x ∈ Ω
(1)
= ∆ −ϵ∆u +
P (u)
∂t
du

ͨͩ͠, ∆ ͸ x ʹؔ͢ΔϥϓϥγΞϯɼ0 < ϵ ≪ 1, P (u) ͸ A ૬ঢ়ଶ, B ૬ঢ়ଶΛఈʹ΋ͭೋॏҪ‫ށ‬
‫਺ؔܕ‬, ‫ڥ‬ք৚݅͸આ໌ུɽ
͜ͷͱ͖ɼ(ద੾ͳ‫ڥ‬ք৚݅ͷ΋ͱͰ) ࣍ͷ 2 ͭͷେҬతੑ࣭͕͋Δɽ
࣭ྔอଘੑ

d
dt

ΤωϧΪʔ‫ݮ‬গੑ



u dx = 0


d
G(u, ∇u) dx ≤ 0
dt
ͨͩ͠ɼG ͸ΤωϧΪʔؔ਺ͱ‫ݺ‬͹ΕΔԼ‫ه‬ͷྔɽ
ϵ
G(u, ∇u) =
|∇u|2 + P (u)
2

－124－

(2)

(3)

(4)
4/19

Cahn–Hilliard ํఔࣜ ֓ཁ (2)
ͱ͘ʹΤωϧΪʔ‫ݮ‬গੑ͸ॏཁͰɼղͷଘࡏੑΛূ໌͢Δ਺ֶతͳΩʔͰ΋͋Δɽͦͷํఔࣜͱͷؔ࿈
ੑ͸໌ྎͰҎԼͷ௨Γɽ
ΤωϧΪʔؔ਺ G ʹରͯ͠ɼͦͷม෼ಋؔ਺͕
δG(u, ∇u)
d
= −ϵ∆u +
P (u)
(5)
δu
du
ͱಋग़Ͱ͖ΔͷͰɼ͜ΕΛ༻͍ͯ Cahn–Hilliard ํఔࣜ͸Լ‫ه‬ͷ‫͑׵͖ॻʹܗ‬ΒΕΔ (͜ͷ‫͕ࣜܗ‬ϙΠ
ϯτ)ɽ

∂u
∂t

=

∆

δG
δu

(6)

ΑͬͯΤωϧΪʔ‫ݮ‬গੑ͸


δG ∂u
d
G(u, ∇u) dx =
dx + ‫ڥ‬ք߲∗
dt
δu ∂t


 
 δG 2
δG δG
 dx + ‫ڥ‬ք߲∗ ≤ 0
=
∆
dx = − 
∇

δu δu
δu 

஫: ཭ࢄม෼ಋؔ਺๏͸͜ͷϖʔδΛ‫ؙ‬ʑ཭ࢄԽ͢Δ

(7)

∗ : ૝ఆ͞Ε͍ͯΔ‫ڥ‬ք৚݅ʹΑͬͯ 0 ͱͳΔɽ
5/19

Cahn–Hilliard ํఔࣜͷ௨ৗͷ਺஋ղ๏
PDE ʹର͢Δී௨ͷղ๏ɼಛघͳղ๏౳͕࢖͑Δɽ
method of line(ۭؒํ޲Λઌʹ཭ࢄԽͯ͠࿈ཱ ODE Ͱۙࣅ͠ɼODE ਺஋ղ๏Ͱղ͘).
͜ͷղ๏ͷղઆ pdf

࣌ؒํ޲ͷ཭ࢄԽ෯ ∆t Λখ͘͢͞Ε͹ҰԠ࢖͑Δ͕ɼ࣭ྔอଘੑɼΤωϧΪʔ‫ݮ‬গੑ͸ͱ͖͓Γ
ഁΕɼͦͷ͍ͤͰ෺ཧతʹΈ͓͔ͯ͠ͳ਺஋ղʹͳΓ͕ͪ (ͱ͘ʹ௕࣌ؒൃలͤͨ͞৔߹)ɽ
ߏ଄อଘ਺஋ղ๏ͱͯ͠ͷ཭ࢄม෼ಋؔ਺๏ (ࠩ෼๏ɼ༗‫ݶ‬ཁૉ๏౳͍Ζ͍Ζ͋Γɽઢ‫ܗ‬εΩʔϜ
΋͋Δ).
௕࣌ؒൃల΋ߟ͑Δͱ‫ࡏݱ‬ͷຊ໋͸ͪ͜Β͔ɽ‫ج‬ຊతʹ͖ͪΜͱಈ͘͠ɼ਺஋ղͷ༷ࢠ͕෺ཧతʹ
͓͔͍͠ͱ͍͏͜ͱ΋͋·Γແ͍ɽ਺஋ղͷҰҙଘࡏੑ΍҆ఆੑ͕ূ໌Ͱ͖Δέʔε΋͋Δɽͨͩ
͠ɼଟ͘ͷ৔߹਺஋εΩʔϜ͕࣌ؒํ޲ʹӄత (ະ஌ͷ਺஋ղʹରͯ͠࿈ཱඇઢ‫ํܗ‬ఔࣜΛղ͍ͯ
਺஋ղΛ‫ٻ‬ΊΔͱ͍͏‫ܗ‬ͷ͜ͱ) Ͱɼ‫ྔࢉܭ‬͸େ͖Ίɽ
ઢ‫ܗ‬εΩʔϜΛઃ‫͢ܭ‬Δํ๏࿦ (ଟஈԽ) ΋͋Δ͕ɼ໰୊ͷඇઢ‫ܗ‬ੑ͕ଟ߲ࣜͷ‫ܗ‬ঢ়Ͱͳ͍ͱ࢖͑
ͳ͍͜ͱɼ4 ࣍Λӽ͑Δͱ਺஋తͳෆ҆ఆੑ͕‫ͳ͘ڧ‬Δ͜ͱͳͲͷ੍‫͕͋ݶ‬Δɽ
ʮଥ౰ͳ਺஋ղΛߴ଎ʹಘΔʯͱ͍͏ํ๏࿦ʹ͍ͭͯ͸ɼॆ࣮͍ͯ͠Δͱ͸·ͩ‫͍ݴ‬೉͍ɽ
6/19

－125－

Cahn–Hilliard ํఔࣜͷ order reduction: ‫ݱ‬ঢ়
ଥ౰ͳ਺஋ղΛಘΔͨΊͷํ๏ͷҰͭʹɼorder reduction (໰୊ΛΑΓ௿࣍‫ݩ‬ͷ໰୊ʹஔ͖‫ͯ͑׵‬ղ͘)
͕͋Δɽ

Alikakos Β͕ ϵ → 0 ͷͱ͖͜ͷૈࢹԽաఔͷ૬ͷ‫ڥ‬ք໘ͷ‫ڍ‬ಈ͕ Hele–Shaw ໰୊ (ྖҬ‫ڥ‬ք໘
ͷ‫ڍ‬ಈ໰୊) ͷղʹऩଋ͢Δ͜ͱΛূ໌ † ɽ
order reduction ͱͯ͠͸্ͷ݁Ռ͸਺ֶతʹ͸΄΅ຬ଍ (໰୊ͷ࣍‫͕ݩ‬จࣈ௨Γ 1 ࣍‫ݩ‬Լ͕Δ)ɽ
ͨͩ͠ɼHele–Shaw ໰୊ͷ਺஋ղੳ͸໽հɽΉ͠ΖʮHele–Shaw ໰୊ͷ਺஋ղ๏ͱͯ͠
Cahn–Hilliard ํఔࣜΛղ͘ํ๏ʯ͕ఏҊ͞ΕΔঢ়ଶɽ
† N. D. Alikakos, P. W. Bates and X. Chen, Convergence of the Cahn–Hilliard equation to the Hele–Shaw model,
Arch. Rat. Mech. Anal. 128 (1994), 165–205.

Hele–Shawﬂow ͷ࣮‫ݧ‬໛ࣜਤͱ࣮‫ݧ‬ը૾ (K.Yoshii and Y.Sumino, arXiv 1904.10673, 2019. ΑΓҾ༻)

7/19

Cahn–Hilliard ํఔࣜͷૈ͍ order reduction Λ࡞Δʹ͸ (method 1)

എ‫ܠ‬ɾΞΠσΟΞͷελʔτ஍఺:
ඍ෼‫ز‬ԿΛϕʔεʹɼཻࢠ๏తͳ particle dynamics model ํఔࣜΛ࡞Δɽͭ·Γɼ֤஍఺Ͱͷ
u(x, t) ͕࣌ؒมԽ͢Δͱ͍͏໰୊Ͱ͸ͳ͘ɼ(ͨͱ͑͹) u ∼
= 1 ͳཻࢠ͕Ҡಈ͢Δɼͱ͍͏໰୊ʹ
ॻ͖‫͑׵‬Δɽ
͜ͷͱ͖ղܾ͢΂͖໰୊͸େ͖͘Θ͚ͯҎԼͷ̎ͭɽ
1
2

ʮղ u ͷ஋͕େ͖͍ ⇐⇒ ཻࢠ͕ͦ͜Ͱີू͍ͯ͠Δʯͱ͍͏ؔ܎ੑΛͲ͏΍ͬͯఆྔԽ͢Δ͔ɽ
ؔ਺ u(x, t) ʹର͢ΔϥϓϥγΞϯ ∆ ͳͲͷඍ෼࡞༻ૉΛɼཻࢠͷҐஔؔ܎ʹΑͬͯൃੜ͢ΔྗͳͲʹ
Ͳ͏΍ͬͯղऍ͠௚͔͢ɽ

࣮͸্ͷ̎ͭͷ໰୊ͱ΋ɼ
ʮཻࢠͷҐஔΛ฼఺, ੜ੒఺ͱۭͨؒ͠෼ׂʯΛ΋ͱʹͨ͠཭ࢄඍ෼Λ༻
͍Δ͜ͱͰۙࣅͱͯ͠ղܾՄೳ =⇒ Voronoi ෼ׂ (࣍εϥΠυ) ͳͲ͕୅දతɽ

8/19

－126－

ϘϩϊΠ෼ׂͱ͸: –‫ے‬ͷྑ͍ۭؒ෼ׂ–
೚ҙͷ఺͸΋ͬͱ΋͍ۙϘϩϊΠྖҬʹ‫·ؚ‬ΕΔͱ͢Δɽͭ·Γɺ฼఺ xi ʹΑΔϘϩϊΠྖҬ
Ωi ΛҎԼͷΑ͏ʹఆٛ͢Δɽ
def

Ωi = {x| ∥x − xi ∥ ≤ ∥x − xj ∥ for any j ̸= i}
j2
rij1
Ǎǡij1

j3

Ǣi

i

j1

lij1
ɾ฼఺, ੜ੒‫ۭؒ = ݩ‬Λ͜ͷ఺Λ΋ͱʹϘϩϊΠ෼
ׂ͢Δ
ɾϘϩϊΠ఺ = ϘϩϊΠྖҬͷ௖఺
ɾΩi ͱ Ωj ͕ྡ઀͍ͯ͠Δͱ͖


def 

=
Ωi ∩ Ωj  ,
rij
def
 l
= ∥x − x ∥.
ij

i

j

฼఺͔Β Voronoi ෼ׂ͕ਐΉ༷ࢠͷ movie

Voronoi ෼ׂ౳ͷ཭ࢄԽʹ൐͏ඍੵ෼ͷ཭ࢄԽ

9/19

ํఔࣜΛཻࢠ๏૬౰ʹਅ໘໨ʹ཭ࢄԽͨ͠ํఔࣜ
Ϟσϧํఔࣜͷશ༰͸ҎԼͷΑ͏ʹͳΔɽ




1
δG 
dxi
nij rij ,
= −
dt
4m
δu i
(Mad )ij is true


2q
δG 
= p(2ui − 1) + r(ui − 1)3 +

δu i
|Ωi |

(Mad )ij is true

ui

=



uj − u i
lij



rij ,

m
.
|Ωi |

ͨͩ͠ɼxi (t) ͸ i-൪ particle ͷҐஔ, m ͸ particle ͷେ͖͞ (໘ੵ), p, q, r ͸ Cahn–Hilliard ํఔࣜ
ͷఆ਺ɽ|Ωi | ͸฼఺ xi ʹΑΓੜ੒͞ΕͨϘϩϊΠྖҬ Ωi ͷେ͖͞Ͱɼߦྻ Mad ͸ϘϩϊΠྖҬͷ
ྡ઀ؔ܎Λྻ‫ྻߦ܎ؔ઀ྡͨ͠ڍ‬. nij , lij , rij ͸ྡ઀͍ͯ͠Δ৔߹ͷΈඇθϩͰɼͦΕͧΕ particle i
ͱ particle j ؒ‫ڥ‬քͷ (i ͔ΒΈͨ) ֎޲͖୯Ґ๏ઢϕΫτϧɼparticle ؒͷ‫཭ڑ‬ɼಉ‫ڥ‬քͷ௕͞Ͱ͋Δɽ
͜ͷϞσϧํఔࣜ͸਺ֶతʹ͸ࠜ‫͕͋ڌ‬Δ߹ཧతͳ΋ͷ͕ͩɼཻࢠͷҐஔ͕͖ۙͮ͗ͨ͢ΓɼྖҬͷ֎
΁ߦͬͯ͠·ͬͨΓ (ີ౓ −→ 0 Λ࣮‫͠ݱ‬Α͏ͱͯ͠) ͢Δͱ Voronoi ෼ׂ͕ෆ҆ఆʹͳΔͨΊɼ࣮ࡍ
ͷ‫ʹࢉܭ‬͸ͦ͏ͱ͏ͷણࡉ͕͞ඞཁͱͳΔɽ
͜ͷϞσϧͷ Julia ϓϩάϥϜ
͜ͷϞσϧʹΑΔ਺஋‫( ྫࢉܭ‬ಈը)

－127－

10/19

਺஋‫݁ࢉܭ‬ՌΑΓ‫ͯ͑͘ݟ‬Δ͜ͱ (method 1)

਺ֶతʹ͸͖ͪΜͱۙࣅ͍ͯ͠ΔͷͰɼ(‫͕ࢉܭ‬ՄೳͳΒ) ‫݁ࢉܭ‬Ռ͸ͦΕͳΓʹ৴པͰ͖ͦ͏ɽ
ີ౓ͷߴ͍ͱ͜Ζʹཻࢠ͕ͨ͘͞Μू·Γۭؒ෼ׂͷղ૾౓͕ߴ͘ͳΔͷͰɼϝογϡௐ੔ͷҰछ
ʹͳ͍ͬͯΔɽ
཭ࢄۙࣅͷࡍʹۭؒ࣍‫ʹݩ‬ґଘ͢ΔՕॴ͕ͳ͍ͷͰɼ̏࣍‫ݩ‬ɼ̐࣍‫Ͱݩ‬΋໰୊ͳ͍ɽ
఺ͷ‫ͳ͕ۙ͘཭ڑ‬Δͱ Voronoi ෼ׂͷ਺஋‫ࢉܭ‬͸ෆ҆ఆʹͳΓ͕ͪɽ
ۭؒͷʮ֎΁ʯग़͍ͯ͜͏ͱ͢Δཻࢠ͕ଘࡏ͠͏ΔͷͰɼ͜ͷରॲ͕ඞཁɽ

ͱ͍͏Θ͚Ͱɼ΋͏গ͠ʮ؆ܿͳϞσϧʯΛߟ͍͑ͨɽ

11/19

Cahn–Hilliard ํఔࣜͷૈ͍ order reduction Λ࡞Δʹ͸ (method 2)
എ‫ܠ‬ɾΞΠσΟΞͷελʔτ஍఺:

Cahn–Hilliard ํఔࣜΛ਺ֶతͳख๏Ͱ order reduction ͢Δͷ͸͍ͬͨΜࢭΊͯʜ
΋ͱ΋ͱͷ૬෼཭‫ݱ‬৅ʹ͓͚ΔʮૈࢹԽաఔ (coarsning process)ʯͷΈΛ࠶‫͢ݱ‬Δૉ๿ͳϞσϧΛ
ߟ͑Δ
ͦͷࡍɼগͳ͘ͱ΋࣭ྔอଘੑΛ࠶‫͢ݱ‬Δ (ͦ͏͠ͳ͍ͱ෺ཧతʹ͔ͳΓ‫͑ݟʹົح‬Δ)
ૈࢹԽաఔΛ‫͢࡯؍‬Δͱ:
શମʹʮྲྀମʯ‫ݱ‬৅ͷΑ͏ͳ‫ڍ‬ಈ
গྖҬಉ࢜ʹओʹ‫ʹ཭ڑ‬ґଘͨ͠Ҿྗతͳ૬‫͕͋༻࡞ޓ‬Γͦ͏
ྲྀମత‫ڍ‬ಈͱߟ͑Δͱɼۙ͘ͷྖҬΛʮ௚઀ӽ͑ͯʯྗ͕ಇ͘Ϟσϧʹ͸ҧ࿨‫ײ‬
ྲྀମͰ͸ɼ෺࣭ಉ࢜ͷ૬‫ʹ༻࡞ޓ‬ΑΓྗ͕ো֐
෺Λӌճͯ͠఻ΘΔ ⇒ ӌճͯ͠ͷྡ઀ؔ܎͕
੒ΓཱͭΑ͏ʹ͍ͨ͠ʜ
ͱ͍͏Θ͚Ͱɼ
12/19

－128－

particledynamics model of Cahn–Hilliard ํఔࣜ (1/2)
u∼
= 1 Ͱ൒‫͕ܘ‬Ұఆͷ‫( ٿ‬particle) ͷҐஔҠಈͷΈΛ‫ه‬ड़͢ΔϞσϧΛߟ͑Δɽ
particle ͸૿‫͍͠ͳ͠ݮ‬ɼେ͖͞΋มΘΒͳ͍͠ɼॏͳΒͳ͍ɽ͜ΕʹΑ࣭ͬͯྔอଘੑΛ֬อɽ
particle ಉ࢜͸ɼ‫͍ޓ‬ͷ‫ʹ཭ڑ‬Ԡͯ͡Ҿ͖߹͏ (ͱΓ͋͑ͣ‫ٯ‬ೋ৐ଇ)
ͨͩ͠ɼparticle ͕ॏͳΒͳ͍Α͏ʹඇৗʹ͍ۙ৔߹͸ମੵഉআޮՌతͳྗ͕ಇ͘.
ͨͩ͠ɼۙ͘ͷ particle ΛԿॏʹ΋ӽ͑ͯͷྗ͸͋·Γಇ͔ͳ͍ (؇࿨ྡ઀ؔ܎)
͜Ε͸ྫ͑͹Լ‫ه‬ͷΑ͏ʹͯ͠ϞσϧԽ͢Δɽ
1

2

্ۭؒͷ particle ΛͦΕͧΕ฼఺ͱۭͯؒ͠Λ Voronoi ෼ׂ͠ɼVoronoi ྖҬͷྡ઀ؔ܎͔Β࡞ΒΕΔ
ྡ઀ߦྻΛ Mad ͱ͢Δɽ
def

Ϟσϧύϥϝʔλ p (ਖ਼੔਺) ʹର͠ɼ؇࿨ྡ઀ߦྻΛ Mp = (Mad )p ͱ͠ɼMp ʹΑͬͯྡ઀ͯ͠
͍Δͱ൑அ͞ΕΔ particle ؒͰͷΈྗ͕ൃੜ͢Δͱ͢Δ.

Խֶ෺࣭ͷқಈ౓ͷ֓೦ͷΑ͏ʹɼparticle ʹಇ͘ྗ ∝ particle ͷ଎౓
ྲྀମ఍߅ͷΑ͏ʹɼparticle ͷ଎౓ʹ͸্‫͋ݶ‬Γ
ʜͱ͍͏ɼ͍Θ͹ re-modeling Λߦ͏
13/19

particledynamics model of Cahn–Hilliard ํఔࣜ (2/2)
Ϟσϧํఔࣜͷશ༰͸ɼ۩ମతʹ͸ҎԼͷ௨Γɽ


tanh(C2 ∥f i (t)∥)
dxi (t)
= C1
f i (t),
dt
C2 ∥f i (t)∥

f i (t) =
f ij (t),
(Mp )ij is true

f ij (t)

=





xj − x i
∥xj − xi ∥3
xj − x i

 −2C3
∥xj − xi ∥3
C3

:

∥xj − xi ∥ ≥ rcri ,

:

∥xj − xi ∥ < rcri .

ͨͩ͠ɼxi (t) ͸ i-൪ particle ͷҐஔ, C1 , C2 , C3 ͸ਖ਼ఆ਺ɽ
def

ߦྻ Mp ͸ Mp = (Mad )p (p ͸ਖ਼੔਺, 3 ఔ౓)ɼͦͷ‫ࢉܭ‬աఔͰඞཁͳਅِ஋‫( ࢉܭ‬ϒʔϧ஋ͷ‫ܭ‬
ࢉ) ͸ a · b = a and b, a + b = a or b ͱͳΔ.
ߦྻ Mad ͸ϘϩϊΠྖҬͷྡ઀ؔ܎Λྻ‫Ͱྻߦ܎ؔ઀ྡͨ͠ڍ‬ɼ฼఺ xi ʹΑΓੜ੒͞ΕͨϘϩϊΠ
ྖҬ Vi ͱ฼఺ xj ʹΑΓੜ੒͞ΕͨϘϩϊΠྖҬ Vj ͱ͕ྡ઀͍ͯ͠Δͱ͖͸ (Mad )ij ͕ਅɼͦ͏Ͱ
ͳ͚Ε͹ِͰ͋Δɽ
14/19

－129－

ϘϩϊΠྖҬͷྡ઀ؔ܎

ਤͷத৺ͷ੺͍ϘϩϊΠྖҬ Vi Λʹରͯ͠ͷྡ઀ؔ܎ߦྻ (Mad )p (p ͸ਖ਼੔਺) ͷҙຯ͸ʜ
̍ ͕ॻ͔ΕͨྖҬ ʜ (Mad )1 ͷ (i, j) ੒෼͕ʮਅʯͱ
ͳΔϘϩϊΠྖҬ Vj .

2 ͕ॻ͔ΕͨྖҬ ʜ (Mad )2 ͷ (i, j) ੒෼͕ʮਅʯͱ
ͳΔϘϩϊΠྖҬ Vj .
3 ͕ॻ͔ΕͨྖҬ ʜ (Mad )3 ͷ (i, j) ੒෼͕ʮਅʯͱ
ͳΔϘϩϊΠྖҬ Vj .
ʜ ஫໨ʂ ͜ͷྡ઀ؔ܎͸ʮӌճ࿏ʯΛ࣋ͪ͏Δ
ˠ ྲྀମͷʮճΓࠐΈʯ‫ڍ‬ಈΛϞσϦϯάͰ͖Δɽ
ۭؒͷϘϩϊΠ෼ׂͷྫɽ
͜ͷղ๏ͷ julia code

>> movie (ॳ‫ظ‬഑ஔ͕ random)

15/19

࣮‫݁ݧ‬Ռ: ૈࢹԽաఔΛͦΕͳΓʹ re-modeling Ͱ͖ͨͷ͔ʁ

ࠨ: particle ϞσϧʹΑΔ਺஋ղ, ӈ: Cahn–Hilliard ํఔࣜʹΑΔ਺஋ղ
>> movie (ॳ‫ظ‬഑ஔ͕ pi ‫ܗ‬ঢ়)

16/19

－130－

਺஋‫݁ࢉܭ‬ՌΑΓ‫ͯ͑͘ݟ‬Δ͜ͱ (method 2) 1/2

ඇৗʹ஗͍‫ڍ‬ಈ͕ൃ‫͍ͯ͠ݱ‬Δɽ͜ͷ‫ڍ‬ಈ͸ྖҬͷҠಈɾม‫ܗ‬աఔͰ‫ݟ‬ΒΕΔಛ௃ͷҰͭͰ͋Δɽ
ྖҬͷҠಈ଎౓͕ͦͷαΠζʹґଘ͢Δ (େ͖͍΄Ͳ஗͍). ͜͏ͨ͠ޮՌΛϞσϧʹ௚઀औΓೖΕ
͍ͯΔΘ͚Ͱ͸ͳ͍ͷͰ‫ڵ‬ຯਂ͍ɽ
࣭ྔอଘଇ͕੒Γཱͭ. ͜Ε͸Ϟσϧͷઃ‫͔ܭ‬Βͯ͠౰વɽ
ܾఆํ๏͕ෆ໌ͳະ஌ύϥϝʔλ͕‫·ؚ‬ΕΔϞσϧͰ͋Δ. ۩ମతʹ͸ɼC1 , C2 , C3 , p ͳͲ͕ະ
஌ύϥϝʔλɽ͜ΕΒΛͲ͏ܾΊΔͱ Cahn–Hilliard ํఔࣜͷղ‫ڍ‬ಈΛۙࣅͰ͖Δͷ͔ɼ͸
future work.
ྖҬͷੜ੒ɾফ໓ϝΧχζϜ͕ཅతʹ͸͜ͷϞσϧʹ͸૊Έࠐ·Ε͍ͯͳ͍ɽͨͩ͠ɼٖࣅతʹ͜
͏ͨ͠ϝΧχζϜ͕࠶‫ݱ‬͸͞Ε͍ͯΔ༷ࢠɽ
ͳ͓ɼ͜ͷܽ೗͸ࠓճͷϞσϧͱͯ͠͸ର৅‫ݱ‬৅͔Β੾Γࣺ͔ͯͨΒ౰વͰ͋Δ͕ɼ૊ΈࠐΊΔͳ
Β͹ͦΕʹӽͨ͜͠ͱ͸ͳ͍ɽ(ˠ future work)

17/19

਺஋‫݁ࢉܭ‬ՌΑΓ‫ͯ͑͘ݟ‬Δ͜ͱ (method 2) 2/2
Cahn–Hilliard ํఔࣜͷղͷ‫ڍ‬ಈΛγϛϡϨʔτ͢ΔͨΊͷγϯϓϧͳૈࢹԽϞσϧͱͯ͠
particle motion model Λઃ‫ͨ͠ܭ‬ɽ
਺஋‫ྫࢉܭ‬ΛΈΔͱ͜ͷϞσϧͰྖҬͷҠಈɾม‫ʹܗ‬஗͍‫ڍ‬ಈ͕ൃ‫͓ͯ͠ݱ‬ΓɼCahn–Hilliard ํ
ఔࣜͷૈࢹԽϞσϧͱͯ͠‫ظ‬଴Ͱ͖Δɽ
Ϟσϧʹ‫·ؚ‬ΕΔະ஌ύϥϝʔλͷܾఆํ๏͸‫ݱ‬ঢ়ෆ໌ɽ

Cahn–Hilliard ํఔࣜͷେҬΤωϧΪʔʹ૬౰͢Δ͜ͷϞσϧͰͷಛ௃ྔΛఆٛɾ‫͠ࢉܭ‬ɼͦͷ‫ڍ‬
ಈΛௐ΂Δඞཁ͕͋ΔɽͳͥͳΒ Cahn–Hilliard ํఔࣜͷେҬΤωϧΪʔ͸࣌ؒൃలʹ൐͍‫ݮ‬গ
͢ΔࢄҳྔͰ͋Δ͜ͱ͕൑໌͍ͯ͠ΔͨΊɽ
3 ࣍‫ྖݩ‬ҬͰͷ͜ͷϞσϧͰͷ਺஋‫ࢉܭ‬ͷ࣮‫ݧ‬౳Λߦ͍͍ͨɽ
ͳ͓ɼԼ‫ه‬ͷΑ͏ʹɼۭؒ࣍‫ݩ‬Λ্͛ͯ΋ɼ‫ࢉܭ‬ίετ͸ࢦ਺తʹ૿େͨ͠Γ͠ͳ͍ɽ
2 ࣍‫ݩ‬໰୊Ͱ͸ۭؒͷϘϩϊΠ෼ׂͷ‫ࢉܭ‬ίετ͸ଟ͘ͷΞϧΰϦζϜͰ O(n log n).
ͨͩ͠ n = # points.
Delaunay ࡾ֯෼ׂΛߦͬͯϘϩϊΠ෼ׂΛ‫͢ࢉܭ‬Δ Bowyer–Watson ΞϧΰϦζϜ͸ɼ໰୊ͷ࣍‫ʹݩ‬Α
Βͣ (͓͓ʂ ), ී௨͸‫ࢉܭ‬ίετ͕ O(n log n).
ͨͩ͠࠷ѱͷέʔε (఺഑ஔ͕ୀԽ͍ͯ͠Δ৔߹) Ͱ͸ O(n2 ) ʹͳΓ͏Δ.
18/19

－131－

·ͱΊ

ෳࡶͳ໰୊ʹରͯ͠΋ɼJulia + jupyter ‫ͱͩڥ؀‬ɼগͣͭ͠ࢼߦࡨ‫͢ޡ‬Δ͜ͱͰҙ֎ʹ‫ʹָؾ‬ϓ
ϩάϥϜͷ‫׬‬੒ʹ͚͗ͭ͜ΒΕΔɽ
૝૾͢ΔΑΓଟ͘ͷϥΠϒϥϦ͕༗Γɼ͔ͳΓʮཔΔʯ͜ͱ͕Ͱ͖Δɽ
ଞͷ‫ʹޠݴ‬͸΋͏໭Εͳ͍ʂ

19/19

－132－

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

－133－

量子古典融合アルゴリズムが拓く計算物質科学 (2023年4月〜)
領域代表: 品岡寛

https://qc-hybrid.github.io/

概要

計算物理は、桁違いに異なる長さスケールの共存を扱う

例 量子物性物理 (低温物理現象 vs 101 eV ~ 105 Kのバンド幅)
宇宙ひも (プランクスケールの細さ vs 宇宙スケールの長さ)
乱流

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

Quantics tensor train

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

スパースモデリング

グリーン関数の

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

SparseIR.jl

MSSTA.jl (仮名)
2009年に応用数学で提案され、ここ数年で流

品岡等によって2017年に提案され、虚時間形

体力学、量子力学等への適用が始まる

式の場の量子論・第一原理計算に応用が進む

－134－

発表の構成

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

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

量子多体系
Bare Coulomb interaction

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

グランドチャレンジ

密度汎関数 (Kohm-Sham), 摂動論など

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

－135－

情報圧縮技術がなぜ必要？ (2/2ページ)
2つの粒子が散乱される過程

桁違いに異なるスケールの混在
→単純メッシュではデータ量が発散

(r′′, t′′)
(r, t)
t

有効相互
作用

→圧縮＆演算する技術が必要

(r′′′, t′′′)
(r′, t′)

流体力学計算等でも同じような問題

テンソルネットワーク、スパースモデリング･･･
→試行錯誤が迅速で、並列計算までカバーする
Juliaが最適

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

－136－

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=4
Bond dimension

D=8

N. Gourianov et al., Nat. Comput. Sci. 2, 30 (2022)

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

José I. Latorre, arXiv:quant-ph/0510031v1

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

2022年8月にアイデア思いつく、2022年10月に投稿 (Juliaだから可能だった)

テンソルネットワーク (1/3)
スカラー

ベクトル

a

行列

3階テンソル

aijk

aij

ai

j

i

i

i

j

k

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

i

テンソルネットワーク

A

j

x

=

∑

Aij xj

j

A
i

α

B
j

β

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

C
k

=

∑
α,β

AiαBjα,βCkβ

上・下付きにも意味はない。

－137－

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

A

B

C

i

j

k

→
i

j

k

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

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

A = UΣV †
A

=

大きさは特異値の数

U

V†

=

ΣV †

U

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

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

x1

xN

≈

x1

xN

⋯

x1

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

－138－

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 k

U

V

i k

j

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 < 1
Binary coding: x = (0.x1x2⋯xR)2 = x1 /2 + x2 /22 + ⋯ + xR /2R

f(x)

(0.0⋯0)2 (0.0⋯1)2

(0.1⋯1)2

f(x1, ⋯, xR) ≈

D1

∑

⋯

xR

⋯

α1=1

≈
x1

x

1

0

x1

DR−1

∑

αR−1=1

α1

－139－

Fx(1),1α Fx(2)
⋯Fx(R),α
1,α1α2
1

⋯

⋯

1

R

αR−1

R−11

xR

QTTのまとめ (2/3ページ)
Same length scale

f(x, x′)

≈
x1

⋯

xR

x′1

x′R

⋯

x1

x = (0.x1x2⋯xn⋯)2 ∈ [0,1)
x′ = (0.x′1x′2⋯x′n⋯)2 ∈ [0,1)

x′1

⋯

xR−1 x′R−1 xR

x′R

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

QTTのまとめ (3/3ページ)
指数関数
2

n

f(x) = e −x = e −x1/2e −x2/2 ⋯e −xn/2 ⋯
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.

多項式
単位行列

M. Ali and A. Nouy, “A. Approximation Theory of Tree Tensor Networks:
Tensorized Univariate Functions”, Constr Approx (2023)

D≤1+p

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

0

=

⨂

⨂

D=1
⨂

⨂

⨂

⋯

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

－140－

Written by Miles Stoudenmire

場の量子論データの圧縮

HS et al., Physical Review X 13, 021015 (2023)
2P vertex functions

1P momentum space

Spectral function
(a) 0.8

Ω(!)

0.6

AF insulator
Doped Mott insulator

0.4
0.2
0.0

°5

0
!

5

(b)

Ω(!)

0.4

Impurity model

0.0
°0.025 0.000 0.025

0.2

0.0

0.2

°5

0
!

5

Bond dimension

(c)

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

Multipolar susceptibility

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

AF insulator
Doped Mott insulator
Impurity model

103

Nonequilibrium system (real-time Green’s function)

102

CeB6

101
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
Bond b

圧縮したまま(場の量子論)計算ができるか？

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
t1
tR
⋯

A

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

B

－141－

δ

t′1

t′R =

t′1′

t′R′

δ

δ

Juliaライブラリの実装 MSSTA.jl (multiscale space-time ansatz)
まだ、実験的なライブラリです (○ユニットテスト、△ドキュメンテーション)
https://github.com/shinaoka/MSSTA.jl

• フーリエ変換などの演算を簡単に実装可能なライブラリ
• バックエンドとして、ITensors.jlを利用
Juliaをなぜ選んだのか？
• forループの使用が避けられない
• 戦略的にJuliaで場の量子論計算のエコシステムを作る
• (ITensors.jlの開発者と交流があった)

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

－142－

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表現の宣伝予定！

－143－

スパースモデリング 疑似多倍精度による積分方程式の解法
「虚時間」に特化した最適な基底の構成

疑似多倍精度演算: MultiFloats.jl

疑似四倍精度

• 倍精度小数点2つの和として、浮動小数点を表現 (例: 1 + 10-100は、有効桁数が16桁程度の倍精度では表現不能)
• 真の四倍精度のソフトウェアエミュレーションよりも速い (でも、倍精度演算の2桁ぐらい遅い)

Float64は桁落ち

桁落ちしていない

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

－144－

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τ

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

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

－145－

虚時間形式の相関関数における情報損失
G(τ) = − ⟨Tτ A(τ)B(0)⟩ =
虚時間の量

∞

∫−∞

dωK(τ, ω)ρ(ω)
実周波数(時間)の量

系に寄らない積分核
−τω

K(τ, ω) =

e
1 + e −βω

フェルミオンの場合

高周波数の情報をτに伝搬しにくい

β = 1/T

虚時間の相関関数は、実時間の相関関数に比べて、高周波数の情報が落ちている。
→ 積分核を使って、虚時間の最適な正規直交基底を作れるか？

(a) 1

ωmax = 1

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

l=0

SlF

S U (τ)Vl(ω)
∑ l l

for

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

Orthonormal intermediate-representation bases

G(τ) =

∑

GlUl(τ) ρ(ω) =

l=0

ρ V (ω)
∑ l l
l=0

Gl = − Sl ρl

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

0.2

10−8
10−12

Number of N
coefficients

K(τ, ω) =

10−4

∞

(b) 0.4

β = 100
β = 1000
β = 10000

UlF(τ )

積分核の特異分解による基底の構成

104

−0.2

0

25

50
l

75

100

0.0

0.2

Matsubara frequency with tailLegendre
O(β)
Chebyshev
Legendre and Chebyshev
IR
O( β)

i!n

103

O(log β)

102
101

0.0

Intermediate representation (IR)
102

104

106

Ø

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

͑ΔͨΊʹ͸ɺ‫ؒ࣌ڏ‬ͷ৔߹ͱಉ༷ʹɺDyson ํఔࣜΛίϯύΫτͳ‫ج‬ఈͷੑ࣭Λར
ͯղ͘ํ๏͕ඞཁͱͳΔɻ

Ҏ্ͷٞ࿦Λཁ໿͢Δͱɺల։܎਺ͷܾఆ͕ՄೳͱͳΔ͚ͩͷ τ ɺiωn ͕ಘΒΕΕ

ଞͷ
τ ɺiωn ʹؔ͢Δ‫ࢉܭ‬͸ෆཁʹͳΔͱ͍͏͜ͱͰ͋Δɻզʑ͸ɺIR ‫ج‬ఈؔ਺ͷಛ
－146－

ͳߏ଄Λ࣋ͭΑ͏ͳҐஔͷ τ ΛબͿ͜ͱͰɺల։܎਺ͷܾఆʹඞཁͳ τ ͷ਺Λɺ‫ج‬

応用の拡がり

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

Analytic continuation

Migdal-Eliashberg theory

• 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

Sparse modeling (SpM)

• 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

Extension to the two-particle level

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)

SparseIR.jl

• 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ライブラリの使用も含み
ます

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)

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

• 基底関数の区分多項式 (Legendre)による展開: 特異分解→特異値分解
• 疑似四倍精度演算による特異値分解 MultiFloats.jl + GenericLinearAlgebra.jl
• 基底の構成が終われば、後は倍精度演算で十分

－147－

Python実装 (sparse-ir)との比較
Pure Python 実装: https://github.com/SpM-lab/sparse-ir
numpyの疑似四倍精度浮動小数点拡張 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のデモンストレーション

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

－148－

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

基底関数の精度は大体10桁

チュートリアル集
• 典型的な計算のサンプルコ
ード

• Jupyter Book
• 手元でのビルド環境は、
VS Code + Dockerで用意
https://spm-lab.github.io/sparseir-tutorial/

－149－

解説 + コード例
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万円)
• 別途、家賃補助、通勤手当あり

－150－

数論における Julia の援用

IMI 共同利用「数学と物理における Julia の活用」

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

横山 俊一 (東京都立大学)
s-yokoyama@tmu.ac.jp
2023年7月11日 ※ 本日38歳になりました
内容の一部は JSPS 科研費・基盤(C) 20K03537

Julia 言語を用いた新しい計算機数論システムの開発とその応用
(代表：横山 俊一) の援助を受けています

代数系（数論）における Julia の立ち位置
• 幾何系や解析系に比べると 活用状況・認知度ともに低い 状況
• 既存の CAS 開発においては、1つの分野（代数的対象物）に特化 したものが

多く、長期間使い込んだユーザが異なるプラットフォームに移行するのは困難
例： Magma（C, 次頁）、Sage（Python, 次々頁）、Pari/GP（C）etc.

• そもそも Julia 言語で書かれたパッケージはまだ少数
• それでも開発速度は徐々に上がってきている：

例： AbstractAlgebra.jl, Gap.jl, GroebnerBasis.jl etc.

• 正しく実装すれば、多項式代数においても高速

－151－

Magma
• 豪 シドニー大学 を中心として開発運営が進められている
計算代数システム。前身は Cayley (1982) で、
正式ローンチは 1993 年
有償だがプロプライエタリライセンス（非営利）

• 主として 数論および計算機代数に特化 した関数群
圏論 に基づいた設計思想

• 2013年から Intel AVX 演算に対応（Linux 64bit OS のみ）
• Bosma-Cannon-Playoust, The Magma algebra system I: the user
language (CANT/London, 1997)

Sage/CoCalc/SageMathCell
• 米 ワシントン大学 を中心として開発運営が
進められて「いた」汎用数式処理システム

• 数論幾何学の研究者 William Stein が開発責任者

1990年代には Magma の HECKE パッケージ（モジュラー形式の計算）の
開発に従事、2005年の正式ローンチの際に Python へ移植

• web ベースの UI 開発の変遷：

Sage Notebook -> salv.us -> SageMathCloud -> CoCalc
SageMathCell は CoCalc の簡易版

• 2019年に大学を退職し、SageMath Inc. CEO/ファウンダーとして運営

－152－

Nemo
• 独 カイザースラウテルン工科大学 を中心として開発運営が
進められている、Julia で開発された数論計算パッケージ
および C/C++ wrapper 集

• Singular の開発本拠地でもある（→ Singular.jl pkg ）
• 中心メンバ： Claus Fieker, William Hart, Tommy Hofmann,
Frederik Johansson, Marek Kaluba, Carlo Sircana

• FHHJ, Nemo/Hecke: Computer Algebra and
Number Theory Packages for the Julia
Programming Language, ISSACʼ17.

Nemo
• 2018年、Nemo は2つのパッケージに分離：

Nemo.jl Flint, Arb, Antic 限定 wrapper
AbstractAlgebra.jl 上記外の generic なパッケージ

• Julia 単独でも

- Z, Q, R, C（Flint/Arb）

- 行列計算, FFT, 並列計算 etc.
などが扱えるが、これに加えて整数論に特化した機能をもつ：
- 代数体（Antic）, 有限体（Flint）
- 1変数／多変数多項式環（generic）
- p進体・関数体（Flint, ﬁnite precision）

－153－

Nemo
• Flint を用いることなく generic に生成されたもの
(例：係数環 R) は AbstractAlgebra.jl に格納

• Flint や Arb 単独で高速処理できるものは、すべて
Julia 言語で実行するだけ

• より高速な実装を目指すならば、まず AbstractAlgebra.jl のリストを確認

• 例： S = Q̲p（p進数体, precision m）

m=p^r 以外の prec が代入されるとエラーを返す

－154－

Abstract-type hierarchy

• 例えば xxx.Field は xxx.Ring に属する

- Magma では異なる category として処理
- Sage では構造が 逆： xxx.ring は xxx.Field の部分集合

Hecke
• Nemo（と AbstractAlgebra パッケージ）に準じて開発されている
整数論 に特化した計算パッケージ

• 1990 年代に William Stein（Sage/CoCalc 開発者）が Magma で

開発していた、モジュラー形式に関するパッケージ HECKE とは別物

• https://github.com/thofma/Hecke.jl（ソースコード）
• https://thofma.github.io/Hecke.jl/latest/（ドキュメント）

－155－

－156－

Magma の Pollardʼs Rho + ECM

－157－

Nemo の Pollardʼs Rho + ECM

注： Magma ではさらに MPQS が使える

－158－

－159－

－160－

Magmaʼs NF-CG-RCF-NF

ベンチマーク例 (from NemoCas page)

－161－

ベンチマーク例 (from NemoCas page)

実はほぼ同一コード

ベンチマーク例 (from NemoCas page)

－162－

ベンチマーク例 (from NemoCas page)

ベンチマーク例 (from NemoCas page)

－163－

NemoCas Project
• http://nemocas.org/（Nemo/Hecke 開発ページ）
• Magma, Sage, Pari/GP と比較したいくつかのベンチマーク結果が閲覧可能
• 最新版は開発者用 Git リポジトリから

Nemo: https://github.com/Nemocas/Nemo.jl
AbstractAlgebra: https://github.com/Nemocas/AbstractAlgebra.jl
Hecke: https://github.com/thofma/Hecke.jl
Singular: https://github.com/oscar-system/Singular.jl/

• 試しに使ってみたいが導入に失敗した方へ：

Nemo は Jupyter Notebook 経由で Sage 上で利用可能

OSCAR

－164－

OSCAR
• 主として NemoCas プロジェクトから派生した、

Julia を用いた新計算代数システム開発プロジェクト

• GAP（Gap.jl）, Singular（Singular.jl）, Polymake（Polymake.jl）,
Antic（Nemo.jl / Hecke.jl）の4ライブラリを主軸に構成

• プロジェクトリーダ：

- Wolfram Decker（TU Kaiserslautern：NemoCas メンバ）
- Claus Fieker（TU Kaiserslautern：NemoCas メンバ）
- Max Horn（TU Kaiserslautern）
- Michael Joswig（TU Berlin）

• DFG, TUK に研究所 TRR-195 を2017年に開設

OSCAR

－165－

OSCAR
OSCAR
OSCAR.jl

Julia

number theory

group and
representation
theory

ANTIC

tropical and
polyhedral
geometry

algebraic
geometry and
commutative
algebra

noncommutative
algebra and
free probability
theory

xGAP.jlx

Polymake.jl

Singular.jl

GAP

polymake

Singular

Cornerstones

CAP
CHEVIE

GB.jl

mptopcom

PLURAL
Letterplace

GMP

FLINT

cdd

bliss

GPI−Space

GBLA

⋯

OSCAR

• 計算代数／数論における連携処理をシームレスに行う
Oscar.jl を提供

• 数論統合システム LMFDB 上のデモプログラムに OSCAR コードが新規採用

－166－

LMFDB

LMFDB

－167－

OSCAR
• すべての情報は OSCAR の web ページから参照：
https://www.oscar-system.org/

• インストールが少し難しいので注意：

- Julia に加えて Xcode ツール が必要（Mac OS の場合）
- Julia を立ち上げてから Pkg ツール を用いる
- インストール・ダウンロード (アーティファクト)・プリコンパイルに
少々時間がかかるので気長に待つ

• 終わったら using Oscar で、あとは楽しむ

ただしパス設定がシビアなので、慣れるまで頑張る
※ Jupyter Notebook 環境を使うと少しは導入がラク

－168－

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

Juliaによる科学技術計算:
大規模並列計算について

国立研究開発法人日本原子力研究開発機構 システム計算科学センター

永井佑紀

1

2

アウトライン
自分野(物性物理学分野)におけるコード開発の状況：みんな大
好きFortran
大規模並列計算
機械学習
まとめ

－169－

3

自己紹介
永井佑紀

専門：物性理論等

日本原子力研究開発機構システム計算科学センター
スパコンを使って、数千-数万CPUコアMPI並列計算を実行している
最近の研究分野：機械学習を物理学に使う
量子多体系へとか格子量子色力学へとか、分子動力学へとか
機械学習を「高速化」に使う話が好き（人間の代わりに頑張れ）
使えるプログラミング言語(得意な順)

Julia, Fortran2008以降, Python, c++

Juliaを知ったきっかけ
2016-2017(MIT客員研究員時代@米国ボストン)
同僚にJuliaを勧められたのでやってみる
ちょうどMITにいるしFortranと同じ配列が1始まりだし

便利なパッケージがたくさんあるし速いし -> Juliaの数値計算で論文を書く

4

自分野(物性物理学分野)に
おけるコード開発の状況
みんな大好きFortran

－170－

5

物性物理学とは
物質の性質を理論的に研究する分野：固体における電子状態の研究が多い
量子力学と統計力学、場の量子論などを駆使する
超伝導体、磁性体、半導体、量子コンピュータetc
現実に存在しうるあらゆるものが対象

more is different (多は異なり)
フィリップ・アンダーソン(1977年、ノーベル物理学賞)

物質中には電子が沢山いる->相互作用によって新しい性質が創発される（超伝導）
電子の数は1023個 ->そもそもまともに計算できない

いかに本質を抽出して計算可能にするかが大事
多彩すぎるため「一つの良い共通コード」がない

6

使われている言語
計算物理学：計算機とともに進展してきた分野
FORTRAN: (formula translation)
有名なFORTRAN

歴史が古い

1957年にコンパイラがリリース

FORTRAN77 1978年
1991-92年
Fortran90

45年前！
31年前！

数値計算でいまだに使われているプログラミング言語
Fortran90がほとんど。FORTRAN77もたまに見かける
大学学部のプログラミング実習でFortranを習う方も多い？

(汎用化できる領域での)汎用パッケージの多くはFortranで書かれている

例：密度汎関数理論に基づく第一原理計算パッケージ
Quantum Espresso
原子の数と種類と位置が与えられれば物性が計算できる
VASP
個別の研究でもFortranが多く使われている
WIEN2k

－171－

7

Fortranが使われている実例
少し古いが C. A. Perroni et al., Phys. Rev. B, 75, 045125 (2007)
“We employ the routine MINIMIZE provided in Ref. 1.”

[1] A. Georges et al., Rev. Mod. Phys. 68, 13 (1996)
動的平均場理論(DMFT)のレビュー論文。FORTRANコードが公開されている

MINIMIZEのコメント欄
(意訳)”このコードは長年使っ
てきた。残念ながら誰が書い
たかわからない。でも気にし
ないで、ちゃんと動くよ”

調べてみた

R.Fletcher, ”FORTRAN SUBROUTINES FOR MINIMIZATION BY QUASINEWTON METHODS", Research group report(United Kingdom Atomic Energy
Authority)

1972年4月のコード

8

この分野での最近の流れ
物性理論分野
ここ10年でc++の公開コードが増えてきた
小さめのコードはPythonで書かれている場合も
近年流行している機械学習+物理学ではFortranではきつい
機械学習系ライブラリはPythonが優位

少数でコードを書いて研究するグループは公開せずにFortranを使っている？

MPI並列計算する場合にはFortranはやはり便利

－172－

9

なぜFortranが(まだ)使われているのか
数百から数万のCPUコアを使ったMPI+Openmpのハイブリッド並列計
算をスーパーコンピュータで行いたい場合

スパコンごとに特別にチューニングされたコンパイラ：FortranとC++しかない

(C++11以降であればまだ良いが)C++よりFortranの方が配列の扱いがしやすい
キャッシュを有効活用し、MPI通信の待ち時間も利用したカリッカリのチューニングが可能
もちろんC++でも可能だが、スパコンではFortran人口が多い

C++と比べるとFortranの方が学習障壁が低い（個人の感想です）
数値計算ではオブジェクト指向はなくても良い
プログラミングよりも物理に着目したい
最初の約束事さえ覚えればすぐに研究
Fortranは数式に近い書き方ができる
先生がFortran(FORTRAN)使い

に取り掛かれる（と予想される）

研究室秘伝のコードがFortran

色々な理由があるが、結局「速い」が重要

「適当に書いても速い」もポイント

10

確かに速い

Fortranに思うところ

基本的には全てのパーツを自分で作らなければならない
LAPACKとBLASが使えるので固有値計算などは”簡単”にできる

例

“2次元シュレーディンガー方程式を円筒座標系で解きたいからBessel関数を使って基底を展開して、Bessel関数の零点を使って境界条件を設
定して対角化して固有値を求める”
「Bessel関数Jn(x)はどうやって計算するの？」
「LAPACKはどうやって使うの？」

“動的平均場理論で松原振動数表示のGreen関数を再現するような有効ハミルトニアンのパラメータを決めるために、非線形関数を
最小化したい”
「非線形最小化問題ってどうやって解くの？」
「非線形最小化問題を解くコードを書きたいんじゃなくて物理の問題を解きたいのに」
「持ってきたコードがFORTRANだけど自分のFortranコードとうまくリンクできない」

スパコンやクラスター計算機で並列計算をするならFortranは使う価値がある

でも、もう少し「軽い」問題なら？

－173－

JuliaならPython規模からFortran規模までカバーできる

11

大規模並列計算

12

大規模並列計算
並列計算とは
一度に複数の計算を行うこと
昨今のCPUであればマルチコアが当たり前

大規模並列計算とは

同時に様々なタスクが走っている

一度に大量の計算を行うこと

科学技術計算分野では
同時に”似た”計算を大量に並列に行う
非常に大量の計算を行なっているので、少しのロスも大きなロスになる
大規模並列計算特有の問題がある

－174－

13

並列計算の種類

大きく分けて二種類ある：スレッド並列とプロセス並列
スレッド並列

手が多い計算

プロセス並列

人が多い計算

共通の脳（メモリ）

それぞれがそれぞれ

にアクセスする

の仕事をする

高速に行うには：

高速に行うには：

それぞれのデータのやり取

いかに手がぶつからないよう
に手を動かすか。脳を同時書
き換えしないようにするか

りは時間がかかるので、通
信を減らす

メモリを共有する
並列計算

ハイブリッド並列

メモリを共有しない
並列計算

並列計算の種類

14

スレッド並列+プロセス並列
スーパーコンピュータ：
複数の物理サーバーが高速なLAN
ケーブルでつながっている
サーバー内：スレッド並列
サーバー間：プロセス並列

最先端物理学研究：スパコン上でFortranやC++でハイブリッド並列計算が行われている

－175－

その他、並列計算の種類
GPU計算

15

マルチGPU並列
特定のタスクでは
非常に強い
if文も含まない同じ計算の集まり
（行列・ベクトル積等）

機械学習などと相性が良い
A100観音

JuliaでマルチGPU：気になるトピックスですが、
自分もよくわかっていないので触れません

GPU:シンプルな演算器
がすごくたくさんある

16

Juliaでの並列計算
スレッド並列:後藤さんのtalk

プロセス並列:このtalk
MPI.jlを使うのが無難

物理の数値計算で「ハイパ
フォーマンスに」行う方法

他言語のMPI使いは
シンプルに移行できる

は正直よくわからない…

先行研究でOpenMPで出たパフォー
マンスと同じくらい出したい

ハイブリッド並列：今のところJuliaで「ハイパフォーマンスに」
行う方法を見つけられていない
->知りたい

－176－

17

Juliaでのプロセス並列
Juliaでのプロセス並列の方法は複数種類ある

1. Distributed.jl

やれ！

はい！

Juliaの標準のプロセス並列用パッケージ

単純で一つ一つが重たい並列計算ならこれが楽
各プロセスに投げたあと、通信はせず、後で結果
を回収するような数値計算

2. MPI.jl

-> pmapを使う

masterとworkerの並列モデル
指示する人とやる人が別

科学技術計算において“業界標準”であるMPIを
使えるパッケージ
スパコンが専用のMPIライブラリを持つ場合も使える
既存コードや論文と”同じように”Juliaで実装できる

同時に走り、それぞれを通信させる

18

Distributed.jlが効果的な例
Krylov部分空間法を用いた手法（共役勾配法等）
Krylov部分空間K(A,x): x,Ax A2x,A3x.,,,

初期ベクトルxが異なると異なるKrylov部分空間ができる

例：A x = biをiを変化させて解きたい
異なるbiは別々に計算できる
a= pmap(i -> goodfunc(i),1:N)のような形
aにはそれぞれのiで計算した値が代入される(要素数Nの配列)
もし各iでの計算が軽い場合は、スレッド並列の方が効率的（速い）

－177－

19

例：RSCG法(縮約シフト共役勾配法)

グリーン関数の行列要素Giiは、eiT (z -H)-1 eiで計算できるので、連立方程式
eiが異なるので異なるGiiはpmapで並列可能
(z-H) x = eiさえ解ければよい

https://github.com/cometscome/RSCG.jl

Julia パッケージ RSCG.jl

ポイント

なるべくJuliaの基本機能を使う
mul!とdot

これによって、mul!とdotが定

義されている型ならどんなもの
でも使える

疎行列を使う

MPI並列

20

MPI(Message Passing Interface)
1

-> CPUもメモリもそれぞれ別に持っている
それぞれのプロセスがそれぞれ自分の仕事をする
時々、他の人の結果をやりとりする

2

0

それぞれの船が時々手旗信号で情報を送り合うような

例

プロセス0番からプロセス1番にデータを渡せ
プロセス0番に全プロセスがデータaを集めろ

3

プロセス0番から全プロセスにデータaを渡せ
全プロセスのデータaを和をとり、結果を全プロセスに戻せ

母艦を指定してもよいしみんな計算してもよい

－178－

21

JuliaでのMPI

MPIはそれぞれのプロセスは独立に走るため、高速化をするには通信部分を

気をつける必要がある
-> 大規模並列計算の場合、並列化技法に関して、それぞれの専門分野で大量
の論文が存在する。それらを使うことができる

Distributed.jlでは、並列モデルが異なるため、自分で一から高速化
技法を考える必要がある -> 計算科学者でない人にとっては辛い

並列化手法を編み出して実装しても実績（論文）にならない…

「大規模並列計算可能な既存の良いFortranコードがあるなら、それを使えばいいじゃない」
->それはそう。
その道のプロならそのまま使えばいい。でも新規参入者はコンパイルできなかったり…
大規模コードの大規模並列化ではなく、自前の小規模コードの大規模並列化に向いている？
コードの拡張がしたい場合は、Juliaの方が読みやすい？

MPI.jlのインストール
1. 何らかの方法でMPIをインストールする

22

しなくてもMPI.jlが自動で入れてくれる

macならbrew install openmpiのようなもので入れられる

2. MPI.jlのドキュメントを見ながら入れる

通常はadd MPIで入るはず

クラスターやスパコンの場合
add MPIPreferencesでMPIPreferences.jlをインストール
using MPIPreferences
MPIPreferences.use̲system̲binary()

これで自分のMPIの環境を自動検出して設定してくれる

あとはadd MPIをするだけ

間違いなく対応しているMPI SGIなど対応していない場合は、OpenMPIを使用すると良いかもしれない

－179－

23

JuliaでのMPI
MPIの初期化と終了
MPIでは、自分が何番目のプロセスなのか、
何個プロセスがあるか、を知る必要がある
“全プロセスのデータwaを和をとり、結果を全プロセスに戻せ”

FortranやCと異なり、waがどんな型なのか、何個の要素を持つのか、ということを指示する必要がない
実行もFortranコードと同じように行えばよい

実行例

日本原子力研究開発機構

HPE SGI8600

GPU演算部：9.739PFLOPS
NVIDIA V100 1080基

Xeon Gold 24コアx2 x272-> 13056コア

CPU演算部:2.801PFLOPS
Xeon Gold 20コアx2 x706-> 28240コア
module loadでOpenMPIをロードし、MPIPreferencesで環境設定

－180－

24
mpirun julia test.jl
Hello world, I am 0 of 40
Hello world, I am 1 of 40
Hello world, I am 2 of 40
Hello world, I am 3 of 40
Hello world, I am 4 of 40
Hello world, I am 5 of 40
Hello world, I am 6 of 40
Hello world, I am 7 of 40
Hello world, I am 8 of 40
Hello world, I am 9 of 40
Hello world, I am 10 of 40
Hello world, I am 11 of 40
Hello world, I am 12 of 40
Hello world, I am 13 of 40
Hello world, I am 14 of 40
Hello world, I am 15 of 40
Hello world, I am 16 of 40
Hello world, I am 17 of 40
Hello world, I am 18 of 40
Hello world, I am 19 of 40
Hello world, I am 21 of 40
Hello world, I am 22 of 40
Hello world, I am 23 of 40
Hello world, I am 24 of 40
Hello world, I am 25 of 40
Hello world, I am 26 of 40
Hello world, I am 27 of 40
Hello world, I am 28 of 40
Hello world, I am 29 of 40
Hello world, I am 30 of 40
Hello world, I am 31 of 40
Hello world, I am 32 of 40
Hello world, I am 33 of 40
Hello world, I am 34 of 40
Hello world, I am 35 of 40
Hello world, I am 36 of 40
Hello world, I am 37 of 40
Hello world, I am 38 of 40
Hello world, I am 39 of 40
Hello world, I am 20 of 40

25

実行例

1次元強束縛模型 差分化された2階微分演算子D+対角要素V
H = D+V
eiT (z -H)-1 eiをRSCG.jlで計算
iのループを並列化

原子力機構のスパコン上で
MPI.jlを使用してみる
Julia 1.9.1で実行

26

実行例

1次元強束縛模型 差分化された2階微分演算子D+対角要素V
H = D+V
eiT (z -H)-1 eiをRSCG.jlで計算
iのループを並列化

ほぼ理想的なスケーリング
（ストロングスケーリング）
則を持つ

－181－

27

並列計算について

MPI.jlがこれまでの並列計算の知見を使えるため有益

スレッド並列＋プロセス並列の効率的な計算方法
はどうすればいいのか？
スレッド並列+プロセス並列をするときに何が問題か
スレッド並列はforループを並列化できる
-> 一つのforループを並列化する時に立ち上がりに若干のオーバーヘッドが存在する
-> 無数のforループを回す場合、無視できないほどのロスになる
FortranでのOpenMPでも同じ問題がある
->OpenMPの場合、スレッドを立ち上げたままにすることで回避できる

Juliaでのやり方はわからない……
OpenMPみたいにスレッド並列ができれ
ば先行研究の技法が使えるのになあ……

28

機械学習

－182－

29

機械学習

そもそも機械学習とは何か？
猫画像判別、自動運転、囲碁、自動お絵描きetc

何かを判断したり、推測したりを機械が行う -> 代わりにやってもらう

A

B

:人間には定式化できない何か

シミュレーション分野での機械学習の使用例

これまで計算が大変だった部分を機械が行う

A

B

定式化済みだけど計算コスト大

A

B
計算コストの軽い機械

シミュレーションの高速化に機械学習を用いる

30

教師あり学習

機械学習とは

猫画像その他

関数f(x)

y

猫?

猫です

yes or no

大量の入力データxを用いて、

学習プロセス

f(x)=yを満たす関数f(x)を決める

x

例：囲碁棋譜データ→勝利の方程式
囲碁プロ棋士に勝利

－183－

31

機械学習とは

教師あり学習

大量の入力データxを用いて、f(x)=yを満たす関数f(x)を決める

y = ax+ b

一番シンプルな例
インプットが複数→xがベクトル

直線で近似（線形回帰）

y = Wx+ b

もっと表現力を高めたい

x

W

y

y = W2 f(W1x+ b1)+b

f:非線形関数

y = W3f(W2 f(W1x+ b1)+b2)+b

：

深層学習

…
猫です

ニューラルネットワーク

32

昨今の機械学習

ChatGPT:生成AI
チャットの文章を入れるといい感じに返答してくれる
お絵描き系AI
文章で絵を説明すると、その文章にあった絵を描いてくれる
さまざまな形のネットワークが提案され続けている
進展が非常に速い…

-> 自分の研究に使いたかったらどうする？

－184－

33

機械学習とは

y = f(θ,x)

パラメータθとインプットxの関数f

“loss関数”L(θ,f)を最小化するθを見つける

最小化するには、パラメータ微分が必要
複雑な非線形関数をどうやって微分する？

-> 自動微分という技術がある

Fortranで機械学習をすることは可能か？
->可能。実際、機械学習分子動力学パッケージaenet等はFortranで書かれている
次々と出てくる新種の構造のネットワークの自動微分をその都度Fortranで実装する？？

Juliaで機械学習

34

自動微分の詳細について見る前に、Juliaで機械学習をする方法について
いくつか機械学習用パッケージがある:Flux.jlやLux.jlなど
->Flux.jlが一番有名
using Flux
function main()
model = Chain(Dense(2,10,relu),Dense(10,10,σ),Dense(10,1),sum)
x = Float32[0.2,0.3]
隠れ層2層の全結合ニューラルネットワーク
println(model(x))
dfdx = gradient(model,x)[1]
println(dfdx)
end
これだけで微分可能なニューラルネットワークが出来上がっている
main()

Chainで次々とレイヤーを繋げていく。普通の関数を間に挟んでもよい

グラディエントの計算ができる

－185－

35

カスタムレイヤー

全結合ではない、自分独自のネットワークを作りたい
using Flux
struct Dense2
W
b
構造体を作る
activation
end
Dense2(in::Integer,out::Integer) = Dense2(randn(out,in),randn(out),x -> x)
Dense2(in::Integer,out::Integer,activation) =
Dense2(randn(out,in),randn(out),activation)
(m::Dense2)(x) = m.activation.(m.W*x .+ m.b) インプットに対してどう応答するか定義する
Flux.trainable(a::Dense2) = (weight=a.W,bias=a.b) 訓練用パラメータを登録
Flux.@functor Dense2
Fluxでレイヤーとして使えるように登録

レイヤー内のそれぞれの計算が自動微分できるならば、標準のレイヤーと同等に使える
自動微分できない関数がある場合には、その関数用に自動微分を定義する（後述）

36

Flux.jlの使用例

“量子力学を機械学習で解く最新論文の結果(の一部)を再現したい”

̂ n(x) = Enψn(x)
Hψ

ℏ ∂2
+ V(x, y, z)
Ĥ = −
2m ∂x 2

En =

̂ n(x)
∫ dxψn(x)†Hψ
∫ dxψn(x)†ψn(x)

エネルギー固有値Enが最小となるような波動関数を機械学習で見つける

ψn(xi) = f(xi)

f(x):ニューラルネットワーク

A simple method for multi-body wave function of ground
and low-lying excited states using deep neural network
Tomoya Naito, Hisashi Naito, Koji Hashimoto, arXiv:2302.08965

やること

1. ハミルトニアンの微分演算子を差分化
2. エネルギーが最小となるような波動関数を探索

論文では、Pythonの機械学習フレームワークTensorFlowを使用しているらしい

－186－

37

Flux.jlの使用例

“量子力学を機械学習で解く最新論文の結果(の一部)を再現したい”
ハミルトニアンの微分演算子を差分化

using Flux
n1 = 4
model = Chain(x ->
[x],Dense(1,n1,Flux.softplus),Dense(n1,1,
Flux.softplus),x -> sum(x)) |> Flux.f64

ニューラルネットワークを構築
隠れ層1層、活性化関数softplusのニューラルネット

あとは訓練するだけ

調和振動子ポテンシャル項

運動エネルギー項

function energy(model,xs)
ψ = model.(xs)
c = ψ'*ψ
E = ψ'*H*ψ/c
return E
end

En =

̂ n(x)
∫ dxψn(x)†Hψ
∫ dxψn(x)†ψn(x)

PyTorchやTensorFlowと違い、
テンソルにする必要がない

エネルギーの定義

38

Flux.jlの使用例

“量子力学を機械学習で解く最新論文の結果(の一部)を再現したい”
Adamというオプティマイザーを使用
ハミルトニアンの微分演算子を差分化

opt = Flux.setup(Adam(), model)
for i=1:50000
energyを最小化するようにmodelを訓練
Flux.train!(energy, model,(xs,),
opt)
ニューラルネットワークを構築
e = energy(model,xs) modelを使ったenergyの計算
if i % 1000 == 0
調和振動子ポテンシャル項
println("i
= $i energy = $e")
function energy(model,xs)
end
̂ (x)
∫ dxψ (x) Hψ
ψ = model.(xs)
E =
∫ dxψ (x) ψ (x)
c = ψ'*ψ
end
これだけで最新の物理学with機械学習
E = ψ'*H*ψ/c
n

return E
訓練はこれだけ
end

運動エネルギー項

†

n

n

n

†

n

の結果を再現できる
エネルギーの定義

－187－

39

自動微分とは？
機械学習で登場する非線形関数の特徴：ブロック化されている

x1 = f1(x)
x2 = f2(x1) + g2(x1)
x3 = f3(x2)
インプットを順番に計算して
いくことで出力が得られる
y = f4(x3)
y = f4( f3( f2( f1(x)) + g2( f1(x))))

例 df/dxを計算したい
a(x) = sin(x)
連鎖律を使う
b(a) = 0.2 + a
∂f
∂c ∂a ∂c ∂c
∂b ∂c
=
=
=
c(b) = 1 + b
∂x ∂x
∂x ∂a ∂a ∂a ∂b
f(x) = c(b(a(x)))

∂c
∂c ∂c
=
∂b ∂b ∂c

誤差逆伝播法と呼ばれている

40

自動微分とは？

例 df/dxを計算したい
a(x) = sin(x)
連鎖律を使う
b(a) = 0.2 + a
∂f
∂c ∂a ∂c ∂c
∂b ∂c
=
=
=
c(b) = 1 + b
∂x ∂x
∂x ∂a ∂a ∂a ∂b
f(x) = c(b(a(x)))

誤差逆伝播と呼ばれている

∂c
∂c ∂c
=
∂b ∂b ∂c

微分を手で書ける
それぞれのブロック（あるいはレイヤー）の微分がわかれば計算できる
プログラミングでどうやって実現するの？
そのブロックがどのような計算をしているかを記録し、
その微分も計算しておけば良い？

－188－

自動微分とは？

41

プログラミングでどうやって実現するの？

そのブロックがどのような計算をしているかを記録し、
その微分も計算しておけば良い？

Python:機械学習分野ではデファクトスタンダードなプログラミング言語
機械学習パッケージ：Tensorﬂow、PyTorch etc.
Deﬁne by Runという技術で自動微分を実現
Tensorという型を導入し、その型を使った計算を記録しておく
-> 微分したい計算は、このTensorという型で表現しなければならない
一方、Juliaはどんな型も自動微分できる。なぜ？？

Juliaでの自動微分

42

Juliaでは、コードは動的にコンパイルされて実行されている
コンパイル：構文を解析し、最適なコードを生成する
Juliaは常に構文を解析している！
-> どのような計算が行われているか、把握している
計算の方法を記録する必要はない！ ->Tensorとかいらない！
だから、Juliaで書いたコードに対して、普通に微分できる
Zygoteという
パッケージを
使えば

function samplefunc(x)
a = sin(x)
b = 0.2+a
c = sqrt(1+b)
return c
end

x = 0.5
dfdx = gradient(samplefunc,x)[1]

－189－

これだけで微分できてしまう

Juliaでの自動微分
計算したい関数

その微分

43

微分計算に必要なもの
∂fi ∂aN
∂x ∂ai
∂fi
∂aN
はxの関数
は伝播してくる数字
∂x
∂ai
∂y
By(ȳ, x) ≡ ȳ 用意すると
∂x
∂a
∂aN
= By( N , ai−1)
∂ai−1
∂ai

と書ける

Juliaでの自動微分
計算したい関数

∂a
∂aN
= Bfi( N , ai−1)
∂ai−1
∂ai

44

各レイヤーiで Bfi(ȳ, x)

がわかれば自動微分できる

y = sin(x) なら By(ȳ, x) = cos(x)ȳ
->このBを「pullback」と呼ぶ
Juliaは構文解析によって基本的な演算の集まりに分割できている
cosやsin、足し算や掛け算など、基本的な演算にBを定義しておけばよい
多重ディスパッチによりrrule(::typeof(sin),x)のようなメソッドを定義する

ChainRulesCoreというパッケージにこれらの演算が定義されている
自作の型の場合は、自分でrruleを作れば微分できるようになる！

－190－

独自モデルの自動微分

45

インプットSに対してHを計算する関数
Ški =

α (S + Ski−l)
∑ l ki+l
l

H(S) = J

function calc_smearing(S,α)
nspin,nsite = size(S)
n = length(α)
Snew = deepcopy(S)
for i=1:nsite
for l=1:n
j = i+l
j += ifelse(j > nsite,-nsite,0)
for k=1:nspin
Snew[k,i] += α[l]*S[k,j]
end
j = i-l
j += ifelse(j <1,nsite,0)
for k=1:nspin
Snew[k,i] += α[l]*S[k,j]
end
end
end
return Snew
end

+ Ški−1)
Š (Š
∑ ∑ ki ki+1
i

k

レイヤーを定義
struct Smearing
α::Vector{Float64}
end
Smearing(n) = Smearing(randn(n))
(m::Smearing)(S) = calc_smearing(S,m.α)
Flux.trainable(a::Smearing) = (α=a.α,)
Flux.@functor Smearing

モデルを定義

独自モデルの自動微分

46

model = Chain(Smearing(n),S -> calc_hamiltonian(S,J))

function calc_hamiltonian(S,J)
nspin,nsite = size(S)
ham = 0.0
for i=1:nsite
j = i+1
j += ifelse(j > nsite,-nsite,0)
for k=1:nspin
ham += J*S[k,i]*S[k,j]
end
j = i-1
j += ifelse(j <1,nsite,0)
for k=1:nspin
ham += J*S[k,i]*S[k,j]
end
end
return ham
end

Ški =

α (S + Ski−l)
∑ l ki+l
l

H(S) = J

+ Ški−1)
Š (Š
∑ ∑ ki ki+1
i

k

モデルをパラメータとモデルに分ける
θ,re = Flux.destructure(model)
H = re(θ)(S)

パラメータ微分を計算したい
dHdθ_a = gradient(θ -> re(θ)(S),θ)[1]
for i=1:length(θ)
println("$i dHdθ = ",dHdθ_a[i])
end

－191－

47

独自モデルの自動微分
Ški =

α (S + Ski−l) をパラメータαで微分できるようにしたい
∑ l ki+l
function ChainRulesCore.rrule(::typeof(calc_smearing),S,α)
l

nspin,nsite = size(S)
n = length(α)
y = calc_smearing(S,α)
function pullback(ybar)
sbar = NoTangent()
fSbar = zero(S)
fabar = zero(α)

∂y
By(ȳ, x) ≡ ȳ を定義すればよい
∂x
ポイント：Bの型はxの型と同じ
∂f ∂Šnj
∑ ∂Š ∂Ski
nj
nj

BSˇki(S¯ki, αi) =

∂f ∂Šnj
∑ ∂Š ∂αi
を定義する
nj
nj

end

function ChainRulesCore.rrule(::typeof(calc_smearing),S,α)
nspin,nsite = size(S)
n = length(α)
y = calc_smearing(S,α)
function pullback(ybar)
sbar = NoTangent()
fSbar = zero(S)
fabar = zero(α)
for i=1:nsite
for l=1:n
j = i+l
j += ifelse(j > nsite,-nsite,0)
for k=1:nspin
fSbar[k, i] += ybar[k, j] * α[l]
fabar[l] += ybar[k, i]*S[k,j]
end
j = i-l
j += ifelse(j <1,nsite,0)
for k=1:nspin
fSbar[k, i] += ybar[k, j] * α[l]
fabar[l] += ybar[k, i]*S[k,j]
end
end
end
return sbar,fSbar,fabar

end

BSˇki(S¯ki, Ski) =

∂f
α
∑ ∑ ∂Š l
n

nj

l

∂f

∑ ∑ ∂Š
n

l

ni−l

Ski

for i=1:nsite
for l=1:n
j = i+l
j += ifelse(j > nsite,-nsite,0)
for k=1:nspin
fSbar[k, i] += ybar[k, j] * α[l]
fabar[l] += ybar[k, i]*S[k,j]
end
j = i-l
j += ifelse(j <1,nsite,0)
for k=1:nspin
fSbar[k, i] += ybar[k, j] * α[l]
fabar[l] += ybar[k, i]*S[k,j]
end
end
end
return sbar,fSbar,fabar

BSˇki(S¯ki, Ski) と BS̄ki(S¯ki, αi) が必要
BSˇki(S¯ki, Ski) =

BSˇki(S¯ki, Ski) =

end
return y, pullback

－192－

end
return y, pullback

Ški =

α (S + Ski−l)
∑ l ki+l
l

BSˇki(S¯ki, Ski) =
BSˇki(S¯ki, Ski) =

48

∂f
α
∑ ∑ ∂Š l
n

l

nj

n

l

ni−l

∂f
Ski
∑ ∑ ∂Š

rruleにpullbackを返す関数
を実装すればよい
dHdθ_a = gradient(θ -> re(θ)(S),θ)[1]
for i=1:length(θ)
println("$i dHdθ = ",dHdθ_a[i])
end

あとはgradientで自動微分できる

独自モデルの自動微分

インプットSに対してHを計算する関数
Ški =

α (S + Ski−l)
∑ l ki+l
l

H(S) = J

先程のrruleを定義するだけで、
∂H(S)
∂H(S)
も
∂Sik
∂αi

+ Ški−1)
Š (Š
∑ ∑ ki ki+1
i

k

49

も計算できる

レイヤーを複数重ねても自動微分ができる
model = Chain(Smearing(2),Smearing(3),Smearing(2),S -> calc_hamiltonian(S,J))

応用例
Attentionレイヤーを
Flux.jlで実装し同変

transformerを実装した

YN and A. Tomiya, arXiv:2306.11527

50

Juliaで機械学習
シンプルに書いて機械学習ができる
どんな型でも自動微分できる -> あとで自動微分を実装することすらできる
rruleはある関数の微分を計算できるようになる
->あえて書くことで高速に計算することも可能

(連鎖律を短くすることができる)

GPUも使える

AMD GPUやMacのMetalもできるようになってきている

θ,re = Flux.destructure(model)
H = re(θ)(S)
パラメータとモデルにバラすと、色々なパッケージで最適化可能

－193－

51

まとめ

52

Juliaでの科学技術計算
並列計算：MPI.jlを使うとよい
他の言語でMPIを使っていた人は簡単に移行可能
MPIはじめてでも、他の言語よりやりやすい
FortranやCのMPIの資料を見ながら勉強し、そのままJuliaで実装できる

機械学習はFlux.jlで便利に使える
Juliaは構文を解析し、コンパイルしている
構文を解析してくれているため、自動微分との相性がよい
Tensorがいらない。どんな型でも自動微分できる

－194－

数学と物理におけるJuliaの活用

https://akio-tomiya.github.io/julia_imi_workshop2023/
https://joint1.imi.kyushu-u.ac.jp/research_chooses/view/2023a015
7/12(水)

IMI in person

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

素粒子物理におけるJuliaの活用
〜格子QCDの大規模計算に向けて〜
Akio Tomiya (IPUT Osaka, Assistant Prof.)
akio_at_yukawa.kyoto-u.ac.jp

Outline
• 素粒子物理学? 格子QCD?
• Julia言語での取り組み
• 何ができれば良いのか

2

－195－

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

How to treat gauge ﬁelds
with neural networks?

3

Applications of neural network in natural science Akio Tomiya
Neural network have been good job

Folding of protein (AlphaFold2, John Jumper+, Nature, 2020+)
Score:
Higher is better

Neural network wave function for many body (Carleo Troyer, Science 355, 602 (2017) )

Variational energy
(lower is better)
# of units ∝

Neural net + Expert knowledge → Best performance

https://horomary.hatenablog.com/entry/2021/10/01/194825

4

https://www.aist.go.jp/aist_j/press_release/pr2020/pr20201111/pr20201111.html

－196－

What is neural networks?

Akio Tomiya

Aﬃne transformation + element-wise transformation
Layers of neural nets l = 2,3,⋯, L, u (1)
⃗ = x⃗

⃗
⃗
z (l)
⃗ = W (l) u (l−1)
+ b (l)

{ u (l) = σ (l)(z (l))
i

Aﬃne transf.

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

(b=0 called linear transf.)

element-wise (local)

i

A fully connected neural net

⃗ )
⃗ ) + b (3)
fθ( x ⃗) = σ (3)(W (3)σ (2)(W (2) x ⃗ + b (2)
θ is a set of parameters: wij(l), bi(l), ⋯

Neural network = map between vectors and vectors
Physicists terminology: Variational transformations
5

What is the neural networks?

Akio Tomiya

Neural network is a universal approximator
Example: Recognition of hand-written numbers
6x6

36 dimension

x ⃗ Input

0.000
0.000
0.8434
0.756
=
0.3456
0.64
0.251
⋮

= x⃗

Neural net

Image is a vector
(6x6=36 dim)

1
10 dimension

Images of “1”

W (2) W (3)

2

“0” = (1,0,0,…)
“1” = (0,1,0,…)
“2” = (0,0,1,…)
…
“9” = (0,0,…,1)

Images of “2”

Fact: Neural network can mimic any function
= A systematic variational function.
In this example, NN mimics image (36-dim vector) and label (10-dim vector)
6

－197－

ML for LQCD is needed
• Neural networks
• Data processing techniques mainly for 2d
image (a picture = pixels = a set of real #)

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

• Neural network helps data processing
e.g. AlphaFold2

• Lattice QCD requires numerical eﬀort

but is more complicated than pictures

• 4 dimension

• Non-abelian gauge d.o.f. and symmetry
• Fermions (Fermi-Dirac statistics)

• Exactness of algorithm is necessary
• Q. How can we deal with neural nets?

http://www.physics.adelaide.edu.au/theory/staﬀ/leinweber/VisualQCD/QCDvacuum/

7

Gauge covariant neural network
= trainable smearing

Akio Tomiya

AT Y. Nagai arXiv: 2103.11965

Smearing = Gauge covariant way of transform gauge conﬁgurations
Covariant sum

Staple
α
̂ ν†(n + μ)̂ + ⋯
Uμ(n) → Uμsmr(n) = 𝒩𝒩 (1 − α)Uμ(n) + Vμ†[U](n) Vμ†[U ](n) = ∑ Uν(n)Uμ(n + ν)U
[
]
6
μ≠ν

𝒩𝒩 [M] =

M

Normalization or projection

M †M

Gauge covariant neural network = General smearing with tunable parameters w

zμ(l)(n) = w1(l)Uμ(l−1)(n) + w2(l)𝒢𝒢(l)[U]

{𝒩𝒩(zμ(l)(n))

θ̄

Train (tune, ﬁtting)

point-wise (local)

Gauge covariant NN: UμNN(n)[U] = Uμ(4)(n)[Uμ(3)(n)[Uμ(2)(n)[Uμ(n)]]]

NN

Gauge covariant variational map: Uμ(n) ↦ Uμ

(n) = UμNN(n)[U]

Alternative realization of gauge symmetric neural net: gauge equivariant neural net -> Tilo Wettig’s talk

－198－

8

Gauge covariant neural network

Akio Tomiya

Schematic illustrations for neural networks (NN)

Neural networks for images

Parametrized function

チュートリアル講演
Julia

Ans
Dog

Cat

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

Dense NN
Parameters θ

Convolution
Parameters θ

Tune by backprop (training)

Neural networks for gauge conﬁgurations
Parametrized function
AT Nagai 2103.11965

Wilson loop

U

Covariant NN
(Tunable smearing)
Parameters θ

U NN[U]

SD[U
Dirac op.
(Functional
of conﬁgurations)

“Ans”

SD[U]

NN

[U]]

Tune by backprop (training)
http://www.physics.adelaide.edu.au/theory/staﬀ/leinweber/VisualQCD/QCDvacuum/

9

Alternative realization of gauge symmetric neural net: gauge equivariant neural net -> Tilo Wettig’s talk

Gauge covariant neural network
= trainable smearing

AT Y. Nagai arXiv: 2103.11965

(convolutional)
Neural network

Gauge Covariant
Neural network

Input

Image
(2d data, structured)

gauge conﬁg
(4d data, structured)

Output

Image
(2d data, structured)

gauge conﬁg
(4d data, structured)

Symmetry

Translation

with Fixed param

Image ﬁlter

Translation, rotation(90°),
Gauge sym.
Gauge sym
(APE/stout …) Smearing

Local operation

Summing up nearest
neighbor with weights

Summing up staples
with weights

Activation function

Tanh, ReLU, sigmoid, …

projection/normalization
in Stout/HYP/HISQ

Formula for chain rule

Backprop

“Smeared force
calculations” (Stout)

Training?

Backprop + Delta rule

AT Nagai 2103.11965

Dictionary

Akio Tomiya

Well-known

(Index i in the neural net corresponds to n & μ in smearing. Information processing with NN is evolution of scalar ﬁeld)
10

－199－

SLHMC = Exact algorithm with ML

Akio Tomiya

SLHMC for gauge system with dynamical fermions

Gauge covariant neural network can mimics gauge invariant functions
-> It can be used in simulation? -> Self learning HMC!

HMC

U

U

U

U

U

U
Eom

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

U

Eom

U

U′
Eom

U

U

Eom

U

ϕ

Metropolis

ϕ

U

π

π′

π′
ϕ

Metropolis

Eom

Eom

Eom

U′

ϕ

Eom

U

Gauss

Self
Learning
HMC

π

Eom

Gauss

U

Eom

Eom

1 / 87

Metropolis

Both use

HHMC =

1
π 2 + Sg + Sf
2∑

Non-conservation of H cancels since
the molecular dynamics is reversible

Metropolis
1
H=
π 2 + Sg + Sf[U]
2∑
Eom

1
π 2 + Sg+Sf[U NN[U]]
2∑
Neural net approximated
fermion action but exact

H=

SLHMC works as an adaptive reweighting!

arXiv: 2103.11965
and reference therein

Application for the staggered in 4d
Problems to solve

11

Akio Tomiya

arXiv: 2103.11965

Mimic diﬀerent action (Dirac operator):

Target action
(Metropolis)

S[U] = Sg[U] + Sf[ϕ, U; m = 0.3],

U

U

U

U

Eom

(Final target: Domain-wall vs overlap)
A toy problem: Staggered (heavy) vs Staggered (light)

U

Eom

U′

ϕ

π′
ϕ

Metropolis

U

Eom

U

Eom

π

Sθ[U] = Sg[U] + Sf[ϕ, UθNN[U]; mh = 0.4],

Eom

Self
Learning
HMC

Gauss

{Action in MD

Metropolis
1
H=
π 2 + Sg + Sf[U]
2∑
Eom

1
π 2 + Sg+Sf[U NN[U]]
2∑
Neural net approximated
fermion action but exact

H=

SLHMC works as an adaptive reweighting!
－200－

mimic

12

Application for the staggered in 4d
Results are consistent with each other

Akio Tomiya

arXiv: 2103.11965

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

Expectation value

Implemented by
13

Julia language on Fugaku,
Lattice code in Julia
(LatticeQCD.jl and GomalizingFlow.jl)

14

－201－

Benchmark

Akio Tomiya

Speed of Julia ~ Clang
1.Open source scientiﬁc language (Just in time compiler/LLVM backend)
2.Fast as C/Fortran (faster sometimes), Practical as Python
3.Machine learning friendly
チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

Compiler

Benchmark
(sec)
single core

Type

Parallelism

GPU

Pros👍👍

Cons👎👎

Column,
row

Julia
(1.8)

JIT,
LLVM

0.0014 Dynamic
& Static

MPI,
others

CUDA

Fast
Practical
ML feiendly

not major

columnmajor

C

Clang
(LLVM)

0.0033

Static

MPI,
others

CUDA

Fast

Long
codes

row-major

Python
+Numba

(CPython)
JIT, LLVM

0.0131 Dynamic

Available

Numba- Practical
Not fully row-major
CUDA ML feiendly supported (Numpy)

Note

(Rosetta2
is used in
benchmark)

C and Julia have similar speed
Benchmarks are performed on m1 mac mini (similar tendency on Xeon)
Benchmark: Multiplications for 12dim vector and 12x12 complex matrix for 10^4 times (repeated 10 times)
15

Benchmark

15

Akio Tomiya

Code comparison
using Random
function main()
T = 10
K = 10^4
N = 12
#
A = zeros(ComplexF64, (N,N))
V = zeros(ComplexF64, N)
W = zeros(ComplexF64, N)
function myprod(A,V,W)
for k = 1:N
for i = 1:N
W[i] += A[i, k]*V[k]
end
end
end
…(cut)…
Attached in backup

#include <stdio.h>
#include <complex.h>
#include <math.h>
#include <time.h>
#include <stdlib.h>
#define T 10
#define K 10000
#define N 12
…(cut)…
void myprod(double complex A[N][N], double complex *V,
double complex *W) {
for (int k = 0; k < N; k++) {
for (int i = 0; i < N; i++) {
W[i] += V[k] * A[k][i];
}
}
}
…(cut)…
Attached in backup

• Complex matrix (12x12) times complex vector (d=12)
• One set= 10^4 times, and repeated 10 times and averaged
• Code of Julia looks like Python (short, simple) but fast as C

Julia: 0.0014 (sec), C: 0.0033 (sec). Single core performance is similar

I thank Taku Izubuchi

16

－202－

16

Benchmark

Akio Tomiya

Why Julia? (My personal opinion)

[1] https://akio-tomiya.github.io/julia_in_physics/
[2] https://qr.ae/prgSG5

• Modern scientiﬁc programming language

• Easy to make codes. Fast as C/C++ (Julia& C use LLVM)
• Fewer compiling/dependency issues.

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

• Many people are potentially interested in. (More than 400 people registered to
“Julia in physics 2022 online workshop” [1]). 4,923 public repo on Github

• No two Language problem. “The fact that while the users are programming in
a high-level language such as R and Python, the performance-critical parts
have to be rewritten in C/C++ for performance”. [2]

• Neural network friendly (Flux.jl). Tensor networks also (iTensor.jl).

• Works on/with

• Xeon, Radeon/Apple silicon/A64FX
• MPI, GPU

17

17

LLVM?

Akio Tomiya

LLVM = common backend for making binaries on multiple architectures
Source

Execution
Frontend

const K = 5
const N = 1000
A = zeros(Float64, (N,N))
B = zeros(Float64, (N,N))
C = zeros(Float64, (N,N))

Julia Compiler
(Just in time)
Optimizing codes

LLVM-IR

Frontend
#include <stdio.h>
#define K 5
#define N 1000
double A[N][N];
double B[N][N];
double C[N][N];

x86

LLVM-optimizer&
LLVM-backend
Binary

LLVM
m1, m2

“Clang” Compiler
(Ahead of Time)
Optimizing codes

https://www.fujitsu.com/jp/about/businesspolicy/tech/fugaku/
https://ja.wikipedia.org/wiki/Apple_M1
https://ja.wikipedia.org/wiki/Xeon
https://ja.wikipedia.org/wiki/Ryzen
https://gigazine.net/news/20200623-japan-fugaku-fastest-supercomputer/

arm

used in Fugaku
See: https://en.wikipedia.org/wiki/LLVM and related pages

18

－203－

Julia is ready on Fugaku(?)
Parallelization with A64FX/Fugaku
Tests of MPI +

M. Giordano, arXiv:2207.12762v1 [cs.DC] 26 Jul 2022

on Fugaku

Send-Recv performance
チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

Fig. 2. Comparison of latency (top panel) and throughput (bottom panel)
of inter-node point-to-point MPI communication between using MPI.jl in
Julia and IMB benchmarks in C (results provided by R-CCS in [18]). Fugaku
scheduler setup: -L "node=2" -mpi "max-proc-per-node=1".

has similar scaling of MPI with C
(no obvious overhead)

19
MPI ranks across 384 nodes using the torus layout, to match the
scheduler configuration of the R-CCS benchmarks. MPI.jl
typically showed very small overhead for messages larger than
Akio Tomiya
1-2 KiB—peak throughput of ping-pong communication with
MPI.jl is within 1% of that reported by R-CCS—, but slightly
larger overhead for messages of smaller sizes. We note that,
AT & Y. Nagai in prep
contrary
to IMB, of
at latency
the present
timeandMPIBenchmarks.jl
Fig. 2. Comparison
(top panel)
throughput (bottom panel)
of inter-node
point-to-point
MPI communicationmechanism,
between usingwhich
MPI.jl
in
does
not implement
a cache-avoidance
may
Open
source
Julia and IMB benchmarks in C (results provided by R-CCS
in [18]).
Fugaku(Julia Oﬃcial package, Now updated to v1.0)
explain
appears
to show
better latency than IMB
schedulerwhy
setup:MPI.jl
-L "node=2"
-mpi
"max-proc-per-node=1".
for messages with size up to 64 KiB, which corresponds to
the size ofMachines:
the L1 cache Laptop/desktop/Jupyter/Supercomputers
of the A64FX CPU. We also observe
MPIcontrary
ranks across
384 we
nodes
the torus
layout, to
match the
that,
to [16],
didusing
not find
a significant
performance
Functions:
SU(Nc)-heatbath,
(R)HMC,
scheduler
configuration
of operation
the R-CCS
MPI.jl
drop
for the
Allreduce
forbenchmarks.
larger message
sizes.Self-learning HMC, SU(Nc) Stout
typically showed very small
overhead forStaggered,
messages largerDynamical
than
Dynamical
Wilson,
Dynamical Domain-wall
Fig. 3.
Comparison of latency of collective MPI operations
B.1-2Type
flexibility
and reduced-precision
with Float16with
KiB—peak
throughput
of
ping-pong
communication
between using MPI.jl in Julia and IMB benchmarks in C (results
Measurements
MPI.jl
is within
1% of applications
that reported by
R-CCS—,
but slightly
Developing
complex
using
Float16
is not provided by R-CCS in [18]): MPI Allreduce (top panel), MPI
larger
of smaller
sizes. We
note that, Gatherv (middle panel), MPI Reduce (bottom panel). Fugaku
easy.
Onoverhead
A64FX, for
evenmessages
the occasional
occurrence
of
subnormals
1. Download Julia scheduler
binary setup: -L "node=4x6x16:torus:strict-io" -L
8the present5 time MPIBenchmarks.jl
toStart
IMB,
Float16
(6 · 10at
to 6 · 10 ) causes a heavy performance "rscgrp=small-torus" -mpi proc=1536.
ofcontrary
LQCD
2.
Add
the
package
through Julia package manager
9
does not
a cache-avoidance
mechanism,
which
may
penalty
butimplement
a compiler-flag
is set to flush
them to zero
instead
in 5 min
explain
why MPI.jl
appears
show better
than
3.6 Execute!
The
available
normal range
of to
Float16,
·latency
10 5 to
65,IMB
504,
for messages with size up to 64 KiB, which corresponds to
9 https://github.com/JuliaLang/julia/issues/40151.
the
size of the L1 cache of thehttps://github.com/akio-tomiya/LatticeQCD.jl
A64FX CPU. We also observe
that, contrary to [16], we did not find a significant performance
drop for the Allreduce operation for larger message sizes.

Lattice QCD code

Open source LQCD code in Julia Language

B. Type flexibility and reduced-precision with Float16
Developing complex applications using Float16 is not
easy. On A64FX, even the occasional occurrence of subnormals
of Float16 (6 · 10 8 to 6 · 10 5 ) causes a heavy performance
penalty but a compiler-flag is set to flush them to zero instead9
The available normal range of Float16, 6 · 10 5 to 65, 504,

Fig. 3.
Comparison of latency of collective MPI operations
between using MPI.jl in Julia and IMB benchmarks in C (results
provided by R-CCS in [18]): MPI Allreduce (top panel), MPI
Gatherv (middle panel), MPI Reduce (bottom panel). Fugaku
-L "node=4x6x16:torus:strict-io" -L
scheduler
setup:
"rscgrp=small-torus" -mpi proc=1536.

9 https://github.com/JuliaLang/julia/issues/40151.

20

－204－

Package structure

Akio Tomiya

Our lattice QCD codes are constructed by following repositories
Dependency (Automatically solved)

QCDMeasurements.jl

Wrapper for LatticeDiracOperators.jl &
Gaugeﬁelds.jl, QCDMeasurements.jl
- Wizard for parameter ﬁles
- HMC/RHMC for SU(Nc)
- Stout + Wilson/Staggered/DW
- Heatbath for SU(Nc)
- Measurements
- Jupyter, Colab/PC/Supercomputers
etc

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

LatticeDiracOperators.jl
Gaugeﬁelds.jl
Wilsonloop.jl

Measurements in LQCD
(Correlator, Flow, Qtop, etc)

CLIME_jll

Fermions (+HMC), Wilson, KS, DW, MPI
PC/Supercomputers

ILDG I/O

Gauge ﬁelds (+HMC/Heatbath), MPI
PC/Supercomputers

Symbolic operations of Wilson/Polyakov loops
See https://github.com/akio-tomiya/LatticeQCD.jl in detail

21

Benchmark of Julia + QCD

Akio Tomiya

Wilson inversion / MPI parallel, Strong Scaling

AT & Y. Nagai in prep
AT & Y. Nagai in prep

Tested on Yukawa-21@YITP

Absolute execution time

Relative speed up

LatticeDiracOperators.jl

# procs

1

2

4

16

32

LatticeDiracOperators.jl

64

# procs

1

2

4

16

32

64

It looks scaling well

We need more contributors!
Please help us

We thank to H. Ohno & Issaku Kanamori

22

－205－

22
Akio Tomiya

Flow based sampling algorithm

Akio Tomiya

Trivializing map realized using neural network

Normalizing ﬂow? = Trivializing map, exact MCMC with ML

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

Change of variable by a neural network (Normalizing ﬂow)

ϕ = F (NN)[φ]

Sampling from Gaussian
→ Inverse trivializing map (neural net)
→ QFT conﬁgurations → Tractable Jacobian (by even-odd strategy)
→ After sampling, Metropolis-Hastings test → exact!
Related talk: Nobuyuki Matsumoto (Feb 16)

arxiv 1904.12072, 2003.06413, 2008.05456

23

Normalizing ﬂow in Julia

Akio Tomiya

A public code in Julia Language

https://arxiv.org/abs/2208.08903

A public Julia code for the ﬂow-based
sampling algorithm for scalar ﬁeld.
This supports not only 2d but also 3d.
CPU/GPU with Docker.

A new type of convolution
improves acceptance rate
(~ shorten the autocorrelation)
I reported in NeurIPS 2022 workshop
https://ml4physicalsciences.github.io/2022/

https://github.com/AtelierArith/GomalizingFlow.jl

24

－206－

Summary

Akio Tomiya

LQCD + Machine learning by Julia language

arXiv: 2103.11965
arXiv: 2208.08903

• Machine learning for LQCD

• Neural net (NN) + expert knowledge -> Best performance
e.g. AlphaFold2, NN wave functions

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

• NN can deal with 4d non-abelian gauge symmetric data now
• Self-learning HMC with NN works for dynamical fermions

• Julia language for LQCD/HPC/ML

• Julia has similar speed with C (w/ & w/o MPI), and machine learning friendly
• Two Julia codes for lattice ﬁeld theory

• LatticeQCD.jl: A suite lattice QCD code, machine learning
• GomalizingFlow.jl: Trivializing map via a neural network
25
KAKENHI: 20K14479, 22H05112,
22H05111, 22K03539

26

－207－

Thanks!

25

MPI benchmark

Akio Tomiya

PingPong

M. Giordano, arXiv:2207.12762v1 [cs.DC] 26 Jul 2022

チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

This section describes the performance of PingPong as a oneto-one communication. PingPong sends data between two
ranks by Send communication from one rank, receives the data
at the other rank, and then sends the data back to the original
rank by Recv communication.

Fig. 2. Comparison of latency (top panel) and throughput (bottom panel)
of inter-node point-to-point MPI communication between using MPI.jl in
Julia and IMB benchmarks in C (results provided by R-CCS in [18]). Fugaku
scheduler setup: -L "node=2" -mpi "max-proc-per-node=1".

Allreduce collects data from all ranks, performs a set
operation, and transfers the result to and transfers
the result to all ranks.

https://jopss.jaea.go.jp/pdfdata/JAERI-Data-Code-2001-010.pdf
27

Introduction

MPI ranks across 384 nodes using the torus layout, to match the
scheduler configuration of the R-CCS benchmarks. MPI.jl
typically showed very small overhead for messages larger than
1-2 KiB—peak throughput of ping-pong communication with
MPI.jl is within 1% of that reported by R-CCS—, but slightly
larger overhead for messages of smaller sizes. We note that,
contrary to IMB, at the present time MPIBenchmarks.jl
does not implement a cache-avoidance mechanism, which may
explain why MPI.jl appears to show better latency than IMB
for messages with size up to 64 KiB, which corresponds to
the size of the L1 cache of the A64FX CPU. We also observe
that, contrary to [16], we did not find a significant performance
drop for the Allreduce operation for larger message sizes.

Akio Tomiya
Fig. 3.
B. Type flexibility and reduced-precision with Float16
Developing complex applications using Float16 is not
easy. On A64FX, even the occasional occurrence of subnormals
of Float16 (6 · 10 8 to 6 · 10 5 ) causes a heavy performance
penalty but a compiler-flag is set to flush them to zero instead9
The available normal
range of Float16, 6 Reference
· 10 5 to 65, 504,
Exact?
Fermion?

Conﬁguration generation with machine learning is developing
Year

Group

2017

AT, Akinori
Tanaka

2018

ML

RBM
+ HMC

Dim.
2d

Theory
Scalar

Gauge sym

-

Fig. 2. Comparison of latency (top panel) and throughput (bottom panel)

K.
Zhou+ point-to-point
Scalar
- in
of inter-node
communication
between using MPI.jl
GAN MPI2d

Julia and IMB benchmarks in C (results provided by R-CCS in [18]). Fugaku

GAN
Pawlowskisetup:
+
2018 J.scheduler
-L "node=2"2d
-mpi "max-proc-per-node=1".
Scalar
+HMC

MIT+
2d
Scalar
Flow
MPI ranks across 384 nodes using the torus layout, to match the
2020
MIT+
2d
Equivariant
U(1)
Flow
scheduler configuration of the R-CCS benchmarks. MPI.jl
very small
for messages
larger than
2020 typically
MIT+ showed
2doverheadSU(N)
Equivariant
Flow
1-2
KiB—peak throughput of ping-pong communication with
AT,
Akinori
2020 MPI.jl is within
4dreportedSU(N)
SLMC
Invariant
1% of that
by R-CCS—,
but slightly
Tanaka +
overhead
for messages
sizes. We note that,
Medvidovic´+
2021 M.larger
2d of smaller
A-NICE
Scalar
contrary to IMB, at the present time MPIBenchmarks.jl
Foreman
2021 S.does
2d
L2HMCa cache-avoidance
not implement
mechanism, which
U(1)
Yesmay
explain why MPI.jl appears to show better latency than IMB
2021
AT+
Covariant
SLHMC
QCD
for messages with size up4d
to 64 KiB,
which corresponds to
Del of the L1 cache of2d
the
A64FX
CPU.
We
also observe
Scalar,
O(N)
2021 theL.size
Flow
Debbio+
that,
contrary to [16], we did not find a significant performance
2021 drop
MIT+
Yukawa
for the Flow
Allreduce2d
operation
for larger message- sizes.
2019

27

Flowed
2021 S.B.Foreman,
2d
Equivariant
U(1)with Float16
Type flexibility
and reduced-precision
AT+
HMC

Neural
complex applications
using Float16
is not
2021 XYDeveloping
Jing
2d
Equivariant
U(1)
neteven the occasional occurrence of subnormals
easy. On A64FX,
Flow
2022 J.ofFinkenrath
U(1)
Float16 (6
· 10 8 to 62d
· 10 5 ) causes
a heavyEquivariant
performance
penalty
but
a
compiler-flag
is
set
to
flush
them
zero instead9
MIT+
Flow
2022
2d, 4d U(1), QCD toEquivariant
5
The available normal range of Float16, 6 · 10 to 65, 504,
AT+
Flow
Scalar
2022
2d, 3d
9

No

No

Yes?

No

9 https://github.com/JuliaLang/julia/issues/40151.

arXiv: 1712.03893

No

arXiv: 1810.12879

Yes
Yes
Yes
Yes
No

No
No
No
No

Partially

arXiv: 1811.03533
arXiv: 1904.12072
arXiv: 2003.06413
arXiv: 2008.05456
arXiv: 2010.11900

No

arXiv: 2012.01442

Yes

YES!

This talk

Yes

Yes

Yes

No

Yes

Yes

Yes
Yes

No
No

Compariso
between using MPI.j
provided by R-CCS
Gatherv (middle pan
scheduler
setup:
"rscgrp=small-tor

Fig. 3.
Comparison
latency of collective
operations
arXiv:MPI
2112.01586
No butofcompatible
between using MPI.jl in Julia and IMB benchmarks in C (results
provided by R-CCS in [18]): MPI Allreduce (top panel), MPI
Gatherv (middle panel), MPI Reduce (bottom panel). Fugaku
scheduler
setup:
-L "node=4x6x16:torus:strict-io"
-L
arxiv: 2201.02216
Yes (diagonalization)
"rscgrp=small-torus" -mpi proc=1536.

https://github.com/JuliaLang/julia/issues/40151.

Yes
Yes
Yrs

arXiv:2202.11712 +

＋…
28

－208－

Benchmark

Akio Tomiya

Code comparison
using Random
function main()
T = 10
K = 10^4
N = 12
#
A = zeros(ComplexF64, (N,N))
V = zeros(ComplexF64, N)
W = zeros(ComplexF64, N)
function myprod(A,V,W)
for k = 1:N
for i = 1:N
W[i] += A[i, k]*V[k]
end
end
end
function test(A,V,W)
for jj=1:T
runtimes=[]
for r=1:K
A .= rand(N,N) + im*rand(N,N)
V .= rand(N) + im*rand(N)
W .= 0
tmp = @elapsed myprod(A,V,W)
push!(runtimes,tmp)
end
println("$(jj-1) $(sum(runtimes)) #W[1] = $(W[1])")
end
end
test(A,V,W)
end
if abspath(PROGRAM_FILE) == @__FILE__
main()
end

#include <stdio.h>
#include <complex.h>
#include <math.h>
#include <time.h>
#include <stdlib.h>
#define T 10
#define K 10000
#define N 12
double urand(){
double m, a;
m = RAND_MAX + 1;
a = (rand() + 0.5)/m;
a = (rand() + a)/m;
return (rand() + a)/m;
}
チュートリアル講演
Julia

⼊⾨

Satoshi Terasaki (AtelierArith)
2023/07/10

1 / 87

void myprod(double complex A[N][N], double complex *V, double complex *W) {
for (int k = 0; k < N; k++) {
for (int i = 0; i < N; i++) {
W[i] += V[k] * A[k][i];
}
}
}
void test(double complex A[N][N], double complex *V, double complex *W) {
for (int jj = 0; jj < T; jj++) {
double runtimes = 0;
for (int r = 0; r < K; r++) {
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) {
A[i][j] = urand() + urand() * I;
}
V[i] = urand() + urand() * I;
W[i] = 0.0 + 0.0 * I;
}
clock_t start = clock();
myprod(A, V, W);
clock_t end = clock();
runtimes += (double)(end - start) / CLOCKS_PER_SEC;
}
printf("%d %f # W[0] = %f %f\n", jj, runtimes, creal(W[0]), cimag(W[0]) );
}
}
int main() {
double complex A[N][N];
double complex V[N];
double complex W[N];
test(A, V, W);
}

return 0;

29

－209－

29

「マス・フォア・インダストリ研究」シリーズ刊行にあたり
本シリーズは，平成 23 年 4 月に設立された九州大学マス・フォア・インダストリ研究所
(IMI)が，平成 25 年 4 月に共同利用・共同研究拠点「産業数学の先進的・基礎的共同研究
拠点」として，文部科学大臣より認定を受けたことにともない刊行するものである．本シ
リーズでは，主として，マス・フォア・インダストリに関する研究集会の会議録，共同研
究の成果報告等を出版する．各巻はマス・フォア・インダストリの最新の研究成果に加え，
その新たな視点からのサーベイ及びレビューなども収録し，マス・フォア・インダストリ
の展開に資するものとする．

令和４年 10 月
マス・フォア・インダストリ研究所
所長 梶原

数学と物理における Julia の活用
マス・フォア・インダストリ研究 No.27, IMI, 九州大学
ISSN 2188-286X
発行日

2024 年 1 月 18 日

編

集

富谷昭夫，横山俊一，永井佑紀，寺崎敏志

発

行

九州大学マス・フォア・インダストリ研究所
〒819-0395 福岡市西区元岡 744
九州大学数理・IMI 事務室
TEL 092-802-4402

FAX 092-802-4405

URL https://www.imi.kyushu-u.ac.jp/
印

刷

城島印刷株式会社
〒810-0012 福岡市中央区白金 2 丁目 9 番 6 号
TEL 092-531-7102

FAX 092-524-4411

健司

シリーズ既刊
Issue

Author ／ Editor

Title

Published

マス・フォア・インダストリ
研究 No.1

穴田 啓晃
安田 貴徳
Xavier Dahan
櫻井 幸一

Functional Encryption as a Social Infrastructure and
Its Realization by Elliptic Curves and Lattices

26 February 2015

マス・フォア・インダストリ
研究 No.2

滝口
藤原

Collaboration Between Theory and Practice in
Inverse Problems

12 March 2015

マス・フォア・インダストリ
研究 No.3

筧

非線形数理モデルの諸相：連続，離散，超離散，
その先
Various aspects of nonlinear mathematical models
: continuous, discrete, ultra-discrete, and beyond

24 March 2015

マス・フォア・インダストリ
研究 No.4

穴田
安田
櫻井
寺西

啓晃
貴徳
幸一
勇

Next-generation Cryptography for
Privacy Protection and Decentralized Control and
Mathematical Structures to Support Techniques

29 January 2016

マス・フォア・インダストリ
研究 No.5

藤原
滝口

宏志
孝志

Mathematical Backgrounds and Future Progress of
Practical Inverse Problems

1 March 2016

マス・フォア・インダストリ
研究 No.6

松谷
佐伯
中川
上坂
濵田

茂樹
修
淳一
正晃
裕康

結晶のらせん転位の数理

10 January 2017

マス・フォア・インダストリ
研究 No.7

滝口
藤原

孝志
宏志

Collaboration among mathematics, engineering and
industry on various problems in infrastructure and
environment

1 March 2017

マス・フォア・インダストリ
研究 No.8

藤原
滝口

宏志
孝志

Practical inverse problems based on interdisciplinary
and industry-academia collaboration

20 February 2018

マス・フォア・インダストリ
研究 No.9

阿部
高島
縫田
安田

拓郎
克幸
光司
雅哉

代数的手法による数理暗号解析
Workshop on analysis of mathematical cryptography
via algebraic methods

1 March 2018

マス・フォア・インダストリ
研究 No.10

阿部
落合
高島
縫田
安田

拓郎
啓之
克幸
光司
雅哉

量子情報社会に向けた数理的アプローチ
Mathematical approach for quantum information
society

26 December 2018

孝志
宏志

三郎

（

）

Issue

Author ／ Editor

Title

Published

マス・フォア・インダストリ
研究 No.11

松谷
佐伯
中川
濵田
上坂

茂樹
修
淳一
裕康
正晃

結晶転位の先進数理解析
Advanced Mathematical Investigation for Dislocations

7 January 2019

マス・フォア・インダストリ
研究 No.12

滝口

孝志

Non-destructive inspection for concrete structures
and related topics

13 February 2019

マス・フォア・インダストリ
研究 No.13

宇波
長野
吉岡
田上
白井

耕一
智絵
秀和
大助
朋之

数理農学における時系列データのモデル化と解析
Modeling and Analysis of Time Series Data in MathAgro Sciences

28 February 2019

マス・フォア・インダストリ
研究 No.14

佐久間 弘文
大津 元一
小嶋 泉
福本 康秀
山本 昌宏
納谷 昌之

ドレスト光子に関する基礎的数理研究

18 March 2019

マス・フォア・インダストリ
研究 No.15

松谷
佐伯
中川
濵田
富安

結晶の界面，転位，構造の先進数理解析

2 December 2019

マス・フォア・インダストリ
研究 No.16

Takuro Abe
Yasuhiko Ikematsu
Koji Nuida
Yutaka Shikano
Katsuyuki Takashima
Masaya Yasuda

Quantum computation, post-quantum cryptography
and quantum codes

17 January 2020

マス・フォア・インダストリ
研究 No.17

河村
津曲
西澤
溝口

彰星
紀宏
弘毅
佳寛

代数・論理・幾何と情報科学―理論から実世界への
展開

10 February 2020

マス・フォア・インダストリ
研究 No.18

Takashi Takiguchi

New technologies for non-destructive and noninvasive inspections and their applications

21 February 2020

マス・フォア・インダストリ
研究 No.19

Hirofumi Sakuma
Motoichi Ohtsu
Masayuki Naya
Izumi Ojima
Yasuhide Fukumoto

Basic mathematical studies on dressed photon
phenomena

19 March 2020

茂樹
修
淳一
裕康
亮子

マス・フォア・インダストリ
研究 No.20

松谷 茂樹
井上 和俊
加葉田雄太朗
佐伯
修
垂水 竜一
内藤 久資
中川 淳一
濵田 裕康

材料科学における幾何と代数 I

24 November 2020

マス・フォア・インダストリ
研究 No.21

Daisuke Sakurai
Shigeo Takahashi
Naoki Hamada
Osamu Saeki
Hamish Carr
Takahiro Yamamoto

Fiber Topology Meets Applications

10 March 2021

マス・フォア・インダストリ
研究 No.22

巴山
中丸
大垣

竜来
啓
真二

数理的生成手法による CG とデジタルファブリケー
ション

15 March 2021

松谷
井上
落合
佐伯
垂水
内藤
中川
濵田
松江

茂樹
和俊
啓之
修
竜一
久資
淳一
裕康
要

材料科学における幾何と代数 II

11 November 2021

櫻井 大督
佐伯
修
高橋 成雄
Hamish Carr
山本 卓宏
濱田 直希

Fiber Topology Meets Applications 2

15 March 2022

佐久間弘文
大津 元一
小嶋
泉
福本 康秀

解析から設計に向けたオフシェル数理科学

28 March 2022

松谷 茂樹
緒方 勇太
落合 啓之
加葉田雄太朗
佐伯
修
濵田 裕康
松江
要

材料科学における幾何と代数 IV

28 November 2023

マス・フォア・インダストリ
研究 No.23

マス・フォア・インダストリ
研究 No.24

マス・フォア・インダストリ
研究 No.25

マス・フォア・インダストリ
研究 No.26

Institute of Mathematics for Industry
Kyushu University

九州大学マス・フォア・インダストリ研究所
〒819 - 0395 福岡市西区元岡744

https://www.imi.kyushu-u.ac.jp