第10回：行列・線形代数

▶ ベクトルの線形結合からなる格子点

複数のベクトルの線形結合とは，それらのベクトルのスカラー倍を加え合わせたものを，それらのベクトルの線形結合という．

二つのベクトル $a_1=\begin{bmatrix} 1 \\\\ 0 \end{bmatrix}, a_2=\begin{bmatrix} 0 \\\\ 1 \end{bmatrix}$ の「整数」係数の線形結合による格子点を描く．

さらに， $b_1=\begin{bmatrix} \dfrac{1}{2} \\\\ \dfrac{1}{2} \end{bmatrix}, b_2=\begin{bmatrix} \dfrac{1}{2} \\\\ -\dfrac{1}{2} \end{bmatrix}$ の「整数」係数の線形結合からなる格子点を重ねる．

using PyPlot
plt.axes().set_aspect("equal")
a1 = [1, 0];
a2 = [0, 1];
for m = -3:3, n = -3:3
   r = m * a1 + n * a2
   plt.plot(r[1], r[2], "bo")
end

b1 = [1 / 2, 1 / 2];
b2 = [1 / 2, -1 / 2];
for m = -3:3, n = -3:3
   r = m * b1 + n * b2
   plt.plot(r[1], r[2], "r.")
end

plt.axhline(0, color = "k", lw = 0.5)
plt.axvline(0, color = "k", lw = 0.5)

どちらも正方格子（cubic lattice）であるが，座標系の取り方が異なる．基底 $b_1, b_2$ で張られる格子点は，基底 $a_1, a_2$ で張られる格子点の中央の点も含んでいることが観察できる．

今度は， $c_1=\begin{bmatrix} 1 \\ 0 \end{bmatrix}, c_2=\begin{bmatrix} -\dfrac{1}{2} \\ \dfrac{\sqrt{3}}{2} \end{bmatrix}$ で張られる格子点を描いてみる．これは，六方格子（hexagonal lattice）と呼ばれる．

using PyPlot
plt.axes().set_aspect("equal")

c1 = [1, 0];
c2 = [-1 / 2, sqrt(3) / 2];
for m = -3:3, n = -3:3
   r = m * c1 + n * c2
   plt.plot(r[1], r[2], "g.")
end

plt.xlim(-3, 3)
plt.ylim(-3, 3)
plt.axhline(0, color = "k", lw = 0.5)
plt.axvline(0, color = "k", lw = 0.5)

■ 平面ベクトルの内積

関数 dot(a,b) は，ベクトル a と b との「内積（inner product）」を返す．

中置き演算子 ⋅ を用いて a⋅b と書くこともできる．「 ⋅ 」は，バックスラッシュ \ に続けて cdots と入力してから，TABキーを押すことによって入力できる．かな漢字変換システムで入力できる「・」（中黒=なかぐろ）とは，別の文字である．

▶ 平面ベクトル同士のなす角を求める

ベクトル $a$ と $b$の内積は，$a$ と $b$ のなす角$\theta$ を用いて，以下のように定義される．

\[a\cdot b = \left\vert{a}\right\vert \left\vert{b}\right\vert \cos\theta\]

これから，$\theta$ を求めるには，次の式を用いればよい．

\[\cos\theta = \dfrac{a\cdot b}{ \left\vert{a}\right\vert \left\vert{b}\right\vert }\]

また，内積の定義から，自分自身の内積は，ノルムの二乗であることも分かる．

\[a\cdot a = \left\vert{a}\right\vert^2\]

▶ 例：ベクトル同士のなす角度を求める

上で出てきたベクトルのうち，a1, a2, c1, c2 のノルムは $1$ である．

julia> using LinearAlgebra

julia> #
       a1 = [1, 0];

julia> a2 = [0, 1];

julia> b1 = [1 / 2, 1 / 2];

julia> b2 = [1 / 2, -1 / 2];

julia> c1 = [1, 0];

julia> c2 = [-1 / 2, sqrt(3) / 2];

julia> norm(a1)
1.0

julia> norm(a2)
1.0

julia> norm(c1)
1.0

julia> norm(c2)
0.9999999999999999

b1 と b2 のノルムは $\dfrac{1}{\sqrt{2}}$ である．自分自身の内積の値と比較しよう

julia> b1 ⋅ b1
0.5

julia> norm(b1)
0.7071067811865476

julia> norm(b1)^2
0.5000000000000001

julia> b1 ⋅ b1
0.5

julia> norm(b2)
0.7071067811865476

julia> norm(b2)^2
0.5000000000000001

内積から算出した $\cos\theta$から角度 $\theta$ を得るには，関数 acos() を用いる．関数 acos(r) は $\cos \theta = r$ となる $\theta$ をラジアンで返す．関数 acosd(r) は，$\theta$ をラジアンで返す．

これらのベクトルのなす角度を算出しよう． a1 と a2,および，b1 と b2 は直交している． a1 と b1 は，45° をなしている． c1 と c2 は，120° をなしている，ことが計算できた．

julia> acosd(a1 ⋅ a2)
90.0

julia> acosd(b1 ⋅ b2 / norm(b1) / norm(b2))
90.0

julia> acosd(a1 ⋅ b1 / norm(a1) / norm(b1))
45.00000000000001

julia> acosd(c1 ⋅ c2)
120.00000000000001

■ タプル

タプル（tuple）は，複数の値をカンマ , で区切って並べ，括弧 () で囲んだものである．ベクトルと似たように使えるが，要素を更新することはできない．

julia> # 要素 1つのタプル
       (1,)
(1,)

julia> # 要素 2つのタプル
       (1, 2)
(1, 2)

julia> # 要素 3つのタプル
       a = (1, 2, 3)
(1, 2, 3)

julia> # タプルの長さ
       length(a)
3

julia> # タプルの要素
       a[2]
2

julia> # 更新はできない
       a[2] = 3 # => MethodError
ERROR: MethodError: no method matching setindex!(::Tuple{Int64,Int64,Int64}, ::Int64, ::Int64)

関数には，複数の値を返すものがある．このとき，タプルが用いられる．

例えば，divrem(x,d) は，div(x,d) と rem(x,d) の二つの値を返す．

julia> divrem(5, 3)
(1, 2)

タプルを右辺において，複数の変数に同時に代入できる．

julia> x, y = (1, 2, 3)
(1, 2, 3)

julia> x
1

julia> y
2

■ 行列

要素を ; で区切って列挙したものを，大かっこ [] で囲むと，行列を作ることができる．

julia> a = [11 12 13 14; 21 22 23 24; 31 32 33 34]
3×4 Array{Int64,2}:
 11  12  13  14
 21  22  23  24
 31  32  33  34

3行 4列の行列

■ ベクトルも，■ 行列も，「配列（array）」として表されている．ベクトルと同じ関数が用いられる．

julia> # 寸法 => タプル
       size(a)
(3, 4)

julia> # 第1軸 = 列の寸法
       size(a, 1)
3

julia> # 第2軸 = 行の寸法
       size(a, 2)
4

julia> # 全要素数
       length(a)
12

■ 行列の転置

transpose(a) は，行列 a を転置する．すなわち，行と列を入れ替える．

julia> a = [11 12 13 14; 21 22 23 24; 31 32 33 34]
3×4 Array{Int64,2}:
 11  12  13  14
 21  22  23  24
 31  32  33  34

julia> size(a)
(3, 4)

julia> b = transpose(a)
4×3 LinearAlgebra.Transpose{Int64,Array{Int64,2}}:
 11  21  31
 12  22  32
 13  23  33
 14  24  34

julia> size(b)
(4, 3)

▶ 行列のスカラー倍・スカラーの和差

以下，しばらく，2x2 の正方行列を例に説明する．

julia> a = [11 12; 21 22]
2×2 Array{Int64,2}:
 11  12
 21  22

julia> a * 2
2×2 Array{Int64,2}:
 22  24
 42  44

julia> a .+ 2
2×2 Array{Int64,2}:
 13  14
 23  24

julia> a .- 2
2×2 Array{Int64,2}:
  9  10
 19  20

▶ 行列に列ベクトルを加減

julia> a = [11 12; 21 22]
2×2 Array{Int64,2}:
 11  12
 21  22

julia> v = [1, -1]
2-element Array{Int64,1}:
  1
 -1

julia> a .+ v
2×2 Array{Int64,2}:
 12  13
 20  21

julia> a .- v
2×2 Array{Int64,2}:
 10  11
 22  23

▶ 行列同士の加減

julia> a = [11 12; 21 22]
2×2 Array{Int64,2}:
 11  12
 21  22

julia> b = a * 2
2×2 Array{Int64,2}:
 22  24
 42  44

julia> a + b
2×2 Array{Int64,2}:
 33  36
 63  66

julia> a - b
2×2 Array{Int64,2}:
 -11  -12
 -21  -22

■ 添字を用いた行列の要素の読み書き

行列の添字は，第1軸（列）と第2軸（行）の番号を，カンマ , で区切って並べ，大かっこ [] で囲んだものである．

ベクトルと同じように，添字で示された要素の読み出し，添字で示された要素の書き換えができる．

julia> # 添字による要素の読み出し
       a[2, 2]
22

julia> # 行列の要素の更新
       a[1, 2] = 30
30

julia> a
2×2 Array{Int64,2}:
 11  30
 21  22

■ 部分行列

julia> a = [11 12 13 14; 21 22 23 24; 31 32 33 34]
3×4 Array{Int64,2}:
 11  12  13  14
 21  22  23  24
 31  32  33  34

julia> # 列を取り出す
       a[:, 2]
3-element Array{Int64,1}:
 12
 22
 32

julia> # 行を取り出す
       a[2, :]
4-element Array{Int64,1}:
 21
 22
 23
 24

julia> # 部分行列
       a[1:2, 1:2]
2×2 Array{Int64,2}:
 11  12
 21  22

julia> a[2:3, 2:3]
2×2 Array{Int64,2}:
 22  23
 32  33

▶ 行列に入れた点座標で図形を描画する

using PyPlot
plt.axes().set_aspect("equal")

xy = [1 2 2 1; 1 1 3 1]
@show xy
plt.plot(xy[1, :], xy[2, :])
xy = xy .+ [1 / 2, 1 / 2]
plt.plot(xy[1, :], xy[2, :])
plt.xlim(0, 4)
plt.ylim(0, 4)

xy = [1 2 2 1; 1 1 3 1]

■ 行列とベクトルの積

julia> a = [11 12; 21 22]
2×2 Array{Int64,2}:
 11  12
 21  22

julia> v = [1, -1]
2-element Array{Int64,1}:
  1
 -1

julia> a * v
2-element Array{Int64,1}:
 -1
 -1

▶ 回転行列とベクトルの積

以下の形の行列を，平面空間に対する回転行列という．

\[R(\theta) = \begin{bmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos \theta \end{bmatrix}\]

回転行列とベクトルの積は，そのベクトルを，原点の周りに反時計方向に角 $\theta$ だけ回転する写像に対応する．

\[x^{\prime} = R(\theta) x\]

julia> # 回転行列
       r15 = [cosd(15) -sind(15); sind(15) cosd(15)]
2×2 Array{Float64,2}:
 0.965926  -0.258819
 0.258819   0.965926

julia> xy = [1, 0]
2-element Array{Int64,1}:
 1
 0

julia> xy = r15 * xy
2-element Array{Float64,1}:
 0.9659258262890683
 0.25881904510252074

julia> xy = r15 * xy
2-element Array{Float64,1}:
 0.8660254037844387
 0.49999999999999994

これらを描こう．軌跡は円を描いた．

using PyPlot
plt.axes().set_aspect("equal")

r15 = [cosd(15) -sind(15); sind(15) cosd(15)]
xy = [1, 0]

for i = 0:20
   global xy
   plt.plot(xy[1], xy[2], "o")
   xy = r15 * xy
end

plt.xlim(-1.2, 1.2)
plt.ylim(-1.2, 1.2)
plt.axhline(0, color = "k", lw = 0.5)
plt.axvline(0, color = "k", lw = 0.5)

原点以外の点 $c$ の周りで回転する場合は，回転の中心をずらして，

\[\begin{aligned} (x^{\prime}-c) & = R(\theta) (x-c) \\ x^{\prime} & = c + R(\theta) (x-c) \end{aligned}\]

とすればよい．

using PyPlot
plt.axes().set_aspect("equal")

r15 = [cosd(15) -sind(15); sind(15) cosd(15)]
xy = [1, 0]
c = [1 / 2, 1 / 2]

for i = 0:20
   global xy
   plt.plot(xy[1], xy[2], "o")
   xy = c + r15 * (xy - c)
end

plt.axvline(c[1], color = "k", lw = 0.5)
plt.axhline(c[2], color = "k", lw = 0.5)

plt.xlim(-1, 2)
plt.ylim(-1, 2)

■ 行列と行列の積

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

julia> b = [5 6; 7 8]
2×2 Array{Int64,2}:
 5  6
 7  8

julia> a * b
2×2 Array{Int64,2}:
 19  22
 43  50

▶ 座標を行列に格納した図形を回転する

using PyPlot
plt.axes().set_aspect("equal")
r15 = [cosd(15) -sind(15); sind(15) cosd(15)]
xy = [1 2 2 1; 1 1 3 1]

for i = 0:20
   global xy
   plt.plot(xy[1, :], xy[2, :])
   xy = r15 * xy
end

plt.xlim(-4, 4)
plt.ylim(-4, 4)
plt.axhline(0, color = "k", lw = 0.5)
plt.axvline(0, color = "k", lw = 0.5)

回転中心をずらしてみる

plt.axes().set_aspect("equal")

r15 = [cosd(15) -sind(15); sind(15) cosd(15)]
xy = [1 2 2 1; 1 1 3 1]
c = [1 / 2, 1 / 2]

for i = 0:20
   global xy
   plt.plot(xy[1, :], xy[2, :])
   xy = c .+ r15 * (xy .- c)
end

plt.axvline(c[1], color = "k", lw = 0.5)
plt.axhline(c[2], color = "k", lw = 0.5)

plt.xlim(-4, 4)
plt.ylim(-4, 4)

■ いろいろな行列の生成

■ 要素が 0 の行列を作る

Base.zeros — Function

関数 zeros は，要素が零 $0$ の行列を作るのに使える．

関数 zeros(n,m) は，要素の型が浮動小数点で，寸法 (n,m) の行列を作る．
関数 zeros(T, n,m) は，要素の型が T で，寸法 (n,m) の行列を作る．

julia> zeros(3, 4) # 要素は浮動小数点
3×4 Array{Float64,2}:
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0

julia> zeros(Float64, 3, 4) # 上と同じ
3×4 Array{Float64,2}:
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0

julia> zeros(Int64, 3, 4) # 要素は整数
3×4 Array{Int64,2}:
 0  0  0  0
 0  0  0  0
 0  0  0  0

行列 a と同じ寸法の 0 ベクトルを作るには，

julia> a = [11 12 13 14; 21 22 23 24; 31 32 33 34]
3×4 Array{Int64,2}:
 11  12  13  14
 21  22  23  24
 31  32  33  34

julia> zeros(Int64, size(a))
3×4 Array{Int64,2}:
 0  0  0  0
 0  0  0  0
 0  0  0  0

julia> zeros(Float64, size(a))
3×4 Array{Float64,2}:
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0

関数 zero.() を以下のように用いれば行列 a の要素の型と同じ要素の型を持ち，a と寸法が等しい 0 ベクトルを作れる．

julia> a = [11 12 13 14; 21 22 23 24; 31 32 33 34]
3×4 Array{Int64,2}:
 11  12  13  14
 21  22  23  24
 31  32  33  34

julia> zero.(a)
3×4 Array{Int64,2}:
 0  0  0  0
 0  0  0  0
 0  0  0  0

■ 要素が 1 の行列を作る

Base.ones — Function

関数 ones は，要素が零 $1$ の行列を作るのに使える．

関数 ones(n,m) は，要素の型が浮動小数点で，寸法 (n,m) の行列を作る．
関数 ones(T, n,m) は，要素の型が T で，寸法 (n,m) の行列を作る．

julia> ones(3, 4) # 要素は浮動小数点
3×4 Array{Float64,2}:
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0

julia> ones(Float64, 3, 4) # 上と同じ
3×4 Array{Float64,2}:
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0

julia> ones(Int64, 3, 4) # 要素は整数
3×4 Array{Int64,2}:
 1  1  1  1
 1  1  1  1
 1  1  1  1

行列 a と同じ寸法の 1 行列を作るには，

julia> a = [11 12 13 14; 21 22 23 24; 31 32 33 34]
3×4 Array{Int64,2}:
 11  12  13  14
 21  22  23  24
 31  32  33  34

julia> ones(Int64, size(a))
3×4 Array{Int64,2}:
 1  1  1  1
 1  1  1  1
 1  1  1  1

julia> ones(Float64, size(a))
3×4 Array{Float64,2}:
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0

関数 one.() を以下のように用いれば行列 a の要素の型と同じ要素の型を持ち，a と寸法が等しい 1 行列を作れる．

julia> a = [11 12 13 14; 21 22 23 24; 31 32 33 34]
3×4 Array{Int64,2}:
 11  12  13  14
 21  22  23  24
 31  32  33  34

julia> one.(a)
3×4 Array{Int64,2}:
 1  1  1  1
 1  1  1  1
 1  1  1  1

■ 対角要素を指定して，正方行列をつくる

julia> using LinearAlgebra

julia> Diagonal([1, 2, 3])
3×3 LinearAlgebra.Diagonal{Int64,Array{Int64,1}}:
 1  ⋅  ⋅
 ⋅  2  ⋅
 ⋅  ⋅  3

■ 疑似乱数を要素とする行列を作る

julia> rand(3, 3)
3×3 Array{Float64,2}:
 0.677818  0.930224    0.789515
 0.354494  0.00405935  0.594474
 0.933992  0.984617    0.932858

■ ベクトルをまとめて行列を作る

julia> a = [1, 2]
2-element Array{Int64,1}:
 1
 2

julia> b = [3, 4]
2-element Array{Int64,1}:
 3
 4

julia> c = [5, 6]
2-element Array{Int64,1}:
 5
 6

julia> [a b c]
2×3 Array{Int64,2}:
 1  3  5
 2  4  6

julia> m = hcat(a, b, c)
2×3 Array{Int64,2}:
 1  3  5
 2  4  6

julia> size(m)
(2, 3)

さらに転置をとろう．

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

julia> size(mt)
(3, 2)

● 内包表現を用いて行列を作る

julia> [[t, t .^ 2] for t in [0, 2, 4]]
3-element Array{Array{Int64,1},1}:
 [0, 0]
 [2, 4]
 [4, 16]

julia> hcat([[t, t .^ 2] for t in [0, 2, 4]]...)
2×3 Array{Int64,2}:
 0  2   4
 0  4  16

▶ 楕円を描く・回転する

using PyPlot
plt.axes().set_aspect("equal")

ts = 0:pi/18:2pi
xs = 2 * cos.(ts)
ys = sin.(ts)
xy = transpose(hcat(xs, ys))
plt.plot(xy[1, :], xy[2, :])

plt.xlim(-2.2, 2.2)
plt.ylim(-2.2, 2.2)

Note

楕円の座標を含む行列 xy を作るのに，以下のように● 内包表記を使ってもよい．

xy=hcat([ [2*cos.(t), sin(t)] for t=0:pi/18:2pi]...)

回転する

using PyPlot
plt.axes().set_aspect("equal")

r15 = [cosd(15) -sind(15); sind(15) cosd(15)]
ts = 0:pi/18:2pi
xy = transpose(hcat(2 * cos.(ts), sin.(ts)))

for i = 0:5
   global xy
   plt.plot(xy[1, :], xy[2, :])
   xy = r15 * xy
end

■ 行列式と逆行列

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

julia> det(a) # 行列式
-1.0

julia> inv(a) # 逆行列
2×2 Array{Float64,2}:
 -3.0   2.0
  2.0  -1.0

julia> a^(-1) # 逆行列
2×2 Array{Float64,2}:
 -3.0   2.0
  2.0  -1.0

■ 行列の商（2x2 行列）

行列 $A$とベクトル $x$ に対して，$x=Ay$ を満たすようなベクトル $y$を，$x$を$A$で割った商という．商を求める演算が用意されている．この演算は，$A$の逆行列を計算するよりも，計算コストの上で有利である．

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

julia> v = [1, 1]
2-element Array{Int64,1}:
 1
 1

julia> a \ v
2-element Array{Float64,1}:
 -1.0
  1.0

julia> w = [3, 5]
2-element Array{Int64,1}:
 3
 5

julia> a \ w
2-element Array{Float64,1}:
 1.0
 1.0

行列 $A$ が正則でない（逆行列が存在しない）場合には，例外が発生する．

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

julia> det(c)
0.0

julia> c \ v # => SingularException
ERROR: LinearAlgebra.SingularException(2)

同様に，行列 $X$ に対して，$X=AY$ を満たすような行列 $Y$ を，$X$ を $A$ で割った商という．

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

julia> a \ b
2×2 Array{Float64,2}:
 -1.0  1.0
  1.0  1.0

▶ 楕円を，逆に回転する

行列の商を用いて，楕円を逆回転してみよう．

using PyPlot
plt.axes().set_aspect("equal")

r15 = [cosd(15) -sind(15); sind(15) cosd(15)]
ts = 0:pi/18:2pi
xy = transpose(hcat(2 * cos.(ts), sin.(ts)))

for i = 0:5
   global xy
   plt.plot(xy[1, :], xy[2, :])
   xy = r15 \ xy
end

行列式が 0 の行列は，正則ではない

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

julia> det(a)
0.0

julia> v = [1, 1]
2-element Array{Int64,1}:
 1
 1

julia> # 例外を発生する
       a \ v
ERROR: LinearAlgebra.SingularException(2)

julia> # 例外を発生する
       inv(a)
ERROR: LinearAlgebra.SingularException(2)

▶ 空間ベクトル：なす角を求める

空間ベクトルとは，寸法 3 のベクトルである．内積が $0$ なら，それらのベクトルは直交である．

例: 以下の３つのベクトルが，互いに直行することを示せ．

julia> using LinearAlgebra

julia> #
       a = [1 / 2, 1 / 2 + sqrt(2) / 4, 1 / 2 - sqrt(2) / 4]
3-element Array{Float64,1}:
 0.5
 0.8535533905932737
 0.1464466094067262

julia> b = [-1 / 2, 1 / 2 - sqrt(2) / 4, 1 / 2 + sqrt(2) / 4]
3-element Array{Float64,1}:
 -0.5
  0.1464466094067262
  0.8535533905932737

julia> c = [1 / sqrt(2), -1 / 2, 1 / 2]
3-element Array{Float64,1}:
  0.7071067811865475
 -0.5
  0.5

julia> norm(a)
1.0

julia> norm(b)
1.0

julia> norm(c)
0.9999999999999999

julia> dot(a, b)
-5.551115123125783e-17

julia> a ⋅ b
-5.551115123125783e-17

julia> b ⋅ c
0.0

julia> c ⋅ a
-2.7755575615628914e-17

二つのベクトルのなす角を求めよ．

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

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

julia> ab = a ⋅ b
-7

julia> na = norm(a)
3.7416573867739413

julia> nb = norm(b)
3.7416573867739413

julia> r = ab / na / nb
-0.5

julia> # ラジアン単位
       acos(r)
2.0943951023931957

julia> # 角度単位
       acosd(r)
120.00000000000001

julia> # 一行で書ける
       acosd(a ⋅ b / norm(a) / norm(b))
120.00000000000001

■ 空間ベクトルの外積

関数 cross(a,b) は，空間ベクトル a と b との外積または「ベクトル積（outer product）」を返す．

中置き演算子 × を用いて a×b と書くこともできる． × は，バックスラッシュ \ に times と入力してから， TABキーを押すことによって入力できる．かな漢字変換システムで入力できる「✕」とは，別の文字である．

外積 $a\times b$ の向きは，$a$ から $b$ へ回転したとき，右ねじが進む方向である．

外積 $a\times b$ の大きさは，$a$ と $b$ のなす角 $\theta$ を用いて，以下のように定義される．これは，ベクトル a と b がなす平行四辺形の面積である．

\[\left\vert a\times b\right\vert = \left\vert{a}\right\vert \left\vert{b}\right\vert \sin\theta\]

▶ 空間座標の基本単位ベクトル

julia> a = [1, 0, 0]
3-element Array{Int64,1}:
 1
 0
 0

julia> b = [0, 1, 0]
3-element Array{Int64,1}:
 0
 1
 0

julia> c = [0, 0, 1]
3-element Array{Int64,1}:
 0
 0
 1

julia> cross(a, b)
3-element Array{Int64,1}:
 0
 0
 1

julia> # a×b = c
       a × b
3-element Array{Int64,1}:
 0
 0
 1

julia> # b×c = a
       b × c
3-element Array{Int64,1}:
 1
 0
 0

julia> # c×a = b
       c × a
3-element Array{Int64,1}:
 0
 1
 0

別の正規直交系の例

julia> a = [1 / 2, 1 / 2 + sqrt(2) / 4, 1 / 2 - sqrt(2) / 4]
3-element Array{Float64,1}:
 0.5
 0.8535533905932737
 0.1464466094067262

julia> b = [-1 / 2, 1 / 2 - sqrt(2) / 4, 1 / 2 + sqrt(2) / 4]
3-element Array{Float64,1}:
 -0.5
  0.1464466094067262
  0.8535533905932737

julia> c = [1 / sqrt(2), -1 / 2, 1 / 2]
3-element Array{Float64,1}:
  0.7071067811865475
 -0.5
  0.5

julia> # a×b = c
       a × b
3-element Array{Float64,1}:
  0.7071067811865475
 -0.5
  0.5

julia> # b×c = a
       b × c
3-element Array{Float64,1}:
 0.5
 0.8535533905932737
 0.14644660940672627

julia> # c×a = b
       c × a
3-element Array{Float64,1}:
 -0.5
  0.14644660940672627
  0.8535533905932737

ベクトル３重積

３つの空間ベクトルに対して，一般に，以下が成り立つ．

\[{a}\times ({b} \times {c}) = ({a}\cdot{c}){b} - ({a}\cdot {b}){c}\]

例：具体的なベクトルの例で，上式が成り立つことを示せ．

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

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

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

julia> # 左辺
       lhs = a × (b × c)
3-element Array{Int64,1}:
  -7
  35
 -28

julia> # 右辺
       rhs = (a ⋅ c) * b - (a ⋅ b) * c
3-element Array{Int64,1}:
  -7
  35
 -28

▶ 行列の商（3x3 行列）

行列 A と行列（またはベクトル） Yに対して，商 A\Y は，$AX-Y$の最小二乗ノルムが最小となる行列（またはベクトル） X を返す．行列 A が正則なら，A の逆行列を左から Y に乗じた行列（またはベクトル）と一致する．

julia> #
       a = [-3, 1, 2]
3-element Array{Int64,1}:
 -3
  1
  2

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

julia> w = [a b]
3×2 Array{Int64,2}:
 -3   2
  1  -3
  2   1

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

julia> v = w \ c
2-element Array{Float64,1}:
 -1.0000000000000004
 -1.0000000000000002

julia> w * v - c
3-element Array{Float64,1}:
  8.881784197001252e-16
  4.440892098500626e-16
 -8.881784197001252e-16

▶ 行列の階数（ランク）

関数 rank(a) は，行列 a の階数（ランク，rank）を返す．

julia> a = [1 / 2, 1 / 2 + sqrt(2) / 4, 1 / 2 - sqrt(2) / 4]
3-element Array{Float64,1}:
 0.5
 0.8535533905932737
 0.1464466094067262

julia> b = [-1 / 2, 1 / 2 - sqrt(2) / 4, 1 / 2 + sqrt(2) / 4]
3-element Array{Float64,1}:
 -0.5
  0.1464466094067262
  0.8535533905932737

julia> c = [1 / sqrt(2), -1 / 2, 1 / 2]
3-element Array{Float64,1}:
  0.7071067811865475
 -0.5
  0.5

julia> v = [a b c]
3×3 Array{Float64,2}:
 0.5       -0.5        0.707107
 0.853553   0.146447  -0.5
 0.146447   0.853553   0.5

julia> rank(v)
3

julia> #
       a = [-3, 1, 2]
3-element Array{Int64,1}:
 -3
  1
  2

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

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

julia> v = [a b c]
3×3 Array{Int64,2}:
 -3   2   1
  1  -3   2
  2   1  -3

julia> rank(v)
2

▼ 行列の固有値・固有ベクトル

\[Ax = \lambda x\]

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

julia> #
       using LinearAlgebra

julia> # 固有値
       lam = eigvals(a)
2-element Array{Float64,1}:
 2.0
 5.0

julia> # 固有ベクトル
       evs = eigvecs(a)
2×2 Array{Float64,2}:
 -0.447214  0.707107
  0.894427  0.707107

julia> evs[:, 1]
2-element Array{Float64,1}:
 -0.4472135954999579
  0.8944271909999159

julia> evs[:, 2]
2-element Array{Float64,1}:
 0.7071067811865475
 0.7071067811865475

julia> # 確認する a x - lambda x == 0 になるべき
       a * evs[:, 1] - lam[1] * evs[:, 1]
2-element Array{Float64,1}:
 0.0
 2.220446049250313e-16

julia> a * evs[:, 2] - lam[2] * evs[:, 2]
2-element Array{Float64,1}:
 0.0
 0.0

★ 今回のまとめ

ベクトルの内積
行列の生成
行列に対する関数
行列とベクトルの演算
行列と行列の演算
部分行列
2次元の回転行列