線形代数の基本 の

線形代数の基本について

「線形代数」とは何でしょうか？

Wikipediaや教科書、3Blue1Brownのような優れたサイトには、技術的な定義がたくさんあります。ただ、ここではもっと気軽に考えてみましょう。

Julia（とPython）のメンテナーはこういう話が大好きです。 ただ、Exercismの運営陣は_それほどでもありません_。

このコンセプトでは、広大なテーマのうち、よりやさしい側面だけを扱います。ただし、先に警告しておきます。これは必然的に、かなり数学的なコンセプトです。

ベクトルとは？

誰に聞くか、また、かなり抽象的なものをどう図で表したいかによって、答えは変わります。

• Juliaでは、Vector{T}はArray{T, 1}の別名にすぎず、eltypeのTにはどんな型でも入ります。
• 物理学では（多くの場合）、ベクトルはN次元空間におけるlengthとdirectionを持つ矢印ですが、位置は固定されていません。
• 線形代数では、物理学者の言う矢印は空間に固定され、その尾は原点にあります。これで矢印の軸は不要になるので、ベクトルはN次元空間のpointとして表せます。N個の要素は、N本の軸それぞれに沿った原点からの距離を表します。

ここでは、文字列や文字のベクトルは無視しましょう。 このコンセプトでは数値型を扱います。Int、Float、Complexです。

行列とは？

ここでも、さまざまな答えがあります。

• 数値を長方形状に並べた2次元配列です（端がぎざぎざしていないもの）。正方行列はよくある特殊ケースです。
• 列ベクトルを横に並べたlinear combinationです。
• functionを他の入力型に適用するのと同じように、ベクトルに適用できるtransformationです。

正方行列の中には、特別な名前が付くほどよく登場するものがあります。

対角行列： 非ゼロの要素はすべてmain diagonal（左上から右下）に並びます。

julia> [1 0 0; 0 2 0; 0 0 3]
3×3 Matrix{Int64}:
 1  0  0
 0  2  0
 0  0  3

# a convenient shortcut
julia> diagm(1:3)
3×3 Matrix{Int64}:
 1  0  0
 0  2  0
 0  0  3

単位行列： 対角成分がすべて1の対角行列です（その理由は後ほど明らかになります）。 しばしばIと略されます。

julia> Matrix{Float64}(I, 3, 3)
3×3 Matrix{Float64}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

上三角行列： 対角成分とその上側に非ゼロの値があり、下側は0です。 下三角行列は想像がつくでしょう。

行列の転置

行と列を本当に入れ替えたいなら、permutedims()という関数がそれを行い、新しい行列を返します。

これはコピーを伴うので、大きな行列を扱うと遅く、メモリも大量に消費します。

線形代数の用途では、transpose()関数のほうが便利です。元の行列を包む遅延ラッパーをすばやく作ってくれます。

さらに便利なのがadjoint()関数で、複素数の虚部の符号も反転させます（なぜそれがそんなに重要かと不思議に思うなら、手っ取り早い答えの一つは量子力学です）。 この操作はよく使われるので、変数名にアポストロフィ'を付けるだけで随伴を作れます。

julia> m
2×3 Matrix{Int64}:
 1  2  3
 4  5  6

# new, full copy
julia> permutedims(m)
3×2 Matrix{Int64}:
 1  4
 2  5
 3  6

# lazy version
julia> transpose(m)
3×2 transpose(::Matrix{Int64}) with eltype Int64:
 1  4
 2  5
 3  6

julia> mc = [1+2im 2+3im; 3+2im 1+2im]
2×2 Matrix{Complex{Int64}}:
 1+2im  2+3im
 3+2im  1+2im

# lazy conjugate transpose
julia> adjoint(mc)
2×2 adjoint(::Matrix{Complex{Int64}}) with eltype Complex{Int64}:
 1-2im  3-2im
 2-3im  1-2im

# syntactic sugar for adjoint
julia> mc'
2×2 adjoint(::Matrix{Complex{Int64}}) with eltype Complex{Int64}:
 1-2im  3-2im
 2-3im  1-2im

乗算

このドキュメントの残りは、LinearAlgebraモジュールの機能を使います。 using LinearAlgebraの行を書けば名前空間に読み込まれますが、例の中では毎回繰り返しません（見た目がごちゃごちゃしすぎるので）。

ベクトル：要素ごとの積

これはベクトルの演算のコンセプトで説明しました。 演算子は.*で、入力のベクトル同士を要素ごとに掛け合わせ、入力と同じサイズ・型の出力を返します。

julia> [1, 2] .* [3, 4]
2-element Vector{Int64}:
 3
 8

ベクトル：内積

この非常によく使われる操作は、教科書ではu ⋅ vと書かれ、内積（ドット積）と呼ばれます。

内積は、要素ごとの積の_合計_に等しいです。

2つのベクトルは、通常の*演算子でも掛け算できますが、左側のベクトルがadjointである場合に限られます。これはu' * vと書くと便利です。 この点は、後ろの行列の乗算のセクションでよりはっきりするでしょう。

上付きの（_中央_の）ドットは、dot()関数の糖衣構文としてJuliaで使えます（\cdotと入力してtabキーを押します）。 随伴を指定する必要はありません。この部分は自動的に処理されるからです。

julia> using LinearAlgebra

# with element-wise syntax
julia> sum([1, 2] .* [3, 4])
11

# with u' * v syntax
julia> [1, 2]' * [3, 4]
11

# with dot()
julia> dot([1, 2], [3, 4])
11

# with \cdot syntax
julia> [1, 2] ⋅ [3, 4]
11

複素数値のベクトルでは、左側のベクトルが共役（虚部の符号を反転したもの）である必要があります。 dot関数はこれを自動で行い、u'の構文は明示的に行いますが、uとvが複素数だとsum(u .* v)は失敗します。

ベクトル：ベクトル積

昔、電磁気学の講義を受けたことのある人には、ちょっとした懐かしさがあるのではないでしょうか。トルクや角運動量のベクトルを計算したことのあるエンジニアにも。

内積が2つのベクトルをscalarに変換するのに対し、ベクトル積は2つの3次元ベクトルを3つ目の3次元ベクトルに変換します。

ベクトル空間の幾何学的な表現では、新しいベクトルは、2つの入力ベクトルを含む平面に_垂直_です。 2つの入力が平行な場合（符号を除いて）、それらは平面を定められないので、出力は[0, 0, 0]になります。

順序が重要です。u × v == -(v × u)です。 これが有名な右手の法則で、私たちの多くは、親指と2本の指を空間でひねりながらじっと見つめてきました（たいていは困惑した表情で）。

# with \times syntax
julia> [1, 2, 3] × [3, 4, 5]
3-element Vector{Int64}:
 -2
  4
 -2
# with cross()
julia> cross([1, 2, 3], [3, 4, 5])
3-element Vector{Int64}:
 -2
  4
 -2

この操作が長さ3のベクトルに限られていることに注目してください（直交するx, y, z軸を持つユークリッド空間と同じです）。 数学者は、他の次元にも喜んでこの操作を定義し、結果として理解しがたい言葉の羅列（高階テンソル、ウェッジ積、外積、マルチベクトルなど）を生み出します。 私たちの多くは、その時点で逃げ出すか、話題をさっと変えてしまいます！

ノルム

ベクトルはどれくらい「大きい」のでしょうか？

normは、ベクトルを適切なスカラーに還元することで、これを捉えようとするものです。

ノルムにはたくさんの種類がありますが、圧倒的によく使われるのは2ノルムで、√(v ⋅ v)です。

この二乗平均平方根の操作は、原点からの_ピタゴラス距離_です（N次元空間において）。 ベクトルを、尾を原点に置いた矢印として思い描くと、2ノルムは矢印の_長さ_です。

任意のpノルムは、pを2番目の引数として渡せば計算できます。 1ノルムはときに便利です。これは単に各要素の絶対値の合計です（つまり、計算がとても速く簡単です）。

# defaults to the 2-norm
julia> norm([1, 2, 3])
3.7416573867739413
# the 1-norm
julia> norm([1, -2, 3], 1)
6.0

行列の乗算

世界中の応用数学者が、この80年の大半を費やしてあらゆる計算を行列の乗算の繰り返しに変えてきたかのように思えることがあります。

これは、コンピューターがとても得意な操作です。

• 非常に繰り返しが多い。
• 効率よく並列化できる。
• より速くするために、多くの専用ハードウェアが開発されてきました。1970年代の8000万ドルもするCray-1から、おそらくノートPCに組み込まれているGPUまで、さまざまです。

細かい仕組みはかなり単純ですが、最初はあまり直感的ではありません。

行列Aにベクトルvを掛けて出力wを得る場合を考えてみましょう（線形代数の慣習で、行列には大文字、ベクトルには小文字を使います）。

Aの一番上の行とvの内積を取るとwの最初の要素になり、2番目の行からは2番目の要素が得られ、以下同様に続きます。

julia> A = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4
julia> v = [5, 6]
2-element Vector{Int64}:
 5
 6
julia> A * v
2-element Vector{Int64}:
 17  # equals [1, 2] ⋅ [5, 6]
 39  # equals [3, 4] ⋅ [5, 6]

行列 * 行列の乗算は、これを右側の行列の各列にわたって広げたものです。

C = A * Bでは、AがBの各列を掛けて、対応するCの列を作ると考えられます。つまり、行列とベクトルの乗算の繰り返しです。

同じことですが、Aの最初の行とBの各列の内積を取るとCの一番上の行になり、2番目の行からは2番目の行が得られ、以下同様に続く、とも言えます。

こうした頭の中でのイメージの持ち方にはいくつかあり、興味のある人はそれらを論じたMITの講義をまるごと見ることもできます。

julia> A
2×2 Matrix{Int64}:
 1  2
 3  4

julia> B = [5 7; 6 8]
2×2 Matrix{Int64}:
 5  7
 6  8

# left column is the same as A*v previously
julia> A * B
2×2 Matrix{Int64}:
 17  23
 39  53

julia> B * A
2×2 Matrix{Int64}:
 26  38
 30  44

上の例でわかるように、行列の乗算は_交換法則が成り立ちません_。A*BとB*Aの間には単純な関係がありません。

行列の乗算は、初めて触れる人にとって、文章を読んだだけでは図で思い描くのが難しいでしょう。 YouTubeには、それを図で示した動画がたくさんあります。「matrix multiplication」で検索して、好みのスタイル・詳しさ・言語のものを選んでください。

次元

2つのベクトルの内積は、それらが同じ長さであることに依存します。

同様に、行列の乗算では、左の行列の_列の数_が、右の行列の_行の数_と一致している必要があります。

サイズを(nrows, ncols)のタプルとして表すと、size(A)の出力のように、(a, b) * (b, c) -> (a, c)となります。 「内側の」次元、ここではbとbは、内積を取るのに適しています。 「外側の」次元、ここではaとcが、出力の次元を決めます。

長方形の行列の例です。

julia> D = reshape(1:6, 2, 3)
2×3 reshape(::UnitRange{Int64}, 2, 3) with eltype Int64:
 1  3  5
 2  4  6

julia> E = reshape(1:12, 3, 4)
3×4 reshape(::UnitRange{Int64}, 3, 4) with eltype Int64:
 1  4  7  10
 2  5  8  11
 3  6  9  12

julia> D * E
2×4 Matrix{Int64}:
 22  49   76  103
 28  64  100  136

# (3, 4) * (2, 3) not possible
julia> E * D
ERROR: DimensionMismatch: matrix A has axes (Base.OneTo(3),Base.OneTo(4)), matrix B has axes (Base.OneTo(2),Base.OneTo(3))

ベクトルの組を行列の流儀で掛ける方法には、2通りあります。

慣習的には、内積と同じ意味でu' * vを使います。 次元は(1, 3) * (3, 1)で、Juliaは(1, 1)の出力をスカラーに簡約します（例えばRとは異なります）。

あるいは、u * v'を使うと、次元は(3, 1) * (1, 3) -> (3, 3)で、すべての要素の組を掛け合わせ、ベクトルを行列に広げられます。

julia> u = [1, 2]
2-element Vector{Int64}:
 1
 2
julia> v = [3, 4]
2-element Vector{Int64}:
 3
 4
julia> u' * v
11
julia> u * v'
2×2 Matrix{Int64}:
 3  4
 6  8

u' * vは、_内積_と呼ばれることもあります。

（あまり使われない）u * v'は、対応して_外積_です。 これはテンソル積と関係がありますが、それはここでの範囲をはるかに超えています。

回転

特によくある行列の乗算の一種に、回転行列を使うものがあります。

2次元では、ベクトルを反時計回りにθラジアン回転させる、比較的単純な行列があります。

julia> rot2d(θ, vec) = [cos(θ) -sin(θ); sin(θ) cos(θ)] * vec
rot2d (generic function with 1 method)

# unit vector in the x direction
julia> i_hat = [1, 0]
2-element Vector{Int64}:
 1
 0

# rotate 45 degrees
julia> rot2d(π/4, i_hat)
2-element Vector{Float64}:
 0.7071067811865476
 0.7071067811865475

# rotate 90 degrees -> unit vector in the y direction, j_hat
julia> rot2d(π/2, i_hat)
2-element Vector{Float64}:
 6.123233995736766e-17  # zero, within numerical error
 1.0

3次元の回転には、2つの角度を使う、より複雑な行列が必要です。 これは、大量のLaTeXを埋め込まない限りMarkdownのドキュメントでは示しにくいので、式についてはWikipediaを確認してください。

OpenGLとその後継のような3次元グラフィックスでは、同次座標を使うのが慣習です。

各頂点は4次元ベクトル（言い換えれば行列の列）で表されます。(x, y, z)にある点なら[x, y, z, 1.0]です。

これにより、より凝った変換行列が作れます。 回転はやはりA[1:3, 1:3]にあり、平行移動はA[1:3, 4]の[Δx, Δy, Δz]、拡大縮小は対角成分にあり、ほかにもせん断や透視などの可能性があります。

パフォーマンス

世界中の計算機科学の学科には、行列の乗算アルゴリズムの小さな漸進的な改善について、1つ以上の章を割いた博士論文がたくさんあります。 これは_本当に_重要で、さまざまな組織が自らの利益のためにその研究に資金を出すことをいとわずにいます。

パフォーマンスは大きなテーマで、深く立ち入れませんが、説明のために、いろいろなサイズのランダムな行列を掛け合わせてみましょう。

以下のコードは、手早くざっくりした見積もりです（よりよい方法にはBenchmarkTools.jlを使ってください）。

使ったシステムは、小型の430ドル（米国）のPCです。Ryzen 9プロセッサ、32GBのRAM、Linux Mint 22.1、Julia 1.11.6です。

julia> mmul(A) = A * A
mmul (generic function with 1 method)
julia> A = rand(Float64, 1_000, 1_000);
julia> @time mmul(A);
  0.011055 seconds (3 allocations: 7.629 MiB)
julia> A = rand(Float64, 10_000, 10_000);
julia> @time mmul(A);
  7.236284 seconds (1.61 k allocations: 763.027 MiB, 17.92% gc time, 0.15% compilation time)

おおよそ、100万要素の行列の組は数ミリ秒、1億要素の行列は数秒かかりました。 さらに何桁も大きくするには、よりよいハードウェアが必要でしょう...