本記事はJulia Advent Calendar 2022の12/23の記事です。

東京大学で働いている松井と申します。

線形代数の講義における演習（実際にコードを書き行列演算を行う）の重要性を感じています。 そのためにjuliaを使えないかと思い至り、pythonとの比較に焦点を当て思っていることを述べます。

線形代数における演習の意義

線形代数は工学全般において重要で基盤的な学問体系ですが、なかなかとっつきにくいものです。その理由の一つは線形代数の諸アルゴリズムは最終的には計算機で実行するにも関わらず、学生は自分の手を動かしてコーディングする機会が少ない点だと感じます。多くの大学のカリキュラムでは大学初年次に線形代数講義があると思いますが、座学がメインであることが多いと思います。本当は、座学と並行して実際にコーディングして行列演算を行う「演習講義」があれば、理解が深まるだろうと感じます。たとえば掃き出し法一つをとっても、実際に自分でコーディングして挙動を確かめることで、本当に行列が変形していくんだな、ということがわかります。

一方で、大学１年や２年の段階で線形代数演習講義を行うことはなかなか難しいとも感じます。すなわち、

  • 線形代数を学ぶころ、学生はプログラミングをそもそも習っていないか習っている途中
  • 仮にプログラミングが可能だとしても、大学１、２年生にc/c++で数値計算はキツい。pythonが現実的な選択肢になる（他カリキュラムとの兼ね合い的にもpythonは有利）
  • しかし後述するように、python（numpy）で行列演算は「クセ」が強い

そこで、pythonのようにクセが無く、かつ初学者にも扱いやすい言語として、juliaに注目しました。juliaは様々な点でpythonに比べて「素直で楽」な面があります。それを以下に紹介していきます。

juliaの簡単実行

まず、juliaを簡単に実行する方法を紹介します。最近GitHub Codespacesが毎月60時間までなら誰でも利用可能になりました。なので、GitHub Codespacesを使ってブラウザベースで環境を作ると簡単です。

まず作業用のリポジトリをGitHub上に作ります。そして以下のようなjsonファイルをリポジトリ直下に.devcontainer.jsonという名前で配置します。あとはgithubの画面からcodespace（vscodeライクなウェブエディタ）を起動するだけで、juliaを利用可能なcodespaceが立ち上がります。

// See https://github.com/julia-vscode/julia-devcontainer/blob/master/Dockerfile for image contents
{
	"name": "Julia (Community)",
	"image": "ghcr.io/julia-vscode/julia-devcontainer",

	// Configure tool-specific properties.
	"customizations": {
		// Configure properties specific to VS Code.
		"vscode": {
			// Add the IDs of extensions you want installed when the container is created.
			"extensions": [
				"julialang.language-julia"
			]
		}
	},

	"postCreateCommand": "/julia-devcontainer-scripts/postcreate.jl",

	"remoteUser": "vscode"
}

その例を以下のリポジトリ置きましたので参考にしてください

実は本当に単にjuliaを使いたいだけ（ソースを保存する必要もない）なのであれば、自分でリポジトリを作らずとも、上記リポジトリに対しcodespaceを作成すればもう使えちゃいます。いずれにせよ、codespaceの無料利用時間を消費しますので、使いすぎには注意してください。

codespacesが無い場合は、大学の計算機に入っていない環境を学生に課すのはかなり難しかったため、気軽にjuliaを講義で使うことは現実的ではありませんでした。ですが、codespaceのおかげで可能になりつつあります。

本章の記述のために、以下の記事の内容を大変参考にさせていただきました。

また、他の簡単実行環境として、以下の記事ではbinderを用いたワンクリック実行が紹介されています。

行列演算表記の比較：julia vs python

さて、それでは行列に対する基本演算について、juliaとpythonを比べていってみましょう。 ここで注目したいのは、プログラミングを習ったばかりの初学者が躓かないかどうかです。

まずはjuliaです。juliaでは、以下のように「行列」、「縦ベクトル」、「横ベクトル」を定義できます。

# 行列
A = [10 20 30
     40 50 60
     70 80 90]
# A = [10 20 30; 40 50 60; 70 80 90]でもいい

# 縦ベクトル
a = [1
     2
     3]
# a = [1; 2; 3]でもいい

# 横ベクトル
b = [4 5 6]

上記のように、自然に定義できます。そして、行列演算は以下のようになります。

# 3x1ベクトルと1x3ベクトルの行列積。すなわち外積。行列になる。
a * b
# 3×3 Matrix{Int64}:
#   4   5   6
#   8  10  12
#  12  15  18

# 1x3ベクトルと3x1ベクトルの行列積。すなわち内積。スカラになる。
b * a
# 1-element Vector{Int64}:
#  32

# 1x3と1x3の行列積や、3x1と3x1の行列積はエラーになる
a * a  # error
b * b  # error

# 3x3行列と3x1ベクトルの行列積。3x1になる。
A * a
# 3-element Vector{Int64}:
#  140
#  320
#  500

# 3x1と3x3だとエラー
a * A  # error

# 1x3ベクトルと3x3行列の行列積。1x3になる。
b * A
# 1×3 Matrix{Int64}:
#  660  810  960

# 3x3と1x3だとエラー
A * b  # error

上記のように、全ての行列演算は自明かつ自然に行うことが出来ました。juliaの場合は、行列積には通常のアスタリスクが演算子として割り振られています。

次に、pythonの場合を見てみましょう。

import numpy as np

# 行列
A = np.array([[10, 20, 30],
              [40, 50, 60],
              [70, 80, 90]])

# 縦ベクトル
a = np.array([[1],
              [2],
              [3]])
a.shape  # (3, 1)

# 横ベクトル
b = np.array([[4, 5, 6]])
b.shape  # (1, 3)

pythonの場合はnumpyを呼び、np.arrayでリストをラップする必要がありますが、とりあえず自然に定義できます。行列演算は以下になります。

# 3x1ベクトルと1x3ベクトルの行列積。すなわち外積。行列になる。
a @ b
# array([[ 4,  5,  6],
#        [ 8, 10, 12],
#        [12, 15, 18]])

# 1x3ベクトルと3x1ベクトルの行列積。すなわち内積。スカラになる。
# （正確にはshapeが(1, 1)の行列だが、とりあえず気にしない）
b @ a
# array([[32]])

# 1x3と1x3の行列積や、3x1と3x1の行列積はエラーになる
a @ a  # error
b @ b  # error

# 3x3行列と3x1ベクトルの行列積。3x1になる。
A @ a
# array([[140],
#        [320],
#        [500]])

# 3x1と3x3だとエラー
a @ A  # error

# 1x3ベクトルと3x3行列の行列積。1x3になる。
b @ A
# array([[660, 810, 960]])

# 3x3と1x3だとエラー
A @ b  # error

pythonでは行列演算には@の演算子を使います。上記の通り、juliaと同様に、全ての演算が自然に定義されました。なので、問題はないように思います。

pythonの難しい点

ところが、問題はここからです。pythonでは、一次元のリストからnp.arrayを作る場合、shapeが(3, 1)でも(3, 1)でもなく、(3, )という状態のベクトルになります。これは行列を考えなければ自然な定義なのですが、このベクトルの存在が話をややこしくします。

# 一次元のリストからは「ベクトル」が出来る
c = np.array([1, 2, 3])
c.shape  # (3, )になる！ (3, 1)でも(1, 3)でもない

# 自分自身と行列演算できてしまう！
c @ c
# 14

# 3x3行列に対して、右からも左からも行列積が計算できてしまう！
c @ A
# array([300, 360, 420])

A @ c
# array([140, 320, 500])

上記のように、ベクトルは「自分自身と行列積を計算できてしまう」「行列に対して右からも左からも行列積を計算できてしまう」という、行列演算的には非常に不思議な挙動をとります。そして、初学者にとっては、ベクトルを表現したいときに上記のように一次元リストから作ろうとするのは自然です。なので、非常にややこしい自体になってしまっています。

さらに話をややこしくするのは、pythonにはただのリストも存在するという点です。

# 一次元のリストそのもの
d = [1, 2, 3]

# リスト同士では@演算子は使えない
d @ d
# エラー

# 3x3行列に対して、右からも左からも行列積が計算できてしまう！
d @ A
# array([300, 360, 420])

A @ d
# array([140, 320, 500])

上記のように、ただのリスト同士だと@は計算できませんが、np.arrayであるAとの演算では利用できます。。。

上記をまとめると、pythonでは「数字３つが並んでいるベクトル」というものを表すというだけで、少なくとも下記のように４通り存在します。

v1 = np.array([[1], [2], [3]])  # 3x1の行列 shape=(3, 1)
v2 = np.array([[1, 2, 3]])      # 1x3の行列 shape=(1, 3)
v3 = np.array([1, 2, 3])        # ベクトル shape=(3, )
v4 = [1, 2, 3]                  # リスト

この状態を初学者に提示するのは、なかなかにハードです。本筋の行列演算に入る前に疲れ切ってしまう可能性があります。この点が、pythonを線形代数講義で使いづらい個人的な理由です。自分自身がpythonを勉強したときを思い返しても、numpyのarrayとpython標準のリストが二種類あることを理解するのは中々大変だった覚えがあります。

本章を執筆するにあたり以下の記事を大変参考にさせていただきました。より包括的で詳細な説明はこちらを参考にしてください。

juliaなら全部OK、というわけでもない

juliaなら全部OKかというと、そういうわけでもありません。実は既にここまでの議論でも、juliaは謎の挙動をとります。縦ベクトル、横ベクトルを作るとき、その型は実は次のようになります。

a = [1; 2; 3]
# 3-element Vector{Int64}:
#  1
#  2
#  3

b = [1 2 3]
# 1×3 Matrix{Int64}:
#  1  2  3

横ベクトルの型が1x3 Matrixであるのは自然ですが、縦の場合は3x1 Matrixではなく3-element Vectorになっています。すなわち、縦ベクトルと横ベクトルで型が違ってしまいます。 実際に計算するうえではこの型の違いは問題ないのですが、それにしろこれは中々に気持ち悪いです。

juliaの思想的には、ユーザはこの部分を気にせずとも多重ディスパッチにより高速な計算が可能だし、可能であるべき、というものがあるように感じます。なので、深入りすると気持ち悪いですが、とりあえず何も考えずに使う場合は気にしないでいいということになるかと思います。

ここまでをまとめると、「juliaもpythonも気持ち悪い部分があるが、線形代数を初学者が学ぶという面では、最初にとっつきやすいのはjuliaであろう」ということになるかと思います。

このあたりは、以下の記事に詳細にまとめられています。

番外編：pythonが支配的すぎる弊害

話は逸れますが、私の専門のcomputer vision分野ではあまりにもpythonがデファクトになりすぎており、色々と弊害も発生しています。 一番の問題は、python記述と数式がゴッチャになってしまい、numpyのブロードキャストを数式に書いてくる人が本当に多いということです。 numpyでは、次のように、スカラを行列に足すとブロードキャストされ、行列の全ての要素に対しスカラを足してくれます。

import numpy as np
a = 10
Y = np.array([[1, 2], [3, 4]])
Y.shape  # (2, 2)
X = a + Y
# array([[11, 12],
#        [13, 14]])

これはOKなのですが、このブロードキャストはあくまでnumpyの機能に過ぎません。ところが、これを式にも書いてくる人が本当に多いです。

\[X = a + Y\]

上記のように書くことはできません。展開すれば、明らかにおかしいことがわかります。

\[10 + \begin{bmatrix} 1 & 2 \\ 3 & 4 \\ \end{bmatrix}\]

このような演算は出来ませんよね？ですが、私がCVPRというCV分野トップの会議論文を5本査読する場合、そのうち13本は上記のようにスカラを行列に足してきます。キビシー！

このように間違った数式を書いてしまうと、その後の全ての数式も意味がない間違ったものになってしまいますので、絶対にこういう表記はしないようにしましょう。こちらも参考にしてください。

このようについついpythonの表記を数式に書いてしまうことを防ぐためには、python以外の言語も普段から使うようにするのがオススメです。juliaの場合は自動でブロードキャストはしてくれません。要素を分配するドットをつける必要があります。

a = 10
Y = [1 2; 3 4]

X = a + Y  # エラー。自動でブロードキャストしない

X = a .+ Y # ドットをつけると分配してくれる
# 2×2 Matrix{Int64}:
#  11  12
#  13  14

これを知っておけば、numpyのブロードキャストはあくまでpython独自の表記にすぎないことがわかります。

pythonは2022現在あまりに便利で支配的な言語になっているので、だからこそ逆に、それ以外の言語を勉強する価値も高いと思います。

（余談ですが、juliaはMATLABから強い影響を受けている言語です。pythonが流行る前はCV分野の人は全員MATLABを使っていました。juliaという新しい言語を使うことはMATLABに回帰する風味もあり、中々に面白いです）

juliaが用いられている講義

最後に、自分が観測した範囲で、juliaが講義で実際に使われている事例を紹介します。他にもありましたら是非教えてください。追記いたします。

国内では、大阪大学の数値計算系の講義でjuliaは使われているようです。以下の講義では積極的に講義資料を公開されており、大変参考になります。

【2023.01.04追記】都立大の菅原先生は「大学新入生にJuliaで演習を行う」と言う、本記事の狙いと大いに重なる講義をされています。（Julia公式ページには教育現場での使用例として都立大（首都大）ロゴがあります！）

また、juliaの創始者の１人であるAlan Edelman先生はMITの教授である（julialab at MIT）とのこともあり、juliaはMITの講義での使用が盛んなようです。以下にいくつかの例をあげます。

また、東京大学ではpythonのnotebookに対する自動採点システムである「PLAGS UT」が、佐藤重幸先生らにより整備・運用されています（「Pythonプログラミング入門 #utpython」講義HPのトップのイントロダクションスライドの中を参照）。PLAGS UTはnotebookに対する自動採点システムであり、本質的にはpythonでなくとも実行可能であるため、juliaに対しても実行可能であると考えます。このへんを仕込んでいけると面白いと考えます（興味ある先生がいらっしゃいましたらご連絡いただけますと。。。）

まとめ

線形代数演習講義にjuliaを導入できないかを考えました。そのために、juliaの簡単実行方法の紹介、pythonとの比較を行い、最後にjuliaを使っている講義を紹介しました。juliaのような新しい言語を考えることは、普段pythonに凝り固まっていしまっている頭に別の視点を導入することになり、とても新鮮な気持ちになりますね。