線形方程式の求解 の

線形方程式の求解について

高校で最初に代数を習ったときのことを覚えていますか？

たいていは、連立方程式を与えられ、xとyを求めなさいと言われたものでした。

4x - 3y = -2
2x + 7y = 16

式が2つ、未知数が2つ。少し代入すれば、すぐにx = 1, y = 2だとわかります。

前の例の左辺をもう一度見てみましょう。 これは、行列[4 -3; 2 7]に未知のベクトル[x, y]を掛けると、ベクトル[-2, 16]になるように見えます。

線形代数の表記では、A x = bです。これは、行列には大文字、ベクトルには小文字を使うという慣例に従っています。

A x = 0 と零空間

重要な特殊ケースは、すべての方程式の右辺がゼロのときです。

4x - 3y = 0
2x + 7y = 0

上の例では、唯一の解はx = y = 0のときだけですが、これは通常あまり面白くありません。

A x = 0の自明でない解は、Aの行列式がゼロのときにのみ存在します。

julia> using LinearAlgebra

julia> A = [4 -3; 2 7]
2×2 Matrix{Int64}:
 4  -3
 2   7

julia> det(A)  # not zero!
34.0

代わりに、次の方程式を考えてみましょう。

4x - 3y = 0
8x + 6y = 0

2番目の方程式は1番目の方程式を2倍しただけなので、新しい情報は何もありません。 x = 0.75yとなる値はどれも解になります。

これらの(x, y)の値は、2次元空間の直線上に並びます。 解は無限にあり、それらは行列の零空間を形作ります。

この場合、行列の行は線形独立ではなく、行列の階数は行数または列数より小さくなります。

julia> A = [4 -3; 8 -6]
2×2 Matrix{Int64}:
 4  -3
 8  -6

julia> det(A)
0.0

julia> size(A)  # total number of rows and columns
(2, 2)

julia> rank(A)  # number of linearly-independent rows (or columns)
1

julia> nullspace(A)
2×1 Matrix{Float64}:
 -0.6
 -0.8

julia> nullspace(A) |> norm  # Julia returns unit vectors when possible
1.0

nullspace()関数は、正規化されて単位ベクトルになった1つの点を返しますが、このベクトルのスカラー倍もすべて零空間に含まれます。

A x = b を解く

さて、1000個の未知数に対して1000本の方程式があるとしたらどうでしょうか？ それがばかげているように聞こえるなら、ワイヤレスのデジタルトランスデューサーが今や安価で多用途であることを思い出してください（ひずみ、風速、3軸加速度など、何でも測定できます）。 吊り橋のような現代的な構造物を1000個のトランスデューサーで監視するのは、まったく合理的です。 データストリームを解釈し、事態が危険になれば警報を発する準備ができた制御システムが必要です。

ここでもA x = bです。ここではA（エンジニアリング設計から）とb（トランスデューサーからの測定値）はわかっていますが、xを見つける必要があります。

これに慣れていない人が最初に思いつくのは、何とかしてAで「割って」右辺に移すことです。

実際、Aには逆行列があり（すべての正方行列ではありませんが、ほとんどの正方行列に存在します。そのためには行列式がゼロでない必要があります）、計算はうまくいきます（遅いです！）。

julia> A = [4 -3; 2 7]
2×2 Matrix{Int64}:
 4  -3
 2   7

julia> b = [-2, 16]
2-element Vector{Int64}:
 -2
 16

# the inverse matrix
julia> inv(A)
2×2 Matrix{Float64}:
  0.205882   0.0882353
 -0.0588235  0.117647

julia> inv(A) * b
2-element Vector{Float64}:
 0.9999999999999999
 2.0

残念ながら、逆行列の計算は小さい行列では遅く、大きい行列では途方もなく遅くなります。 計算中に橋が崩落してしまうとしたら、理想的ではありませんね！

幸いなことに、他のアルゴリズムは劇的に速いです（この場合はガウスの消去法です）。

Juliaは（Matlabをまねて）ソルバーにバックスラッシュを使うだけです（専門的には「左除算」です）。

julia> x = A \ b
2-element Vector{Float64}:
 1.0
 2.0

これは簡単な例ですが、注意すべき細かい点がいくつかあります。

• 係数は列をそろえる必要があるので、欠けている項には必要に応じてゼロを挿入します。
• 方程式（つまり行列の行）は線形独立である必要があります。3行目が1行目と2行目の和にすぎないなら、新しい情報を加えないので、問題は劣決定になります。

rank()コマンドは線形独立な行／列の数を返すので、これが未知変数の数と等しくなることを目指します。

julia> rank(A)
2

長方行列

10代になる前のころに、N個の未知数を解くにはN本の方程式が必要だと教わったのではないでしょうか。

数学のほとんどのことと同じように、現実はもう少し微妙です（前のセクションで触れた線形独立性の詳細だけではありません）。

行が「少なすぎる」場合

N-1本の方程式（行列Aの行）では、問題は劣決定です。 絶望して諦めるべきでしょうか？

状況によります！

N-1本の方程式にも、まだ多くの情報が含まれています。 幾何学的に言えば、N次元空間で解を_点_として定めるにはNが必要ですが、N-1でも解が_直線_上のどこかにあることはわかります。

すると解は無限にありますが、解が_ない_パラメータ空間も無限に存在することになります。

解の直線は、まともなエンジニアなら心配すべきパラメータ空間の領域を通っているでしょうか？ あるいは、構造物を一般の立ち入り禁止にすることを正当化するような領域でしょうか（吊り橋は強い横風で揺れやすくなります）？ それを調べるために、もう少し計算してみる価値はあるかもしれません！

julia> A3 = [4 -3 1; 2 2 2; 3 -4 -3]
3×3 Matrix{Int64}:
 4  -3   1
 2   2   2
 3  -4  -3

julia> b3 = [1, 12, -14]
3-element Vector{Int64}:
   1
  12
 -14

# fully-determined solution
julia> A3 \ b3
3-element Vector{Float64}:
 1.0
 2.0
 3.0

# remove row 3 from A3 and B3
julia> A2 = A3[1:2, :]
2×3 Matrix{Int64}:
 4  -3  1
 2   2  2

julia> b2 = b3[1:2]
2-element Vector{Int64}:
  1
 12

# an under-determined solution
julia> z1 = A2 \ b2
3-element Vector{Float64}:
 1.594594594594594
 2.445945945945947
 1.9594594594594592

バックスラッシュのソルバーが解を返してくれます！ 非公式には、これはどうやらノルムが最小の解のようです（幾何学的には原点に最も近い）。

他の解を得るには、A2のnullspaceが必要です。

元の解に零空間の任意のスカラー倍を加えると、別の有効な解が得られます。

# get the nummspace of A2
julia> N = nullspace(A2)
3×1 Matrix{Float64}:
 -0.46499055497527725
 -0.3487429162314578
  0.813733471206735

# add a random multiple of the nullspace to the earlier solution
julia> z = z1 + rand() * N
3×1 Matrix{Float64}:
 1.274331860650258
 2.205748895487695
 2.5199192438620472

# our random z is a valid solution, recovering b2 = [1, 12]
julia> A2 * z
2×1 Matrix{Float64}:
  0.9999999999999947
 12.0

同様に、N-2本の方程式は（無限にある）解を平面に制約し、同じ原理が当てはまります。

行が「多すぎる」場合

逆の状況は、未知数より_多く_の方程式があるのに、それでもなぜか線形独立である場合に起こります。

現実のエンジニアリングでは、これはまったく普通のことで、良いことだと見なされています！

トランスデューサーには精度の限界があり、bベクトルにも精度の限界があり、計算にはノイズが含まれます。 この場合、解はノイズのあるデータに対する最小二乗近似になります。

最も単純な手法は行列のpseudoinverseを使うもので、Juliaではpinv()関数として実装されています。

これは（非特異な）正方行列の逆行列とほぼ同じように使え、厳密な解ではなく変数の最小二乗推定値を与えます。

# create 5-row A and b for 2 variables
julia> A5 = [4 -3; 2 7; -1 2; 1 -1; 3 1]
5×2 Matrix{Int64}:
  4  -3
  2   7
 -1   2
  1  -1
  3   1

julia> rank(A5)
2

# add some random-normal noise to b5
julia> b5n = [-2, 16, 3, -1, 5] + randn(5) * 0.1
5-element Vector{Float64}:
 -1.9337349223581917
 16.024913008059457
  3.0087646177373992
 -0.8675748721176717
  4.914568047308805

# solve to get x ≈ 1, y ≈ 2, using the pseudoinverse
julia> pinv(A5) * b5n
2-element Vector{Float64}:
 1.006117942769543
 1.9962973764508272

固有値と固有ベクトル

前のセクションでは、A x = bの解について説明しました。これは、おなじみの代数と等価です。

このセクションでは、A x = λ xの解を扱います。ここでλはスカラーです。

これは線形代数に特有の概念ですが、幾何学的に解釈できます。

適切なxとλの値に対して、正方行列Aは非ゼロベクトルxの長さをλ倍に伸縮させ、その向きは変えません（負のλが向きを反転させる場合を除く）。

ニッチに聞こえますが、実際には_ばかげているほど_役に立ちます！

用語: λの有効な値はAの_固有値_で、対応するxの値はAの_固有ベクトル_です。 残念ながら、途中で言語が切り替わる言葉と付き合っていくしかありません。

学生は通常、2×2行列の固有値／固有ベクトルを手計算する方法を教わりますが、コンピューターを使うとはるかに簡単です（ただしそれでもかなり遅く、一般的な場合、n×n行列ではO(n^3)でスケールします）。

julia> A = rand(-9:9, 2, 2)
2×2 Matrix{Int64}:
 9  5
 3  2

julia> F = eigen(A)
Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
2-element Vector{Float64}:
  0.27984674554472466
 10.720153254455274
vectors:
2×2 Matrix{Float64}:
 -0.497417  0.945605
  0.867511  0.325317

# eigenvalues
julia> F.values
2-element Vector{Float64}:
  0.27984674554472466
 10.720153254455274

# each column is an eigenvector, normalized to a unit vector
julia> F.vectors
2×2 Matrix{Float64}:
 -0.497417  0.945605
  0.867511  0.325317

一般に、n×n行列にはn個の固有値がありますが、常に_異なる_値とは限りません。 それらはn次多項式（特性多項式と呼ばれます）の根だと考えてください。根は重複することもあり、実数値行列であってもしばしば複素数になります。

各固有ベクトルは1つの方向を表し、その任意のスカラー倍も有効な固有ベクトルです。 以降の計算を便利にするため、Juliaはノルムが1の単位ベクトルを返します。

応用

固有ベクトルがなぜ重要かを説明するのは、数段落ではなく、理想的には500ページの教科書の仕事です。 この話題は、現代の応用数学の_非常に多く_に浸透しています。

大まかに言えば、固有ベクトルはデータセット内で「最も重要」な軸（方向）を表します。 固有値は各軸の「相対的な重要度」を示します（入力の正規化に関するいくつかの仮定の下で）。

これが実際に_何を意味するか_は、応用によって異なります。

主成分分析

典型的なデータサイエンスの問題では、データの列として保存された100以上の「特徴量」があるかもしれません。 ほぼ避けられないことに、ノイズ、冗長性、不要な相関が存在します。

次元削減を行う必要がありますが、PCAは混沌に秩序をもたらす1つの方法です：

• データの共分散行列を計算します。
• これで正方行列が得られたので、次に固有値と固有ベクトルを計算します。
• 固有ベクトルを固有値の降順（より正確には絶対値|λ|）に並べ替えます。
• 上位k個の固有ベクトルがprincipal componentsになります。ここでkは元の特徴量の数よりかなり小さくなります。
• データセット全体をこれらのk軸に射影し、興味深いパターンを探し始めます。

この形のPCAは、創薬で分子特性を改善しようとする医学研究者から、世論調査データから有権者の選好を理解しようとする政治運動家まで、実に多様なグループに使われてきました。 おそらく、何かを売り込もうとするマーケティンググループにも使われているでしょう。しかし、どんな技術も善にも悪にも使えます。

PCAは最小構成のJuliaには含まれていませんが、（Exercismの外では）MultivariateStatsパッケージに必要なものが入っています。

画像処理

PCAのアイデアをさらに進めると、デジタル画像はピクセル値の行列にすぎず、その主成分を計算できます。

これは何十年もの間、画像圧縮に使われてきました。そこでは、ファイルサイズを小さくするときに何を保持するのが最も重要かをPCAが教えてくれます。

ますます、PCAは顔認識などの画像分類器の不可欠な部分になっています。 次に空港を歩くとき、ビッグ・ブラザーはただ見ているだけでなく、見たものを理解するために線形代数も使っています！

機械工学

それほど物議を醸さない例では、機械部品の主軸は、その慣性モーメントテンソル I（少し違う名前の行列）の固有ベクトルです。

車の各ホイールには、このテンソルIの非対角要素をゼロにするように調整されたバランスウェイトが1つ以上あるでしょう。 町の整備士は間違いなく計算を避けています（航空機の設計者やロケットエンジニアとは対照的です）が、「ぶれ」はテンソルのクロス項を表す日常的な言葉にすぎません。この場合、クロス項があると乗り心地が悪くなり、ベアリングの機械的摩耗が増えます。