Juliaの LinearAlgebra.nullspace() 関数:初心者向け解説と応用

2025-05-27

零空間とは何か?

行列 mathbfA の零空間とは、mathbfAmathbfx=mathbf0 を満たす全てのベクトル mathbfx の集合のことです。ここで、mathbf0 はゼロベクトルを表します。言い換えれば、行列 mathbfA を左から掛けるとゼロベクトルになるようなベクトル mathbfx 全体の集まりが零空間です。

LinearAlgebra.nullspace() 関数の働き

LinearAlgebra.nullspace(A) 関数に正方行列または非正方行列 A を入力すると、この関数の出力は、A の零空間を張る**基底ベクトル(きていベクトル、basis vectors)**を列ベクトルとして持つ行列です。

基底ベクトルとは?

零空間を構成する上で、線形独立であり、かつ零空間内の全てのベクトルをそれらの線形結合で表現できるベクトルの集まりのことです。LinearAlgebra.nullspace() 関数は、通常、これらの基底ベクトルを正規直交化(せいきちょっこうか、orthonormalization)したものを返します。

関数の使い方と出力例

using LinearAlgebra

# 例として、以下の行列を定義します
A = [1 2 3; 2 4 6; 3 6 9]

# nullspace() 関数を使って零空間の基底を計算します
N = nullspace(A)

# 結果を表示します
println(N)

上記の例では、行列 A の各行は線形従属であるため、零空間は非自明なベクトルを含みます。nullspace(A) を実行すると、零空間の基底ベクトルが列ベクトルとして格納された行列 N が出力されます。この例の場合、出力される行列 N の列ベクトルは、mathbfAmathbfx=mathbf0 を満たす線形独立なベクトルを表しています。

  • 線形代数の理論的な理解
    零空間を計算することで、行列のランク(rank)や退化性(degeneracy)といった重要な性質を理解するのに役立ちます。
  • 特異値分解との関連
    零空間は、行列の特異値分解(SVD)とも深く関連しています。特に、特異値がゼロに対応する特異ベクトルの空間が零空間となります。
  • 線形方程式系の解の分析
    mathbfAmathbfx=mathbfb という線形方程式系において、mathbfb=mathbf0 の場合の非自明な解(自明な解 mathbfx=mathbf0 以外の解)を見つけるために使用できます。


よくあるエラーとトラブルシューティング

  1. 引数の型のエラー (MethodError)

    • エラー内容
      LinearAlgebra.nullspace() 関数に、行列(AbstractMatrix のサブタイプ)ではない型の引数を渡すと発生します。例えば、ベクトルやスカラー値を渡した場合などです。
    • トラブルシューティング
      • 引数が正しい行列の型であることを確認してください。typeof(your_variable) で変数の型を調べることができます。
      • データを適切に二次元配列(行列)として構成しているか確認してください。必要であれば、reshape() 関数などを使用して形状を変換します。

    <!-- end list -->

    using LinearAlgebra
    
    v = [1, 2, 3]
    # nullspace(v) # これはエラーになります
    
    A = [1 2; 3 4]
    N = nullspace(A) # これは正常に動作します
    
    • 警告内容
      内部的に nullspace() 関数は特異値分解(SVD)を利用しています。非常に病的な行列(条件数が大きいなど)の場合、SVD のアルゴリズムが指定された許容誤差内で収束しないことがあります。この場合、警告メッセージが表示されることがあります。
    • トラブルシューティング
      • 多くの場合、この警告が出ても結果は返されますが、精度が低い可能性があります。
      • 許容誤差を調整する必要があるかもしれません(tol キーワード引数)。ただし、通常はデフォルト値で十分です。
      • 行列の状態を調べ、数値的に不安定でないか確認してください。
      • 問題の性質によっては、他の数値計算手法を検討する必要があるかもしれません。
    using LinearAlgebra
    
    A = [1.0 1.0000000000000001; 1.0 1.0] # 病的な行列の例
    N = nullspace(A) # 警告が表示される可能性があります
    
    N_tol = nullspace(A, tol=1e-8) # 許容誤差を調整
    
  2. 零空間が空の場合(出力が空の行列に近い)

    • 状況
      与えられた行列がフルランク(列フルランクまたは行フルランク)である場合、零空間は自明な解(ゼロベクトルのみ)を持ちます。この場合、nullspace() 関数は、非常に小さな値(数値的ゼロに近い)を持つ空の行列または非常に小さい次元の行列を返すことがあります。
    • トラブルシューティング
      • 行列のランクを確認してください。rank(A) 関数でランクを計算できます。もし rank(A) が行列の列数と等しい(列フルランク)場合、零空間は自明であり、意味のある基底ベクトルは存在しません。
      • 期待される零空間の次元が本当にゼロであるか、それとも計算上の誤差で小さくなっているのかを検討してください。
    using LinearAlgebra
    
    A = [1 0; 0 1] # フルランクの正方行列
    N = nullspace(A)
    println(N) # ほぼ空の行列が出力されることがあります
    
    rank_A = rank(A)
    cols_A = size(A, 2)
    if rank_A == cols_A
        println("行列は列フルランクであり、零空間は自明です。")
    end
    
  3. 計算時間の問題

    • 状況
      非常に大きなサイズの行列に対して nullspace() 関数を実行すると、特異値分解の計算に時間がかかることがあります。
    • トラブルシューティング
      • 行列のサイズを削減できるか検討してください。
      • 他の線形代数ライブラリや、大規模データ向けの特殊なアルゴリズムの利用を検討するかもしれません。
      • 計算資源(CPU、メモリ)が十分であるか確認してください。
  4. 複素数行列の扱い

    • 状況
      入力行列が複素数要素を持つ場合、nullspace() 関数は複素数ベクトルを基底として持つ零空間を返します。
    • トラブルシューティング
      • 複素数演算を意図しているか確認してください。実数のみを扱う場合は、入力データに複素数が混入していないか確認します。
      • 後続の処理で複素数ベクトルを適切に扱えるように注意してください。

一般的なトラブルシューティングのヒント

  • ドキュメントを参照する
    Julia の公式ドキュメントや ?LinearAlgebra.nullspace で関数の詳細な仕様や引数を確認しましょう。
  • 簡単な例で試す
    問題のある行列の代わりに、既知の零空間を持つ小さな行列で nullspace() 関数を試してみることで、関数の基本的な動作を確認できます。
  • 入力データの確認
    size(), typeof(), eltype() などを使って、入力行列の形状、型、要素の型を確認しましょう。
  • エラーメッセージを注意深く読む
    Julia のエラーメッセージは、問題の原因に関する重要な情報を提供してくれます。


基本的な使い方

まず、最も基本的な使い方として、与えられた行列の零空間の基底を計算する例です。

using LinearAlgebra

# 例1: 線形従属な行を持つ行列
A1 = [1 2 3;
      2 4 6;
      0 0 0]

N1 = nullspace(A1)
println("行列 A1 の零空間の基底:")
println(N1)
println("N1 の列ベクトルは A1 * x = 0 の解 (近似的に):")
println(A1 * N1) # ほぼゼロベクトルになるはず

# 例2: 正方行列で特異な行列 (行列式がゼロ)
A2 = [1 1;
      1 1]

N2 = nullspace(A2)
println("\n行列 A2 の零空間の基底:")
println(N2)
println("N2 の列ベクトルは A2 * x = 0 の解 (近似的に):")
println(A2 * N2) # ほぼゼロベクトルになるはず

この例では、nullspace() 関数に与えられた行列の零空間を張る正規直交基底ベクトルが列ベクトルとして返されます。A * N がほぼゼロ行列になることで、得られた基底ベクトルが mathbfAmathbfx=mathbf0 の解空間を張っていることが確認できます。

許容誤差 (tol) の指定

nullspace() 関数には、計算の際の許容誤差を指定する tol キーワード引数があります。これは、特異値分解に基づいて零空間を計算する際に、どの程度の小さな特異値をゼロとみなすかを制御するために使用します。

using LinearAlgebra

# 非常に小さな値を持つ行列
A3 = [1.0 0.0000000001;
      0.0000000001 1.0]

N3_default = nullspace(A3)
println("\n行列 A3 の零空間 (デフォルトの許容誤差):")
println(N3_default)

N3_tight = nullspace(A3, tol=1e-10)
println("\n行列 A3 の零空間 (より小さな許容誤差):")
println(N3_tight)

# ランク落ちがわずかな行列
A4 = [1.0 1.0;
      1.0 1.00000000001]

N4_default = nullspace(A4)
println("\n行列 A4 の零空間 (デフォルトの許容誤差):")
println(N4_default)

N4_loose = nullspace(A4, tol=1e-5)
println("\n行列 A4 の零空間 (より大きな許容誤差):")
println(N4_loose)

この例では、tol の値を変更することで、零空間の次元や基底ベクトルがどのように影響を受けるかを示しています。許容誤差を小さくすると、より厳密にゼロに近い特異値のみが無視されるため、零空間の次元が小さくなる可能性があります。逆に、大きくすると、より多くの特異値がゼロとみなされ、零空間の次元が大きくなる可能性があります。

線形方程式系の解との関連

零空間は、非同次線形方程式系 mathbfAmathbfx=mathbfb の一般解を理解する上でも重要です。もし mathbfx_p が一つの特定解であれば、一般解は mathbfx_p+mathbfx_h と表されます。ここで、mathbfx_h は同次方程式 mathbfAmathbfx=mathbf0 の一般解であり、mathbfA の零空間の任意の線形結合で表されます。

using LinearAlgebra

# 例: 解を持つ線形方程式系
A5 = [1 2; 2 4]
b5 = [3; 6]

# 特定解を見つける (ここでは擬似逆行列を使用)
xp = pinv(A5) * b5
println("\n線形方程式 A5 * x = b5 の特定解:")
println(xp)

# 同次方程式 A5 * x = 0 の一般解 (零空間)
Nh5 = nullspace(A5)
println("\n同次方程式 A5 * x = 0 の一般解 (零空間の基底):")
println(Nh5)

# 一般解の例: 特定解 + 零空間の基底の線形結合
c = 2.5 # 任意の係数
xh = Nh5 * [c]
x_general = xp + xh
println("\n一般解の例:")
println(x_general)
println("A5 * x_general (近似的に b5):")
println(A5 * x_general) # ほぼ b5 になるはず

この例では、線形方程式系の特定解と零空間の基底を組み合わせることで、一般解を表現できることを示しています。

応用例: 特異値分解との関係

nullspace() 関数は内部で特異値分解(SVD)を利用しています。特異値分解の結果からも零空間の情報を得ることができます。

using LinearAlgebra

A6 = [1 0 0; 0 0 0; 0 0 2]

U, S, V = svd(A6)
println("\n行列 A6 の特異値分解:")
println("U =")
println(U)
println("S =")
println(S)
println("V =")
println(V)

# 特異値がゼロに近い列に対応する V の列ベクトルが零空間の基底
# (許容誤差に依存)
tolerance = 1e-8
zero_singular_values_indices = findall(s -> s < tolerance, S)
nullspace_basis_from_svd = V[:, zero_singular_values_indices]
println("\n特異値分解から求めた零空間の基底:")
println(nullspace_basis_from_svd)

nullspace_direct = nullspace(A6)
println("\nnullspace() 関数で直接求めた零空間の基底:")
println(nullspace_direct)

この例では、特異値分解の結果の特異値行列 S において、ゼロに近い特異値に対応する特異ベクトル(V の列ベクトル)が、nullspace() 関数が返す基底ベクトルと対応することを示唆しています。



特異値分解 (Singular Value Decomposition, SVD) を利用する方法

LinearAlgebra.nullspace() 関数は内部的に特異値分解を利用しています。したがって、直接 svd() 関数を使って特異値分解を行い、その結果から零空間の基底を求めることができます。

using LinearAlgebra

function nullspace_svd(A; tol=1e-8)
    U, S, V = svd(A)
    # 特異値が許容誤差よりも小さい列のインデックスを見つける
    zero_singular_value_indices = findall(s -> s < tol, S)
    # 対応する V の列ベクトルが零空間の基底
    return V[:, zero_singular_value_indices]
end

A = [1 2 3;
     2 4 6;
     0 0 0]

N_svd = nullspace_svd(A)
println("SVD を用いて計算した零空間の基底:")
println(N_svd)

N_builtin = nullspace(A)
println("\n組み込み関数 nullspace() で計算した零空間の基底:")
println(N_builtin)

# 結果が(符号などを除いて)一致することを確認
# all(abs.(N_svd) .≈ abs.(N_builtin))

この方法では、svd(A) で行列 A の特異値分解を行い、特異値ベクトル S の要素が許容誤差 tol より小さいインデックスに対応する右特異ベクトル(V の列ベクトル)を取り出すことで、零空間の基底を得ています。

利点

  • 許容誤差の制御をより細かく行えます。
  • 特異値分解の結果を他の目的にも利用できます。
  • nullspace() 関数の内部動作を理解するのに役立ちます。

ランク落ち検出と基本解の探索

行列のランクを計算し、ランク落ちがある場合に、ガウスの消去法などを応用して基本解を求める方法です。これは、特に手動で零空間を理解するのに役立ちますが、プログラミングで実装する場合はやや複雑になることがあります。

using LinearAlgebra

function nullspace_rank_reduction(A; tol=1e-8)
    m, n = size(A)
    rank_A = rank(A, tol=tol)

    if rank_A == n
        return Matrix{Float64}(undef, n, 0) # 零空間は自明
    end

    # より複雑なガウスの消去法やLU分解などを応用して基本解を求める必要がある
    # これは高度な内容になるため、ここでは概要のみを示す

    # (実際のコードは、ピボット選択、行の入れ替え、変数の依存関係の特定などを含む)
    # ...
    # return basis_vectors
    println("ランク落ち検出と基本解の探索は、より複雑な実装が必要です。")
    return nothing
end

A = [1 2; 2 4]
N_rank = nullspace_rank_reduction(A)
# println(N_rank)

この方法は、行列のランクが列数よりも小さい場合に零空間が存在することを利用します。ランク落ちの度合いに応じて、自由変数を特定し、それらに対する基本解を求める必要があります。

利点

  • 特定の構造を持つ行列に対して、より効率的なアルゴリズムを設計できる可能性があります。
  • 線形代数の基本的な概念を深く理解できます。

欠点

  • 数値的な安定性の問題に対処する必要がある場合があります。
  • 実装が複雑になることが多いです。

固有値と固有ベクトルを利用する方法 (特定のケース)

正方行列の場合、固有値がゼロに対応する固有ベクトルは、その行列の零空間のベクトルとなります。ただし、これは正方行列に限定され、かつ固有値が正確にゼロである場合に限ります。数値計算においては、固有値が非常に小さい場合に近似的な零空間ベクトルとして利用できることもあります。

using LinearAlgebra

function nullspace_eigen(A; tol=1e-8)
    if size(A, 1) != size(A, 2)
        error("固有値分解は正方行列に対してのみ定義されます。")
    end
    eigen_results = eigen(A)
    eigenvalues = eigen_results.values
    eigenvectors = eigen_results.vectors

    zero_eigenvalue_indices = findall(λ -> abs(λ) < tol, eigenvalues)
    return eigenvectors[:, zero_eigenvalue_indices]
end

A = [1 1; 1 1]

N_eigen = nullspace_eigen(A)
println("固有値分解を用いて計算した零空間の基底:")
println(N_eigen)

N_builtin = nullspace(A)
println("\n組み込み関数 nullspace() で計算した零空間の基底:")
println(N_builtin)

この方法では、eigen(A) で固有値と固有ベクトルを計算し、絶対値が許容誤差よりも小さい固有値に対応する固有ベクトルを零空間の基底としています。

利点

  • 固有値と固有ベクトルの概念との関連を理解できます。
  • 代数的重複度が 1 より大きいゼロ固有値の場合、複数の線形独立な固有ベクトルが存在しないことがあります(幾何学的重複度が代数的重複度より小さい場合)。
  • 数値的な誤差により、正確にゼロの固有値が得られない場合があります。
  • 正方行列にしか適用できません。