「線形代数」とは何でしょうか?
技術的な定義は、[Wikipedia][linalg-wiki]や教科書、[3Blue1Brown][3blue1brown]のような優れたウェブサイトにたくさんあります。ただ、ここでの目的には、もっとくだけた説明で十分です。
線形代数とは、ベクトルと行列でできる、たくさんのおもしろい(そしてとても役に立つ)ことについての学問です。
Julia(とPython)のメンテナーはこういう話が大好きです。 Exercismの運営陣は、それほどでもありません。
このコンセプトは、巨大なテーマの中でもより簡単な側面だけに限定します。ただし、警告しておきます。これは必然的に、かなり数学的なコンセプトです。
誰に聞くか、そしてかなり抽象的なものをどう可視化したいかによって変わります。
Vector{T}はArray{T, 1}の別名にすぎず、eltypeであるTにはどんな型でも入ります。lengthとdirectionを持つ矢印ですが、位置は固定されていません。pointとして表せます。N個の要素が、N本の各軸に沿った原点からの距離を表します。ここでの目的のためには、文字列や文字のベクトルは無視してください。
このコンセプトでは、数値型を扱います。Int、Float、Complexです。
ここでも、さまざまな答えがあります。
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
上三角行列: 対角成分とその上側にゼロでない値があり、下側はゼロです。 下三角行列は想像がつくでしょう。
本当に行と列を入れ替えたいなら、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の行を書けば名前空間に読み込まれますが、例の中ではいちいち繰り返しません(見た目がごちゃごちゃしすぎるので)。
これは[Vector Operations][vector-ops]コンセプトで説明しました。
演算子は.*で、入力ベクトルに対して要素のペアごとに演算し、入力と同じサイズと型の出力を返します。
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)は失敗します。
ベクトルはどれくらい「大きい」のでしょうか?
normは、ベクトルを適切なスカラーに還元することで、これを捉えようとするものです。
ノルムにはさまざまな種類が定義されていますが、圧倒的によく使われるのは2ノルムで、√(v ⋅ v)です。
この二乗平均平方根の演算は、(N次元空間における)原点からの_ピタゴラス距離_です。 ベクトルを、尾を原点に置いた矢印として思い描くなら、2ノルムは矢印の_長さ_です。
pを2番目の引数として渡せば、任意のpノルムを計算できます。
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年の大半を、あらゆる計算を行列の乗算の連続に置き換えることに費やしてきたように思えることがあります。
これはコンピューターがとても得意とする演算です。
細かい仕組みはかなりシンプルですが、一見するとあまり直感的ではありません。
行列Aがベクトルvを掛けて出力wを生み出す状況を考えましょう(線形代数の慣習に従い、行列には大文字、ベクトルには小文字を使います)。
Aの1番上の行と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の1番上の行とBの各列のドット積がCの1番上の行になり、2番目の行が2番目の行になり、以下同様に続く、と言うこともできます。
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つのベクトルのドット積は、それらの長さが等しいことに依存します。
ひいては、行列の乗算では、左の行列の_列の数_が右の行列の_行の数_と一致している必要があります。
サイズをsize(A)が返すような(nrows, ncols)タプルで表すと、(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とは異なります)。
別の方法として、次元(3, 1) * (1, 3) -> (3, 3)のu * v'を使えば、可能なすべての要素のペアを掛けて、ベクトルを行列に広げられます。
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は、_内積_と呼ばれることもあります。
特に一般的な行列の乗算の一種に、回転行列を使うものがあります。
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
ロボット工学のスタートアップで働いています。その会社は、概念実証として簡単なロボットを開発しています。ここでは、そのロボットの動きを制御する機能をいくつか実装する任務を任されました。
ロボットの向きと広がりを把握するために、ロボットには中心から1単位の距離にある3つのマーカーがあります。その位置を初期化するには、3つの方向ベクトルを正規化して、行列にまとめる必要があります。
orientrobot(vectors)関数を実装してください。この関数は、3つのベクトルからなるベクトルを受け取ります。正規化したベクトルを列に持つ2x3行列を返してください。
julia> orientrobot([[-1,1],[1,0],[-1,-1]])
2×3 Matrix{Float64}:
-0.707107 1.0 -0.707107
0.707107 0.0 -0.707107
次に、進む方向を変える機能が必要です。そのためには、ロボットが進みたい方向を向くように回転させる必要があります。
rotaterobot(orientation, θ)関数を実装してください。この関数は、ロボットの向きを表す行列と、反時計回りに回転させる角度θを受け取ります。新しい向きの行列を返してください。
julia> orientmatrix = initialize([[-1,1],[1,0],[-1,-1]]);
julia> rotaterobot(orientmatrix, π/2)
2×3 Matrix{Float64}:
-0.707107 6.12323e-17 0.707107
-0.707107 1.0 -0.707107
ロボットをある場所から別の場所へ動かすには、まず動かす前に向きが正しいかどうかを確認する必要があります。向きを表す行列の2列目が、前を向く方向を表しています。
robotoriented(orientation, direction)関数を実装してください。この関数は、向きを表す行列と相対位置ベクトルを受け取ります。ロボットが相対位置ベクトルと同じ方向を向いていれば(丸め誤差の範囲内で)trueを返してください。
julia> orientmatrix = initialize([[-1,1],[1,0],[-1,-1]]);
julia> robotoriented(orientmatrix, [5, 0])
true
julia> robotoriented(orientmatrix, [0, 5])
false
julia> robotoriented(orientmatrix, [-5, 0])
false
向きを表す行列はロボットの形状も記録しているので、ロボットを動かしたあとに、これらの点が原点から見てどこにあるのかを知る必要があります。これは、ロボットが移動するときに他の物体との衝突を避けるのに役立ちます。
bodylocation(orientation, position)関数を実装してください。この関数は、向きを表す行列と、ロボットの中心の現在位置を受け取ります。平行移動した向きの行列を返してください。
julia> orientmatrix = initialize([[-1,1],[1,0],[-1,-1]]);
julia> bodylocation(orientmatrix, [5, 3])
2×3 Matrix{Float64}:
4.29289 6.0 4.29289
3.70711 3.0 2.29289