Juliaと線形代数

3. Juliaと線形代数#

第2章では、用意された減衰計算のコードをファイルに保存し、実行結果をログに残した。本章では、そのコードを読み、自分で条件や計算内容を変えられるようになることを目指す。その後、複数の数をまとめたベクトル・行列を扱う方法へ進む。

なぜ数値解析で線形代数が必要なのだろうか。例えば、複数の測定値を並べるとベクトルになり、複数の未知数の関係をまとめると行列が現れる。後の章では、連立方程式や回帰、微分方程式の近似などを、ベクトルと行列を使って計算する。本章で、小さな例の意味を手計算とコードの両方で確かめておこう。

3.1. この章の目標#

第1・2章の更新式と、Juliaの変数・関数・繰り返しを対応づけられる。
ベクトルと行列を作り、成分と大きさを確認できる。
行列積と要素ごとの演算の違いを、手計算の例で説明できる。
ベクトルの長さや内積を計算し、連立方程式の答えを代入して確認できる。
後の章のコードに現れるキーワード引数・無名関数・名前付きタプルなどの書き方を読める。

線形代数の基礎を復習したい場合は、AI・機械学習論1の線形代数の章も参照してほしい。ここでは実数のベクトルと行列を中心に扱い、必要な式を具体例とともに確認する。

3.2. コードを試す準備#

第2章で使ったJuliaを引き続き利用する。本章の基本部分は、Julia本体と、同梱されている標準ライブラリ LinearAlgebra で実行できる。追加パッケージは必要ない。未導入の場合はJulia公式インストール案内を参照する。

実行方法	本章のコードを試すには
`.jl` ファイル	コードを上から順に一つのファイルへ保存し、第2章と同様に `julia ファイル名.jl` で実行する
Juliaの対話画面（REPL）	ターミナルで `julia` を起動し、`julia>` の後へJuliaのコードを入力する
Jupyter Notebookなど	Juliaを実行するカーネル・実行環境を選び、セルを上から順に実行する

ターミナルのシェルに対するコマンドと、Juliaのコードは区別しよう。本章のコードセルはJulia用である。複数行の関数などは、途中の1行だけでなく、対応する end までまとめて入力する。

Notebookでは、実行済みの変数や関数が残っている。上のセルを書き換えたら、そのセルと、それに依存する後のセルを実行し直す。結果がおかしい場合は、カーネルを再起動して先頭から実行すると、古い変数への依存に気づけることがある。Colabなどを利用する場合も、コードを実行する言語がJuliaであることを確認する。

本章では println で結果を明示的に表示する。対話画面では式の値が自動表示される場合があるが、ファイルとして実行した場合も出力が分かるように、この書き方を使う。

3.3. 変数と計算式：まず1ステップを読む#

第1章のEuler法は、現在の量 $N$ から短時間後の量を

\[ N_{\mathrm{next}}=N-h\lambda N \]

で求めた。第2章と同じく、$\lambda$ は入力しやすい名前 lambda を使って表す。$N=1$、$h=0.5$、$\lambda=0.4$ なら、次の量は $1-0.5\times0.4\times1=0.8$ になる。

lambda = 0.4
h = 0.5
N = 1.0
Nnext = N - h*lambda*N
println("現在の量 = ", N)
println("次の量 = ", Nnext)

現在の量 = 1.0
次の量 = 0.8

= は右辺の値を左辺の名前に代入する記号である。数学の等式とは役割が異なる。例えば N = N - h*lambda*N と書けば、右辺を更新前の値で計算してから、N が指す値を更新する。上の例では Nnext という別の名前に結果を入れたので、元の N は1のままである。等しいかを調べる比較には == を使う。

掛け算には *、累乗には ^ を使う。割り算 / の結果は、整数同士でも通常は浮動小数点数になる。Pythonで累乗に使う ** との違いに注意しよう。

3.3.1. 型は、どのような値かを表す#

整数 3 と小数 3.0 は数学的には等しい値だが、Juliaでは異なる型を持つ。typeof で型を確認できる。

count = 3
value = 3.0
label = "decay"
println("3の型 = ", typeof(count))
println("3.0の型 = ", typeof(value))
println("文字列の型 = ", typeof(label))
println("3/2 = ", 3/2)
println("3^2 = ", 3^2)

3の型 = Int64
3.0の型 = Float64
文字列の型 = String
3/2 = 1.5
3^2 = 9

一般的な64ビット環境では、型はそれぞれ Int64（整数）、Float64（倍精度浮動小数点数）、String（文字列）となる。通常、変数を使うたびに型を宣言する必要はなく、代入する値から決まる。Float64 も有限桁で数を表すため、第1章で見た丸めが起こることを思い出そう。

3.4. 関数と繰り返し：入力を変えて同じ計算を使う#

同じ更新を何度も書く代わりに、入力を受け取って計算結果を返す関数にまとめる。例えば、1ステップ分の計算は次のように書ける。

function decay_step(N, lambda, h)
    return N - h*lambda*N
end

println("1ステップ後 = ", decay_step(1.0, 0.4, 0.5))
println("もう1ステップ後 = ", decay_step(0.8, 0.4, 0.5))

1ステップ後 = 0.8
もう1ステップ後 = 0.64

function の後ろに関数名と入力の名前（引数）を書き、return で結果を返す。end は関数の定義の終わりである。関数を定義することと、decay_step(1.0, 0.4, 0.5) のように呼び出して計算することは別の操作である。

この場合、関数内の N、lambda、h は呼び出し時に渡した値を使う。外側に同じ名前の変数があっても、それを使っているわけではない。必要な値を引数として渡すと、計算条件が明確になり、同じ関数を別の条件でも使える。

3.4.1. forで更新を繰り返す#

第2章の main は、関数の中で初期量やステップ数を固定していた。今度はそれらを引数として受け取り、条件を変えて使えるようにする。次は第1章と同じEuler法の関数を、条件分岐も明示して書いたものである。

function euler_decay(lambda, N0, T, n::Integer)
    if n <= 0
        error("ステップ数nは正の整数にしてください")
    end
    if T <= 0
        error("最終時刻Tは正にしてください")
    end
    h = T/n
    N = N0
    for k in 1:n
        N = decay_step(N, lambda, h)
    end
    return N
end

for n in (5, 10, 20)
    result = euler_decay(0.4, 1.0, 5.0, n)
    println("n=", n, ": N(T)=", result)
end

n=5: N(T)=0.07776
n=10: N(T)=0.10737418240000003
n=20: N(T)=0.12157665459056928

結果は第1章と同じく、約0.077760、0.107374、0.121577となる。関数内の各部分を読み解こう。

記述	意味
`n::Integer`	この引数として整数を受け取る指定
`if n <= 0 ... end`	条件が成り立つときだけ中を実行する。ここでは `error` で理由を表示して処理を止める
`h = T/n`	比較する最終時刻をそろえるため、ステップ数から刻み幅を決める
`for k in 1:n ... end`	`k` を1から `n` まで変えながら、中の更新を `n` 回行う
`return N`	最後に得られた量を関数の結果として返す

外側の for n in (5, 10, 20) は、三つの条件を順に試すための繰り返しである。Juliaでは function、for、if などを end で閉じる。字下げだけで範囲が決まるPythonと異なるが、どこに何が含まれるか分かるように字下げして書く。

この例のように、途中の値を書き換えながら進む計算を関数内にまとめると、変数の使われる範囲が明確になる。ファイルの最上位やNotebookに直接、外側の変数を書き換えるループを置くと、実行環境による変数の扱いの違いに戸惑うことがある。まずは関数の入力と出力を明らかにする書き方に慣れよう。

3.5. 複数の数をまとめる：ベクトルと配列#

これまでは量 $N$ を一つの数で表した。複数地点の温度や複数の未知数を扱うには、数を順に並べたベクトルが便利である。例えば

\[\begin{split} \mathbf{x}=\begin{bmatrix}3\\4\end{bmatrix} \end{split}\]

を、Juliaでは x = [3.0, 4.0] と書く。この x は、順序を持つ数の入れ物である1次元配列で、線形代数の計算では列ベクトルとして扱える。

x[1] は第1成分、x[2] は第2成分である。Juliaの標準的な配列の添字は1から始まる。Pythonの配列のように最初を x[0] とすると、この配列では範囲外になる。

x = [3.0, 4.0]
println("x = ", x)
println("成分数 = ", length(x))
println("第1成分 = ", x[1])
println("第2成分 = ", x[2])
println("成分の和 = ", sum(x))

x = [3.0, 4.0]
成分数 = 2
第1成分 = 3.0
第2成分 = 4.0
成分の和 = 7.0

length(x) は成分数で、この場合は2である。「ベクトルの幾何学的な長さ」とは異なり、そちらは後で扱うノルムで測る。sum(x) は成分の和 $3+4=7$ を求める。

成分を書き換えるには、例えば x[1] = 10.0 と代入する。元の配列を残したまま試すには copy を使う。

x_copy = copy(x)
x_copy[1] = 10.0
println("元の配列 = ", x)
println("変更したコピー = ", x_copy)

元の配列 = [3.0, 4.0]
変更したコピー = [10.0, 4.0]

x_copy = x と代入するだけでは配列の中身はコピーされず、二つの名前で同じ配列を参照することになる。元データを残したい場合には、この違いに注意する。ここで扱った数値の配列なら copy(x) で独立した配列を作れる。

3.6. 行列の成分と形を確認する#

行列は、数を行と列に並べたものである。例えば

\[\begin{split} \mathbf{A}=\begin{bmatrix}1&2\\0&1\end{bmatrix} \end{split}\]

は2行2列の行列で、Juliaでは A = [1.0 2.0; 0.0 1.0] と書く。同じ行の成分は空白、行の区切りはセミコロンで表す。ベクトルの [3.0, 4.0] に使ったカンマと役割が異なる。

A[i,j] は第 $i$ 行、第 $j$ 列の成分である。size(A) は「行数、列数」を返す。次の出力を先に予想してから実行しよう。

A = [1.0 2.0; 0.0 1.0]
println("A = ")
display(A)
println("Aの形 = ", size(A))
println("A[1,2] = ", A[1,2])
println("Aの第1行 = ", A[1,:])
println("Aの第2列 = ", A[:,2])
println("転置した行列 = ")
display(transpose(A))

A = 
2×2 Matrix{Float64}:
 1.0  2.0
 0.0  1.0
Aの形 = (2, 2)
A[1,2] = 2.0
Aの第1行 = [1.0, 2.0]
Aの第2列 = [2.0, 1.0]
転置した行列 = 
2×2 transpose(::Matrix{Float64}) with eltype Float64:
 1.0  0.0
 2.0  1.0

size(A) は (2, 2)、A[1,2] は2となる。: はその方向の成分をすべて選ぶ指定なので、A[1,:] は第1行、A[:,2] は第2列である。transpose(A) は行と列を入れ替えた転置行列を返す。display は行列を見やすい形で表示するために使っている。

Juliaのベクトル x の size(x) は (2,) であり、2行1列の2次元配列 (2, 1) とは区別される。ただし、次に見る A*x では数学の列ベクトルと同じように計算できる。

3.6.1. 行列をベクトルにかける#

行列積は、各行とベクトルの対応する成分を掛けて足す計算である。上の $\mathbf{A}$ と $\mathbf{x}=(3,4)^\mathsf{T}$ では、

\[\begin{split} \mathbf{A}\mathbf{x} =\begin{bmatrix}\textcolor{#0072B2}{1}&\textcolor{#0072B2}{2}\\0&1\end{bmatrix} \begin{bmatrix}3\\4\end{bmatrix} =\begin{bmatrix}\textcolor{#0072B2}{1\times3+2\times4}\\0\times3+1\times4\end{bmatrix} =\begin{bmatrix}\textcolor{#0072B2}{11}\\4\end{bmatrix} \end{split}\]

となる。青色は「行列の第1行 → 掛けて足す計算 → 第1出力」の対応を示す。第1出力は二つの入力を $1:2$ の重みで足したもの、第2出力は第2入力をそのまま取り出したものである。行列は、このように複数の数を組み合わせて新しい数の組へ変換する働きを持つ。

println("A*x = ", A*x)
println("第1成分を手順通りに計算 = ", A[1,1]*x[1] + A[1,2]*x[2])

A_rect = [1.0 2.0 0.0; 0.0 1.0 1.0]
x_three = [3.0, 4.0, 5.0]
println("長方形の行列の形 = ", size(A_rect))
println("A_rect*x_three = ", A_rect*x_three)

A*x = [11.0, 4.0]
第1成分を手順通りに計算 = 11.0
長方形の行列の形 = (2, 3)
A_rect*x_three = [11.0, 9.0]

2行3列の行列と3成分のベクトルなら、結果は2成分となり、上の例では $(11,9)^\mathsf{T}$ である。一般に $m$ 行 $n$ 列の行列は、$n$ 成分のベクトルにかけると $m$ 成分のベクトルを返す。内側の大きさが一致するか、結果の大きさはいくつかを、計算前に確認しよう。

3.7. 行列積と要素ごとの演算を区別する#

行列同士の積も、左の行列の各行と、右の行列の各列の成分を掛けて足す計算である。例えば

\[\begin{split} \mathbf{B}=\begin{bmatrix}1&2\\3&4\end{bmatrix},\qquad \mathbf{C}=\begin{bmatrix}2&0\\-1&1\end{bmatrix} \end{split}\]

なら、積の第1行第1列は $1\times2+2\times(-1)=0$ になる。すべての成分を計算すると

\[\begin{split} \mathbf{B}\mathbf{C}=\begin{bmatrix}0&2\\2&4\end{bmatrix} \end{split}\]

である。Juliaでは B*C と書く。一般に、$m$ 行 $n$ 列と $n$ 行 $p$ 列の行列の積は $m$ 行 $p$ 列になる。

一方、同じ位置の成分だけを掛けたい場合は B .* C と書く。上の例なら結果は

\[\begin{split} \begin{bmatrix}1\times2&2\times0\\3\times(-1)&4\times1\end{bmatrix} =\begin{bmatrix}2&0\\-3&4\end{bmatrix} \end{split}\]

となり、行列積とは異なる。同じ形のデータ同士を成分ごとに処理する場合などに用いる。

B = [1.0 2.0; 3.0 4.0]
C = [2.0 0.0; -1.0 1.0]
println("行列積 B*C = ")
display(B*C)
println("要素ごとの積 B.*C = ")
display(B .* C)

行列積 B*C = 
2×2 Matrix{Float64}:
 0.0  2.0
 2.0  4.0
要素ごとの積 B.*C = 
2×2 Matrix{Float64}:
  2.0  0.0
 -3.0  4.0

3.7.1. ドットを付けて各要素に適用する#

配列の各要素に同じ操作を行う書き方をbroadcastという。x .^ 2 は各成分を二乗し、2 .* x .+ 1 は各成分に「2倍して1を足す」を適用する。関数の場合は exp.(配列) のように関数名の後ろにドットを付ける。

例えば、複数の時刻で減衰の解析解 $e^{-0.4t}$ を計算してみよう。times に時刻を並べ、各時刻へ同じ式を適用する。

println("各成分の二乗 = ", x .^ 2)
println("各成分を2倍して1を足す = ", 2 .* x .+ 1)
times = [0.0, 1.0, 2.0]
exact_values = exp.(-0.4 .* times)
println("時刻 = ", times)
println("解析解 = ", exact_values)

各成分の二乗 = [9.0, 16.0]
各成分を2倍して1を足す = [7.0, 9.0]
時刻 = [0.0, 1.0, 2.0]
解析解 = [1.0, 0.6703200460356393, 0.44932896411722156]

最初の二つの結果は [9.0, 16.0] と [7.0, 9.0] である。解析解はおよそ [1.0, 0.670320, 0.449329] となり、第1章の値とも対応する。複数時刻の値をまとめて処理する場面でも、配列が役立つことが分かる。

なお、B .* C のような同じ形同士だけでなく、上の例の 0.4 .* times のように一つの数を全成分へ使うこともできる。まずは「行列積を行うのか、各成分に同じ操作を行うのか」を数式から判断しよう。

3.8. ベクトルの内積と長さを測る#

二つの実ベクトルの対応する成分を掛けて足したものを内積という。

\[ \mathbf{x}^\mathsf{T}\mathbf{y}=\sum_{i=1}^n x_i y_i \]

例えば $\mathbf{x}=(3,4)^\mathsf{T}$、$\mathbf{y}=(-4,3)^\mathsf{T}$ なら、内積は $3\times(-4)+4\times3=0$ である。平面の矢印と考えると、この二つは直交している。

また、$\mathbf{x}$ の長さは三平方の定理から $\sqrt{3^2+4^2}=5$ である。これを一般の $n$ 成分へ広げた

\[ \|\mathbf{x}\|_2=\sqrt{\sum_{i=1}^n x_i^2} \]

を**2ノルム（ユークリッドノルム）**という。複数の成分を持つベクトルの大きさを、一つの非負の数で表せる。

これらを計算する関数は、標準ライブラリ LinearAlgebra に用意されている。using LinearAlgebra はその機能を読み込む命令で、追加インストールは不要である。

using LinearAlgebra

x = [3.0, 4.0]
y = [-4.0, 3.0]
println("内積 = ", dot(x, y))
println("成分数 length(x) = ", length(x))
println("長さ norm(x) = ", norm(x))
println("三平方の定理で計算 = ", sqrt(x[1]^2 + x[2]^2))

内積 = 0.0
成分数 length(x) = 2
長さ norm(x) = 5.0
三平方の定理で計算 = 5.0

内積は0、成分数は2、ノルムは5となる。dot(x,y) はここで定義した内積、ベクトルに対する norm(x) は2ノルムを計算する。

二つのベクトルの差にノルムを使えば、点の間の距離を測れる。真の値を $\mathbf{x}$、近似値を $\hat{\mathbf{x}}$ とすると、$\|\hat{\mathbf{x}}-\mathbf{x}\|_2$ は誤差ベクトルの大きさである。

x_near = [3.1, 3.9]
error_vector = x_near - x
println("成分ごとの差 = ", error_vector)
println("差の大きさ = ", norm(error_vector))

成分ごとの差 = [0.10000000000000009, -0.10000000000000009]
差の大きさ = 0.14142135623730964

各成分の差はおよそ $(0.1,-0.1)^\mathsf{T}$ で、その大きさは $\sqrt{0.1^2+(-0.1)^2}\approx0.1414$ となる。ここでも表示される小数の末尾には丸めの影響がある。

ベクトルの大きさの測り方は2ノルム以外にもあり、行列を測るノルムもある。それらの定義や誤差評価での使い分けは第4章で扱う。特に、Juliaの norm(A) と opnorm(A) は行列に対して異なる量を測るので、ベクトルの norm(x) の説明をそのまま行列へ当てはめないようにしよう。本章ではベクトルの長さと距離を理解できればよい。

3.9. 連立一次方程式を解き、答えを確かめる#

次の連立方程式を考える。

\[\begin{split} \begin{aligned} 2x_1+x_2&=5,\\ x_1+3x_2&=5. \end{aligned} \end{split}\]

第1式から $x_2=5-2x_1$ として第2式へ代入すれば、$x_1=2$、$x_2=1$ と分かる。この問題は、係数、未知数、右辺をまとめて

\[\begin{split} \underbrace{\begin{bmatrix}2&1\\1&3\end{bmatrix}}_{\mathbf{A}} \underbrace{\begin{bmatrix}x_1\\x_2\end{bmatrix}}_{\mathbf{x}} =\underbrace{\begin{bmatrix}5\\5\end{bmatrix}}_{\mathbf{b}} \end{split}\]

と書ける。行列とベクトルの積を成分ごとに展開すると、元の二つの式に戻る。

Juliaでは、係数行列 A_system と右辺 b を用意し、A_system \ b で解を求められる。バックスラッシュ \ は連立方程式を解く演算であり、成分ごとの割り算ではない。この例では係数行列が正則であり、解が一つに定まる。

A_system = [2.0 1.0; 1.0 3.0]
b = [5.0, 5.0]
xhat = A_system \ b
residual = b - A_system*xhat
println("解 = ", xhat)
println("左辺に代入した値 = ", A_system*xhat)
println("残差ベクトル = ", residual)
println("残差の大きさ = ", norm(residual))

解 = [2.0, 1.0]
左辺に代入した値 = [5.0, 5.0]
残差ベクトル = [0.0, 0.0]
残差の大きさ = 0.0

手計算と同じく、解はおよそ $(2,1)^\mathsf{T}$ になる。求めた解を左辺へ代入した値が右辺 $(5,5)^\mathsf{T}$ に近いかを確認するため、

\[ \mathbf{r}=\mathbf{b}-\mathbf{A}\hat{\mathbf{x}} \]

を計算した。これが残差ベクトルである。各成分は、各方程式の右辺と左辺の差を表す。norm(residual) は、それらの差を一つの大きさにまとめている。

比較のため、少し違う候補 $\mathbf{x}_{\mathrm{trial}}=(2.1,0.9)^\mathsf{T}$ も代入してみよう。

x_trial = [2.1, 0.9]
trial_residual = b - A_system*x_trial
println("候補を左辺に代入 = ", A_system*x_trial)
println("候補の残差 = ", trial_residual)
println("候補の残差の大きさ = ", norm(trial_residual))

候補を左辺に代入 = [5.1000000000000005, 4.8]
候補の残差 = [-0.10000000000000053, 0.20000000000000018]
候補の残差の大きさ = 0.22360679774997935

左辺は約 $(5.1,4.8)^\mathsf{T}$、残差は約 $(-0.1,0.2)^\mathsf{T}$ となる。その大きさは約0.224であり、この候補が方程式を満たしていないことが分かる。

数学的には、正則な行列に対する解を $\mathbf{x}=\mathbf{A}^{-1}\mathbf{b}$ と書ける。しかし、コードで逆行列を明示的に計算する必要はない。A_system \ b は、行列に応じた方法を使って方程式を解く。計算の中で行われる消去法や分解は第5章で学ぶ。

本章では、手計算で分かる解と照合し、さらに代入して残差を確認した。真の解が分からない場合にも残差は計算できるが、残差が小さいことと解そのものの誤差が小さいことの関係は、問題によって変わる。この区別を第4章の条件数の学習につなげる。

3.10. 後の章のコードを読むために#

第5章以降のコードでは、本章までに紹介していない書き方もいくつか使う。すべてを今すぐ覚える必要はないが、「見たことがある」状態にしておくと、アルゴリズムの本体に集中しやすい。分からない記号が出てきたら、この表に戻ろう。

書き方	意味	主に使う章
`f(x) = cos(x) - x`	1行で関数を定義する。`function f(x) ... end` と同じ	第7章以降
`x -> x^3`	名前を付けない関数（無名関数）。関数を引数として渡すときに使う	第8章以降
`function g(a, b; tol=1e-10)`	`;` より後ろはキーワード引数で、既定値を持つ。呼び出すときは `g(1, 2; tol=1e-6)` と名前で指定する	第7章以降
`while 条件 ... end`	条件が成り立つ間、繰り返す。反復回数が事前に分からない計算で使う	第7章以降
`条件 ? A : B`	条件が真ならA、偽ならBを値とする	第1章以降
`条件 \|\| error("...")`	条件が偽のときだけ右側を実行する。入力の確認によく使う	第1章以降
`条件 && return x`	条件が真のときだけ右側を実行する	第7章以降
`(root=x, iterations=k)`	名前付きタプル。複数の結果をまとめて返し、`result.root` のように取り出す	第7章以降
`v = Float64[]`、`push!(v, x)`	空の配列を作り、末尾へ値を追加する。反復の履歴を記録するときに使う	第4章以降
`"x = $(x)"`	文字列の中へ値を埋め込む。`*` で文字列同士をつなぐこともできる	第1章以降
`range(0, 1, length=5)`、`collect(...)`	等間隔の点を作る。`collect` で通常の配列に変換する	第8章以降
`let ... end`	中で使う変数の範囲を区切る	第1章以降

コードには λ、≤、≈、α のような文字も現れる。Juliaの対話画面、Jupyter、VS Codeなどでは、\lambda、\le、\approx、\alpha と入力してTabキーを押すと変換できる。a ≤ b は a <= b、a ≈ b は isapprox(a, b) と同じ意味なので、入力しにくければこちらを使ってもよい。

次のコードは、これらの書き方をまとめて使った小さな例である。値を半分にし続け、tol 以下になるまでの回数と履歴を返す。実行前に、halve_until(1.0) が何回で止まるか予想してみよう。

function halve_until(x; tol=0.1)
    x > 0 || error("正の値を入力してください")
    history = Float64[]          # 空の配列を用意する
    while x > tol
        push!(history, x)        # 更新前の値を記録する
        x = x/2
    end
    return (value=x, steps=length(history), history=history)
end

result = halve_until(1.0)
println("最終値 = ", result.value, ", 回数 = ", result.steps)
println("途中の値 = ", result.history)
println("tol=0.01 の回数 = ", halve_until(1.0; tol=0.01).steps)

square = x -> x^2
println("無名関数 = ", square(3.0), ", 各成分へ適用 = ", square.([1.0, 2.0, 3.0]))
println("条件演算子 = ", result.steps > 3 ? "4回以上" : "3回以下")
println("0.1 + 0.2 ≈ 0.3 ? ", 0.1 + 0.2 ≈ 0.3)

最終値 = 0.0625, 回数 = 4
途中の値 = [1.0, 0.5, 0.25, 0.125]
tol=0.01 の回数 = 7
無名関数 = 9.0, 各成分へ適用 = [1.0, 4.0, 9.0]
条件演算子 = 4回以上
0.1 + 0.2 ≈ 0.3 ? true

$1\to0.5\to0.25\to0.125\to0.0625$ と進み、4回目の更新後に初めて0.1以下となるので、回数は4である。キーワード引数 tol を変えると、関数の中身を書き換えずに停止条件だけを変えられる。第7章の二分法やNewton法の関数も、同じ形で「許容誤差」「最大反復回数」を受け取り、結果と履歴を名前付きタプルで返す。

3.11. Python/NumPyとの対応を振り返る#

ここまでに使った操作を、Pythonを使ったことがある人向けに整理する。Python側では import numpy as np を実行しているものとする。Pythonの通常のリストとNumPy配列では演算の意味が異なるため、この表は数値計算用のNumPy配列との比較である。Pythonを使ったことがない場合は、Juliaの欄を復習に使えばよい。

操作	Python/NumPy	Julia
ベクトルを作る	`np.array([3.0, 4.0])`	`[3.0, 4.0]`
最初の成分	`x[0]`	`x[1]`
2行2列の行列	`np.array([[1., 2.], [0., 1.]])`	`[1.0 2.0; 0.0 1.0]`
行列の形	`A.shape`	`size(A)`
行列積	`A @ B`	`A * B`
要素ごとの積	`A * B`	`A .* B`
各成分の二乗	`x ** 2`	`x .^ 2`
各成分に指数関数を適用	`np.exp(x)`	`exp.(x)`
転置	`A.T`	`transpose(A)`
ベクトルの2ノルム	`np.linalg.norm(x)`	`norm(x)`
連立方程式を解く	`np.linalg.solve(A, b)`	`A \ b`
独立した数値配列を作る	`x.copy()`	`copy(x)`

Juliaの norm や dot には using LinearAlgebra が必要である。実数行列では A' も転置になるが、複素行列では共役転置を表す。複素数を扱う章ではこの違いも確認する。

3.12. つまずいたときに確認すること#

コードで思った結果が出ないときは、まず小さな入力で途中の値を表示し、手計算と比べよう。

状況	確認すること
`UndefVarError`：名前が定義されていない	変数・関数を定義するセルを実行したか。名前の綴りが一致するか。`norm` なら `using LinearAlgebra` を実行したか
`BoundsError`：配列の範囲外	`length(x)` や `size(A)` を確認する。この章の配列は添字1から始まる
`DimensionMismatch`：大きさが合わない	行列積の内側の大きさが一致しているか。行列積と要素ごとの演算を取り違えていないか
コードを変更したのに出力が同じ	ファイルなら保存してから再実行したか。Notebookなら変更箇所以降を実行し直したか

ライブラリを呼び出せたことと、意図した計算ができたことは別に確認する。特に、計算する前に結果の大きさやおおよその値を予想しておくと、取り違えに気づきやすい。

3.13. $\clubsuit$ 計算時間が気になったら#

計算が大きくなったときには、同じ結果をより少ない時間やメモリで得るための工夫が必要になる。本章では、まず式とコードの対応、入力と出力、結果の検証を優先する。以下は、その先で計算時間を調べたくなったときの補足である。

Juliaでは、関数を新しい引数の型で初めて呼び出す際などに、実行用コードを生成するコンパイルの時間がかかることがある。そのため、初回の時間を毎回の計算時間とみなさず、同じ条件で複数回調べる。例えば、Julia本体の @time で時間やメモリ確保の情報を表示できる。

@time euler_decay(0.4, 1.0, 5.0, 10_000)
@time euler_decay(0.4, 1.0, 5.0, 10_000)

10_000 は読みやすく区切った整数10000である。この章を順に実行していれば、同じ引数の型でこの関数をすでに使っているため、ここで必ず初回だけ遅くなるわけではない。また、短い計算の時間は揺らぎやすく、一度の測定で速さを断定しない。

処理を関数にまとめて必要な値を引数で渡すことや、数値配列の要素の型をそろえることは、コンパイラが計算内容を判断する助けにもなる。ただし、型の注釈を付けるだけで自動的に高速になるわけではない。測定で時間のかかる箇所を確かめ、変更後も同じ計算結果が得られることを確認しよう。

より詳しい測定に用いる BenchmarkTools は追加パッケージであり、本章の実行には不要である。詳しい考え方はJulia公式のPerformance Tipsを参照できる。高速化の細かな技法は、本章の演習では要求しない。

3.14. まとめ#

本章のまとめ

変数への代入、関数、繰り返しを使うと、数式の計算手順をコードで表せる。
ベクトルと行列では、成分と大きさを確かめる。行列積と要素ごとの演算は意味が異なる。
ベクトルの2ノルムは長さを、差の2ノルムは距離を測る。
連立方程式は A \ b で解き、手計算や残差との照合によって結果を確かめる。

次の第4章「誤差・条件数・計算量」では、ここで用いた誤差や残差の大きさを、問題の尺度や入力への敏感さと結びつけて考える。ノルムの使い分けや、必要な演算回数・メモリ量の見積もりも学ぶ。

3.15. 演習#

本文のコードを見ながら取り組んでよい。実行前の予想と、実行結果から分かったことを短く書こう。新しい構文を調べたり、計算時間を競ったりする必要はない。

演習

入力を変えると何が変わるか：euler_decay(0.4, 1.0, 5.0, 10) の初期量だけを2.0に変えると、最終値はどうなると予想できるか。更新式を見て理由を考え、実行して確認せよ。
積の意味を確かめる：本文の $\mathbf{B},\mathbf{C}$ について、B*C と B .* C の第1行第1列を手計算せよ。両者の値が違う理由を説明し、コードの出力と照合せよ。
答えを代入してみる：本文の連立方程式について、候補 $(1,2)^\mathsf{T}$ は解だろうか。まず二つの式に代入して考え、次にコードで残差を確認せよ。どの方程式をどの程度満たさないか、残差の成分を使って説明せよ。

3.16. 関連資料#

AI・機械学習論1：線形代数：ベクトル・行列・線形写像などの数学的な背景。
Julia公式：Getting Startedと配列のマニュアル：コードの実行方法、配列の作り方と演算。
Fundamentals of Numerical Computation：§2.2 Computing with matrices：数値計算での行列の扱い。