📏

Juliaの行列・ベクトルを完全に理解すっぞ！！

2020/12/22に公開

Julia Advent Calendar 2020/12/23

小ネタと思ってましたが、調べたら意外と非自明でした。なので中ネタくらいです。
Juliaでの行列やベクトルの扱いについて書いていきます。

事の発端

@ceptreeさん・@taketo1024さんとのリプライで、Juliaでの行列やベクトルについて話して、Juliaの行列とベクトルの詳細は紛らわしい部分があるなーと思った次第です。
なので色々調べてまとめました！

行列とベクトルの簡単な導入

Juliaでの行列やベクトルの扱いについて、よくある説明を復習しようと思います。

Juliaでの1次元配列は縦ベクトル

Juliaでは、1次元配列は縦ベクトルとして扱われます。

v = [1,2,3]

と書いても

v = [1
     2
     3]

と書いても

v = [1;2;3]

と書いても同じものとして解釈されます。
なぜこのような記法がいくつもあるかと言うと...

[1,2,3]は普通の配列の書き方で、"普通の配列"は縦ベクトルと解釈して欲しい
[1(改行)2(改行)3]は縦ベクトルに見えるので、これも縦ベクトルとして解釈して欲しい
[1;2;3]において、;は改行の代わりの記号として解釈して欲しい

といった背景からだと思います。

ちなみに、[1,2,3]の型はArray{Int64,1}で、「要素が64ビット整数Int64の1次元配列」という意味になります。Vector{Int64}と書いても同じ型を意味します。

他にも、等間隔に並んだ実数1:2:7などもベクトルとして解釈されます。(初項1, 間隔2, 末項7の有限列)
1:2:7の型はStepRange{Int64,Int64}で、AbstractVector{Int}の部分型(subtype)です。

Juliaでの2次元配列は行列

Juliaでは、2次元配列は行列として扱われます。

M = [1 2 3;4 5 6]

と書いても

M = [1 2 3
     4 5 6]

と書いても同じです。
複数の表記法のある理由はベクトルのときと同様ですね。

またJuliaでは、インデックスの向きがrow-major orderではなく、column-major orderです。
以下のコードは1から6までの列を2×3にreshapeしているのですが、その結果は[1 2 3;4 5 6]とは異なっています。

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

転置・随伴

転置を計算するにはtransposeを使います。

julia> M
2×3 Array{Int64,2}:
 1  2  3
 4  5  6

julia> transpose(M)
3×2 LinearAlgebra.Transpose{Int64,Array{Int64,2}}:
 1  4
 2  5
 3  6

随伴行列を計算するにはadjointを使います。

julia> adjoint(M)
3×2 LinearAlgebra.Adjoint{Int64,Array{Int64,2}}:
 1  4
 2  5
 3  6

adjoint(M)の変わりにM'とも書けます。

julia> M'
3×2 LinearAlgebra.Adjoint{Int64,Array{Int64,2}}:
 1  4
 2  5
 3  6

随伴とは、転置+複素共役でした。
以下のコードで型までちゃんと合ってることが確認できます。

julia> conj(transpose(M))
3×2 LinearAlgebra.Adjoint{Int64,Array{Int64,2}}:
 1  4
 2  5
 3  6

非自明なところ

ここからが本題です！！

横ベクトルや3×1行列は直接宣言できない

コード例

冒頭で述べたように、Juliaでは通常のベクトルは縦ベクトルです。
では、横ベクトルを直接に定義する方法はあるのでしょうか？
行列が

[1 2 3;4 5 6]

で定義できたのですから、[1 2 3]は横ベクトルになる気がしますね。

julia> [1 2 3]
1×3 Array{Int64,2}:
 1  2  3

しかし結果は1×3行列(Array{Int64,2})でした。

では[1;2]を実行すれば2×1行列になるでしょうか？
いいえ、冒頭で見たようにこれは2要素のベクトルArray{Int64,1}(==Vector{Int64})になります。

Julia1.7からは直接に2×1行列を宣言できるようになりました！

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

セミコロン;を並べた数に応じて、配列のどのインデックスで改行するかを決めることができます。

julia> [1;2;;3;4;;5;6;;;2;2;;2;2;;2;2;;;3;3;;3;3;;3;3;;;4;4;;4;4;;4;4]
2×3×4 Array{Int64, 3}:
[:, :, 1] =
 1  3  5
 2  4  6

[:, :, 2] =
 2  2  2
 2  2  2

[:, :, 3] =
 3  3  3
 3  3  3

[:, :, 4] =
 4  4  4
 4  4  4

julia> [3;;;;]
1×1×1×1 Array{Int64, 4}:
[:, :, 1, 1] =
 3

なぜこうなっているのか？

横ベクトルを直接宣言できない理由

行儀の悪いユーザーが、横ベクトルを基本としたコードを書くのを減らすためでしょう。
実用上は特にこれで困ることは無いはずです。

線形変換は

\begin{pmatrix} x' & y' \end{pmatrix} =\begin{pmatrix} x & y \end{pmatrix} \begin{pmatrix} a & b \\ c & d \end{pmatrix}

よりも

\begin{pmatrix} x' \\ y' \end{pmatrix} =\begin{pmatrix} a & b \\ c & d \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix}

の方が気持ち良いですよね。
なので縦ベクトルを推します。

3×1行列を直接宣言できない理由

「3×1行列」と「縦ベクトル」を混同させないためでしょう。
実用上は、特にこれで困ることは無いでしょう。
どうしても定義したい場合は、以下のようにreshapeを使えばOKです。

julia> reshape([1,2,3],3,1)
3×1 Array{Int64,2}:
 1
 2
 3

2つのベクトルの積

コード例

v = [1,2,3]
w = [2,3,4]

として2つのベクトルを定義したとしましょう。
ここでこれらのベクトル同士の積v*wを計算したいと思います。
Juliaでは以下のうち、どの積として解釈されるでしょうか？

内積(ドット積)
外積(クロス積)
アダマール積(要素ごとの積)

正解は..v*wはエラーになります！

julia> v*w
ERROR: MethodError: no method matching *(::Array{Int64,1}, ::Array{Int64,1})
Closest candidates are:
  *(::Any, ::Any, ::Any, ::Any...) at operators.jl:538
  *(::LinearAlgebra.Adjoint{var"#s828",var"#s8281"} where var"#s8281"<:(AbstractArray{T,1} where T) where var"#s828"<:Number, ::AbstractArray{var"#s827",1} where var"#s827"<:Number) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/adjtrans.jl:283
  *(::LinearAlgebra.Transpose{T,var"#s828"} where var"#s828"<:(AbstractArray{T,1} where T), ::AbstractArray{T,1}) where T<:Real at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/adjtrans.jl:284
  ...
Stacktrace:
 [1] top-level scope at REPL[12]:1

では上記3つの積をどうやって計算するのかと言うと、以下のようにすればOKです。

内積(ドット積)
- LinearAlgebra.dot(v,w)
外積(クロス積)
- LinearAlgebra.cross(v,w)
アダマール積(要素ごとの積)
- v.*w

内積や外積の計算のためには

using LinearAlgebra

が必要で、アダマール積を計算するには単にbroadcastをすれば良いだけでした。

ちなみに、内積を計算したいだけならLinearAlgebraを使わずに、v'*wとしても計算できます。

なぜこうなっているのか？

ベクトル同士のナマの掛け算が定義されない理由

理由は、Juliaでは「記号*は標準的な積を表すから」になります。
「行列同士の積」や「実数とベクトルの積」などは標準的なものが決まっていますが、「ベクトル同士の積」には(少なくとも)上記の3つがあり、どれが標準的かは決まっていません。
そのためv*wがエラーになった訳です。

*が使えないと不便に感じる方もいるかと思いますが、ご安心ください！
JuliaではUnicode文字を使ったコーディングを推奨しているので(using LinearAlgebraした後で)以下のようにドット積とクロス積が計算できます！

julia> [1,2,3]⋅[4,8,-3]
11

julia> [1,2,3]×[4,8,-3]
3-element Array{Int64,1}:
 -30
  15
   0

これらの記号⋅、×はそれぞれ\cdot、\timesで出すことができます。

横ベクトルはベクトルではない

ちょっと何を言ってるかわかりませんが、こういうことです。

julia> [1,2,3]' isa AbstractVector
false

やっぱり、ちょっとよくわかりませんね。

もう少し確認してみましょう。

julia> [1,2,3] isa AbstractVector
true

julia> [1,2,3]' isa AbstractVector
false

julia> [1,2,3]'' isa AbstractVector
true

どうやらJuliaのベクトルは、随伴を取る度にAbstractVectorになったりならなかったりするようです。(transposeでも同様です)

julia> [1,2,3]''
3-element Array{Int64,1}:
 1
 2
 3

偶数回の随伴(転置)を取れば元に戻ってきてますね。

結局どうなってるかと言うと… 横ベクトルはベクトルではなく行列でした！

julia> [1,2,3] isa AbstractMatrix
false

julia> [1,2,3]' isa AbstractMatrix
true

なぜこうなっているのか？

横ベクトルが`AbstractVector`ではない理由

意外でしたが、JuliaにおいてはAbstractVector{T}どうしの足し算+は(要素数が合う限り)必ず定義されている必要があります。
例えば[1,2,3] + (3:2:7)は[4,7,10]として評価されます。
もし仮に、縦ベクトル[1,2,3]と横ベクトル[1,2,3]'の足し算が定義されていたとすれば、その結果は縦ベクトルでしょうか？横ベクトルでしょうか？
これが決められないから、そもそも「縦ベクトルと横ベクトルの足し算」を定義したくないという訳です。
そのためには、横ベクトルはAbstractVectorであってはなりません。
JuliaでのVector(AbstractVector)は1次元配列よりも縦ベクトルと考えた方が良いです。

少し省略していますが、行列やベクトルの型の階層関係は以下のようになっています。
このような階層関係はsubtypesやsupertypeで調べることができます。

横ベクトルが`AbstractRowVector`ではない理由

新たに抽象型AbstractRowVectorを用意して、以下のような階層関係にした方が良いと思う方もいるかも知れません。

しかし、これだと問題があります。

AbstractMatrixやAbstractVectorは単にAbstractArray{T,2} where TやAbstractArray{T,1} where Tのエイリアスになっています。
これによって型の包含関係AbstractMatrix <: AbstractArrayが従っているという訳です。
じゃあその場合、AbstractRowVectorをどこに入れば良いでしょうか？
少なくともAbstractRowVectorをAbstractArrayのsubtypeにするのは難しそうですね。
supertype(AbstractArray)はAnyなので、もはやAbstractRowVectorの入る場所は無さそうです。

詳しくはこちらのissueのコメントをご覧ください↓

横ベクトルが`AbstractMatrix`である理由

ここは少し汚いと思いますが、妥当な折衷案のようにも思えます。

既に見たように、AbstractVector(==AbstractArray{T,1} where T)は使えません。
既に見たように、AbstractRowVectorを定義するのはあまり良い方法では無さそうです。
ここでAbstractMatrixです！
- 「行数を1に制限した行列」のことを「横ベクトル」と考える訳です。
- ここでの横ベクトルとは、Adjoint{T,Vector{T}}やTranspose{T,Vector{T}}などのことです。
数学的には、これで良いんですか？
- そもそも、初等的な線形代数での横ベクトルの導入というのもかなり怪しいものです。
- 縦ベクトル全体の空間は $\mathbb{R}^3$ だとして、横ベクトル全体の空間は何でしょうか？ $\mathbb{R}^3$ であれば縦ベクトルと横ベクトルを区別できないということになります。
- 定義の面倒を避けるために、初等的には "数を横に並べたもの" として横ベクトルが導入されることが多いです。
- 数学の寛大さに倣って、「計算が同型なら横ベクトルと解釈してOK」と考えましょう。(まあ中身は行列なんですけれども)

添字の上下を考えない理由

でもやっぱり、「横ベクトルは行列」と言われても、かなり気持ち悪いですよね。
一番"正しい"方法は、「添字の上下や、成分の共変・反変を全部考えること」でしょう。
そうなると、

縦ベクトル((1,0)型テンソル)の成分 $v^i$
横ベクトル((0,1)型テンソル)の成分 $v_i$
(2,0)型テンソルの成分(行列) $A^{ij}$
(1,1)型テンソルの成分(行列) $A^{i}_{j}$
(0,2)型テンソルの成分(行列) $A_{ij}$
...

のように、全部を区別できて嬉しい！という訳です。

でも、これって本当に私たちが望んでいたことでしょうか？
$n$ 階テンソルには $n+1$ 種類があることになりますし、実装の度に毎回悩むのは大変です。
しかも、素朴な行列計算がしたい時に共変・反変とかは考えたくないですよね。
「通常のベクトルは縦ベクトル」という思想も崩れてきそうな雰囲気があります。

そう、私たちに必要な計算は、本来は素朴な行列計算だったはずです。
なので添字の上下に拘らずに、横ベクトルとしてのTransposeやAdjointを有り難く受け入れましょう。実用上は、中身がどうなっているかを気にしなくて良い場合が殆どなはずです。

どうしても添字の上下区別をしたい場合は、自分で型を定義して、計算させることも可能です。(やっていきましょう！💪)

番外: なんで横ベクトル？

「こんなに紛らわしいのに何故、転置とか随伴を考えてるんですか？」
「もう横ベクトルとか要らなくないですか？」

まあそうかも知れないですね。

でもベクトルv, wの内積を計算するときに

v'*w

って書けると便利じゃないですか？

行列に対する随伴 ' や転置transposeは必要だとして、これを使って「縦ベクトルの随伴」や「縦ベクトルの転置」も定義したくないですか？

私はあった方が便利だと思います。

collect∘adjoint、collect∘transposeの挙動

(∘は写像の合成を表す記号です。\circで出せます)
(縦)ベクトルの随伴(adjoint)を取って、Arrayに変換(collect)してみましょう。
つまり「横ベクトルを通常の配列に変換するとどうなるか」ということを考えます。

julia> (collect∘adjoint)([2,3])
1×2 Array{Int64,2}:
 2  3

横ベクトル(Adjoint{T,<:AbstractVector{T}} where T, Transpose{T,<:AbstractVector{T}} where T)が行列だというのは既に見たので当たり前の結果に思えますね。

では次の例はどうでしょうか？

julia> [1,2]'*[1,2]
5

julia> collect([1,2]')*[1,2]
1-element Array{Int64,1}:
 5

collectを使う前後で結果が変わりますね。
嫌な感じがしますが、何が理由なのでしょうか？

なぜこうなっているのか？

`collect`有りで1要素のベクトルが得られる理由

Juliaでは、「行列 × ベクトル」は「ベクトル」を返します。
これは「(1×n)行列 × (n要素)ベクトル」にも適用されるルールで、型安定性のために必須となります。
なのでcollect有りのケースで1要素のベクトルが得られたという訳です。

`collect`無しでスカラーが得られる理由

Adjoint{Int64,Array{Int64,1}}はAdjoint{T,<:AbstractVector{T}} where Tの部分型(subtype)なので、計算の実行前に「行数が1であること」が予め分かっています。
なので、Juliaは多重ディスパッチを使って「(1×n)行列 × (n要素)ベクトル」に対してスカラーを返すことができる訳です。

これによって内積v' * wの結果がちゃんとスカラーになっていた訳ですね。
「スカラーを返す親切さ」と「掛け算*の統一性」を天秤にかけた結果、前者が選ばれた訳です。
(個人的には「掛け算*の統一性」が失われたのは、少し気をつけるとカバーできる程度なので無問題だと思います。初見で驚きこそしましたが。)

Transpose、Adjointの型

コード例

随伴adjointを計算してみましょう

julia> [1 2;3 4]'
2×2 Adjoint{Int64,Array{Int64,2}}:
 1  3
 2  4

ここでのAdjoint{Int64,Array{Int64,2}}って何を表している型なのでしょうか？
Int64が2回入ってるって冗長じゃないですか？

なぜこうなっているのか？

型の意味

第一印象としてAdjoint{Int64,Array{Int64,2}}はややこしそうですが、そんなに難しくは無くて

Adjoint: 随伴を表す型
Int64: 要素の型(パラメトリック型に必要)
- Adjoint{T,S} <: AbstractMatrix{T}として定義されてます。
Array{Int64,2} 随伴を取られる側の型

というようになっています。

型パラメータが冗長な理由 (冗長に見えるだけの理由)

結論から言えば、上記の「Adjoint{T,S} <: AbstractMatrix{T}として定義されてます」が理由です。
でもやっぱり、Sに要素の型の情報が入ってるなら不要な気がしますよね。

しかし、本当に必要なのです！
具体的には以下のような状況が考えられます。

自作の型MyMatrix <: AbstractMatrixを定義したとします。
そしてMyMarixは具象型で、その「要素の型」の型推論が困難だとします。

この場合に、型パラメータに「MyMatrixの要素の型」の情報が入っていないとAdjointのための型推論が難しくなるのです。

実はJuliaでは、"型を実行"したときに、「そのような型が存在し得るか」は評価してくれません。
つまり例えば

julia> AbstractVector
AbstractArray{T,1} where T

のように型のエイリアスを調べることはできますが

julia> Adjoint{Int64, Vector{Float64}}
Adjoint{Int64,Array{Float64,1}}

のように、「存在しない型への警告」は一切出してくれません。

より極端な例を出すなら、以下もエラーになりません。

julia> Array{'a','b'}
Array{'a','b'}

エラーが欲しい気持ちもありますが、パラメトリック型の都合上、仕方ないと思われます。

関数transposeとコンストラクタTranspose

Juliaでは、関数とコンストラクタが同じように振る舞うことがあります。

julia> transpose([1 2;3 4])
2×2 Transpose{Int64,Array{Int64,2}}:
 1  3
 2  4

julia> Transpose([1 2;3 4])
2×2 Transpose{Int64,Array{Int64,2}}:
 1  3
 2  4

しかし、実はtransposeとTransposeは全く同じ振る舞いをする訳ではありません。
これは下の例から分かると思います。

julia> transpose(transpose([1 2;3 4]))
2×2 Array{Int64,2}:
 1  2
 3  4

julia> Transpose(Transpose([1 2;3 4]))
2×2 Transpose{Int64,Transpose{Int64,Array{Int64,2}}}:
 1  2
 3  4

なぜこうなっているのか？

関数とコンストラクタの違い

実装を確認してみましょう。
Juliaでの実装を確認するには@lessマクロが超便利です。

julia> @less transpose([1,2,3])

julia> @less transpose(transpose([1,2,3]))

結果を一部抜粋します↓

transpose(A::AbstractVecOrMat) = Transpose(A)

(中略)

transpose(A::Transpose) = A.parent

(中略)

struct Transpose{T,S} <: AbstractMatrix{T}
    parent::S
    function Transpose{T,S}(A::S) where {T,S}
        checkeltype_transpose(T, eltype(A))
        new(A)
    end
end

そうです！多重ディスパッチで、「2回目のtransposeで元にに戻る」性質が担保されているようになっていました。
コンストラクタにはこの機能を入れずにTransposeインスタンスを作っているだけという訳ですね。

Matrix{<:Number}以外に対する転置

文字列を要素に持つ行列の転置はできるでしょうか？

julia> A = ["hoge" "fuga";"piyo" "teke"]
2×2 Array{String,2}:
 "hoge"  "fuga"
 "piyo"  "teke"

julia> transpose(A)
2×2 Transpose{Union{},Array{String,2}}:
Error showing value of type Transpose{Union{},Array{String,2}}:
ERROR: MethodError: no method matching transpose(::String)
Closest candidates are:
  transpose(::Missing) at missing.jl:100
  transpose(::Number) at number.jl:168
  transpose(::Transpose) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/adjtrans.jl:165
  ...
Stacktrace:
 [1] getindex at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/adjtrans.jl:190 [inlined]

(略)

エラーになりましたね。

なぜこうなっているのか？

数を要素に持つ行列しかtransposeできない理由

transposeは行列のためのものです。
単なる多次元配列で添字を入れ替えたい場合はpermutedimsを使いましょう。

julia> permutedims(A)
2×2 Array{String,2}:
 "hoge"  "piyo"
 "fuga"  "teke"

broadcastの挙動

julia> [1,2]'*3
1×2 Adjoint{Int64,Array{Int64,1}}:
 3  6

julia> [1,2]'.*3
1×2 Array{Int64,2}:
 3  6

えっこれらは等価じゃないんですか？？

なぜこうなっているのか？

[1,2]'にbroadcacstを適用する場合に、先にArrayに変換されているようです。
だとしても、Adjoint{Int64,Array{Int64,1}}の型を保ってくれてる方が便利じゃないですか？？
分からない！！
(知ってる方、教えて下さい)

一応issue立ってるようでした↓↓

将来的にはもっと分かりやすくなってるかもですね。

参考文献など

JuliaCon 2017での解説動画↓

当時のissue↓

まとめ

行列・ベクトルに関する雑多な話題でしたが、私が疑問に思ったところは回収できた気がします。(最後を除いてですが)

ここで書いたように、内部的な振る舞いまで色々考えるとややこしい部分もありますが、実用上は@genkurokiさんの以下のツイートの通りだと思います。

「n次元横ベクトルのつもりで1×n行列を作らない！」と覚えておけば大抵の場合に自然にうまく行くと思います。

Juliaは他の言語に比べて、この辺りの行列・ベクトルの処理を上手く扱っていると思います。
他言語との具体的な比較まで書こうと考えてましたが間に合いませんでした。

おしまい

P.S.
間違い・疑問点・他の面白い例などあればコメントよろしくおねがいします！

GitHubで編集を提案

Julia Advent Calendar 2020/12/23

事の発端

行列とベクトルの簡単な導入

Juliaでの1次元配列は縦ベクトル

Juliaでの2次元配列は行列

転置・随伴

非自明なところ

横ベクトルや3×1行列は直接宣言できない

コード例

なぜこうなっているのか？

横ベクトルを直接宣言できない理由

3×1行列を直接宣言できない理由

2つのベクトルの積

コード例

なぜこうなっているのか？

ベクトル同士のナマの掛け算が定義されない理由

横ベクトルはベクトルではない

なぜこうなっているのか？

横ベクトルがAbstractVectorではない理由

横ベクトルがAbstractRowVectorではない理由

横ベクトルがAbstractMatrixである理由

添字の上下を考えない理由

番外: なんで横ベクトル？

collect∘adjoint、collect∘transposeの挙動

なぜこうなっているのか？

collect有りで1要素のベクトルが得られる理由

collect無しでスカラーが得られる理由

Transpose、Adjointの型

コード例

なぜこうなっているのか？

型の意味

型パラメータが冗長な理由 (冗長に見えるだけの理由)

関数transposeとコンストラクタTranspose

なぜこうなっているのか？

関数とコンストラクタの違い

Matrix{<:Number}以外に対する転置

なぜこうなっているのか？

数を要素に持つ行列しかtransposeできない理由

broadcastの挙動

なぜこうなっているのか？

参考文献など

まとめ

Discussion

横ベクトルが`AbstractVector`ではない理由

横ベクトルが`AbstractRowVector`ではない理由

横ベクトルが`AbstractMatrix`である理由

`collect`有りで1要素のベクトルが得られる理由

`collect`無しでスカラーが得られる理由