【Julia】LinearAlgebra.svd()のエラー解決から高速化まで

Julia言語のLinearAlgebra.svd()関数は、行列の特異値分解 (Singular Value Decomposition, SVD) を計算するために使われます。SVDは、線形代数における非常に強力なツールであり、データ分析、画像処理、機械学習など、幅広い分野で利用されています。

特異値分解 (SVD) とは？

任意の実数または複素数の行列 A (mtimesn 行列) は、以下のように分解することができます。

$A = U S V^\*$

ここで、各要素は以下の通りです。

• V∗: ntimesn のユニタリー行列 V の共役転置 (または転置)。V の列ベクトルは、A∗A の固有ベクトルであり、「右特異ベクトル」と呼ばれます。
• S: mtimesn の対角行列。対角成分には、非負の実数である「特異値 (singular values)」が降順に並んでいます。特異値は、A の「スケール」または「重要度」を表します。対角成分以外はすべて0です。
• U: mtimesm のユニタリー行列 (または直交行列)。U の列ベクトルは、$A A^\*$ の固有ベクトルであり、「左特異ベクトル」と呼ばれます。

簡単に言うと、SVDは任意の行列を、回転 (U)、スケーリング (S)、*別の回転 (V)** の3つの操作の組み合わせとして表現できることを示しています。

LinearAlgebra.svd() の使い方と返り値

JuliaでLinearAlgebra.svd()関数を行列 A に適用すると、SVD型のオブジェクトが返されます。このオブジェクトには、特異値分解の結果である U, S, V (または V の転置) が含まれています。

基本的な使い方は以下の通りです。

using LinearAlgebra

A = [1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 9.0] # 例として3x3行列を定義

F = svd(A)

F は SVD オブジェクトなので、以下のフィールドにアクセスして各成分を取り出すことができます。

• F.V: 右特異ベクトルを行列として取得します。
• F.S: 特異値をベクトルとして取得します（対角行列ではありません）。
• F.U: 左特異ベクトルを行列として取得します。

例

using LinearAlgebra

A = [1.0 2.0; 3.0 4.0]

F = svd(A)

println("U行列:")
println(F.U)

println("\n特異値ベクトル:")
println(F.S)

println("\nV行列:")
println(F.V)

# 分解された成分を使って元の行列を再構築する (近似)
# 対角行列Sを作るには、F.SベクトルをDiagonal()でラップします
A_reconstructed = F.U * Diagonal(F.S) * F.V'
println("\n再構築されたA行列:")
println(A_reconstructed)

svd() のオプション

svd() 関数には、いくつかのオプション引数があります。

• alg: SVDを計算するためのアルゴリズムを指定できます。通常はデフォルトで問題ありませんが、特定のケース（例えば、非常に大規模な行列や特定の特性を持つ行列）でパフォーマンスを最適化するために変更することがあります。利用可能なアルゴリズムは LinearAlgebra.DivideAndConquer() や LinearAlgebra.QRIteration() などです。

  F = svd(A, alg=LinearAlgebra.DivideAndConquer())
  

• full::Bool:

  • true (デフォルト): U は mtimesm 行列、 V は ntimesn 行列として返されます。つまり、完全な（「full」）ユニタリー行列が計算されます。
  • false: 「コンパクトSVD」または「薄いSVD (thin SVD)」が計算されます。この場合、U は mtimesmin(m,n) 行列、V は ntimesmin(m,n) 行列として返されます。これは、特に一方の次元がもう一方の次元よりもはるかに大きい場合（例えば、非常に「細長い」行列）に、計算コストとメモリ使用量を削減するために役立ちます。

  <!-- end list -->

  A = randn(5, 3) # 5x3の行列
  F_full = svd(A, full=true) # Uは5x5, Vは3x3
  F_thin = svd(A, full=false) # Uは5x3, Vは3x3
  

特異値分解は、以下のような様々な応用があります。

• 擬似逆行列の計算 (Pseudoinverse Calculation): 逆行列が存在しない行列に対しても、SVDを使って擬似逆行列を計算できます。これは線形方程式の最小二乗解を求める際などに役立ちます。
• 推薦システム: ユーザーとアイテムの評価行列に対してSVDを適用することで、ユーザーの好みを予測したり、アイテムを推薦したりできます。
• 画像圧縮: 画像データをSVDで分解し、特異値の大きい上位の成分のみを保持することで、画像を圧縮できます。
• ノイズ除去: 特異値の小さい成分（通常、ノイズに相当する）を無視することで、データからノイズを除去できます。
• 次元削減 (Dimensionality Reduction): 主成分分析 (PCA) の基盤として利用され、データの主要な変動を捉えることができます。特異値の大きい成分ほど重要度が高いと解釈できます。

---

LAPACKException

SVDの計算は、内部的にLAPACK (Linear Algebra PACKage) という高性能な線形代数ライブラリに依存しています。LAPACKException は、LAPACKルーチンがエラーを報告したときに発生します。これにはいくつかの原因が考えられます。

考えられる原因とトラブルシューティング

• NaN または Inf を含む行列: 行列に非数 (NaN) や無限大 (Inf) の値が含まれている場合、SVDはエラーを発生させます。

  • トラブルシューティング:
    • 行列の要素をチェックする: any(isnan, A) や any(isinf, A) を使って、行列に問題のある値が含まれていないか確認します。
    • 前処理を行う: データの前処理ステップで、NaN や Inf を適切な値（例えば、ゼロ、平均値、中央値など）に置き換えるか、該当する行/列を削除することを検討してください。

• 特異な（Singular）またはほぼ特異な（Nearly Singular）行列: 行列が特異（逆行列が存在しない）であるか、非常に小さい特異値を持つ（ほぼ特異）場合、数値的な問題が発生しやすくなります。

  • トラブルシューティング:
    • 行列のランクを確認する: rank(A) を使って行列のランクを調べ、期待通りのランクがあるか確認します。
    • 小さな摂動を加える: ごく稀に、行列のごくわずかな要素を修正する（例: A + 1e-15 * I のように単位行列に微小な値を加える）ことで、数値的な安定性が向上し、エラーが回避されることがあります。ただし、これは問題の根本的な解決にはならず、結果に影響を与える可能性があるため注意が必要です。
    • SVDのアルゴリズムを切り替える: svd() 関数はデフォルトで gesdd! を使用しますが、場合によっては gesvd! の方が安定していることがあります。 <!-- end list -->

    using LinearAlgebra
    A = rand(3, 3) # 問題のある行列を想定
    try
        F = svd(A, alg=LinearAlgebra.DivideAndConquer()) # デフォルトのgesdd!
    catch e
        if isa(e, LAPACKException)
            println("LAPACKExceptionが発生しました。別のアルゴリズムを試します。")
            F = svd(A, alg=LinearAlgebra.QRIteration()) # gesvd!に対応する
        else
            rethrow(e)
        end
    end
    

    これは、特に大規模なSVDを多数実行する際に、特定の病的な行列で発生することが報告されています。

メモリ不足エラー (OutOfMemoryError)

非常に大きな行列のSVDを計算しようとすると、メモリが足りなくなることがあります。特に、full=true (デフォルト) の場合、大きな U および V 行列が生成されるため、メモリ消費量が大きくなります。

考えられる原因とトラブルシューティング

• 行列のサイズ: 行列のサイズが非常に大きい場合、SVDの計算に必要なメモリがシステムの物理メモリやスワップ領域を超過することがあります。
  • トラブルシューティング:
    • コンパクトSVD (full=false) の使用: もし、すべての特異ベクトルが必要ではなく、最も重要な特異値とそのベクトル（コンパクトSVD）だけで十分な場合は、full=false オプションを使用します。これにより、生成される U と V の次元が min(m,n) に削減され、大幅にメモリを節約できます。 <!-- end list -->

    A = rand(10000, 1000) # 大規模な行列
    F_thin = svd(A, full=false) # メモリ消費を抑える
    

    これは、特に「細長い」行列 (m >> n または n >> m) で非常に有効です。
    • 近似SVD/部分的SVDライブラリの利用: 全ての特異値や特異ベクトルが必要ない場合（例えば、上位数個の特異値とそのベクトルで十分な場合）、近似SVDを計算するライブラリを検討します。これにより、計算時間とメモリ使用量を劇的に削減できます。
      • TSVD.jl: トランケートSVD (Truncated SVD) を計算します。
      • KrylovKit.jl: 厳密なSVDではないですが、特定の大きな固有値/特異値と対応するベクトルを見つけるのに適しています。
      • IncrementalSVD.jl: データがストリーミングされる場合や、メモリに収まらない場合に有用な増分SVDを提供します。
    • アウトオブコアSVD: 行列全体をメモリにロードできない場合は、アウトオブコアSVD (Out-of-core SVD) を検討します。これは、行列をチャンクに分割し、ディスクからデータを読み込みながらSVDを計算する手法です。Juliaには直接的な標準ライブラリのサポートはありませんが、Mmap.jl などと組み合わせて自分で実装したり、外部ライブラリを探したりする必要があります。
    • データ型の変更: Float64 の代わりに Float32 を使用するなど、よりメモリ効率の良いデータ型を使用することで、メモリ消費を抑えることができます。ただし、これには精度のトレードオフが伴います。

MethodError またはサポートされていない型

svd() 関数は、すべての数値型や抽象行列型に対して直接サポートされているわけではありません。

考えられる原因とトラブルシューティング

• 特殊な行列型: Diagonal, Symmetric, Hermitian など、特殊な構造を持つ行列型は、svd() が内部的に最適化されたメソッドを使用することがあります。しかし、独自に定義した抽象行列型や、特定のパッケージで提供される特殊な行列型の場合、直接 svd() が適用できないことがあります。

  • トラブルシューティング:
    • 具体的な行列型に変換: Matrix(A) のように、その行列を標準の Matrix 型に変換してから svd() を適用することを検討します。ただし、これにより特殊な構造による最適化の恩恵が失われたり、大きな行列の場合はメモリ消費が増えたりする可能性があります。
    • ドキュメントの確認とパッケージの調査: 使用している特殊な行列型のドキュメントを確認し、SVDをサポートしているか、または推奨されるSVDの計算方法があるかを確認します。関連するパッケージがSVDを提供する場合があります。

• BigFloat やカスタム数値型: BigFloat のような高精度数値型や、LinearAlgebra が直接サポートしないカスタム数値型の場合、svd() が MethodError を発生させることがあります。

  • トラブルシューティング:
    • GenericLinearAlgebra.jl の使用: GenericLinearAlgebra.jl パッケージは、Juliaの標準ライブラリである LinearAlgebra を拡張し、より汎用的な数値型（例: BigFloat）に対する線形代数演算を提供します。 <!-- end list -->

    using LinearAlgebra
    using GenericLinearAlgebra # これをインポートすると、svdがBigFloatに対応する
    A_big = BigFloat.(rand(5, 5))
    F_big = svd(A_big) # GenericLinearAlgebra.jlがsvdを拡張して対応
    

数値精度に関する問題

SVDの計算は、浮動小数点数の精度に影響されることがあります。特に、条件数が非常に悪い行列（小さな特異値と大きな特異値の差が大きい行列）では、計算された特異値や特異ベクトルの精度が低下する可能性があります。

考えられる原因とトラブルシューティング

• 条件数の悪い行列: 特異値の比率（最大特異値 / 最小特異値）が大きい行列は「条件が悪い」とされ、数値的に不安定になりやすいです。
  • トラブルシューティング:
    • cond(A) で条件数を確認: cond(A) を使って行列の条件数を計算します。非常に大きな値（例: 1e15 以上）であれば、計算精度が問題になる可能性があります。
    • データのスケーリング: 行列の要素のスケールが大きく異なる場合、データの前処理としてスケーリングを検討します。これにより、条件数が改善されることがあります。
    • 高精度演算: BigFloat を使用して計算精度を上げることを検討します。ただし、計算コストは大幅に増加します。
    • 問題の定式化の見直し: 元の数学的な問題を、SVDの計算精度が影響しにくい形に再定式化できないか検討します。例えば、逆行列を直接計算する代わりに、線形方程式を解くためにSVDベースのアプローチ（擬似逆行列など）を使用する方が安定することがあります。

• MWE (Minimal Working Example) の作成: エラーが発生した場合、問題の行列を最小限のコードで再現できるMWEを作成し、それを使ってデバッグを行うのが効果的です。これにより、複雑なコードの中から問題の原因を切り分けやすくなります。
• エラーメッセージを読む: Juliaのエラーメッセージは非常に役立ちます。ERROR: の後に続くメッセージや、Stacktrace を注意深く読み、問題の根本原因を特定するためのヒントを探しましょう。

---

SVDに関する基本的な説明は以前に行いましたが、ここではより実践的なコード例に焦点を当てます。

基本的なSVDの計算と再構成

最も基本的な使い方は、与えられた行列のSVDを計算し、その結果を使って元の行列を再構成する例です。

using LinearAlgebra

# 例として行列Aを定義
A = [1.0 2.0 3.0;
     4.0 5.0 6.0;
     7.0 8.0 9.0;
     10.0 11.0 12.0] # 4x3の行列

println("元の行列 A:")
println(A)

# SVDを計算
F = svd(A)

# SVDオブジェクトからU, S, Vを取得
U = F.U       # 左特異ベクトル (4x4行列)
S = F.S       # 特異値ベクトル (長さmin(m,n)=3のベクトル)
V = F.V       # 右特異ベクトル (3x3行列)

println("\nU行列:")
println(U)

println("\n特異値ベクトル S:")
println(S)

println("\nV行列:")
println(V)

# 元の行列を再構成
# Sはベクトルなので、Diagonal()を使って対角行列に変換する必要がある
A_reconstructed = U * Diagonal(S) * V'

println("\n再構成された行列 A_reconstructed:")
println(A_reconstructed)

# 元の行列と再構成された行列がほぼ同じであることを確認
println("\n元の行列と再構成された行列の差のノルム:")
println(norm(A - A_reconstructed)) # 非常に小さい値になるはず

解説

• norm(A - A_reconstructed) は、元の行列と再構成された行列の差が非常に小さいことを示し、SVDが正しく機能していることを確認できます。
• V' は V の共役転置（実数の場合は単なる転置）です。
• S は特異値のベクトルであるため、行列積を行う際には Diagonal(S) を使って対角行列に変換する必要があります。
• F.U、F.S、F.V でそれぞれの成分にアクセスします。
• svd(A) は SVD 型のオブジェクト F を返します。

次元削減（主成分分析 - PCA）

SVDは主成分分析（PCA）の基礎としてよく使われます。データの分散の大部分を説明する少数の「主成分」を特定するために、上位の特異値と対応する特異ベクトルを使用します。

ここでは、データ行列 X (各行がサンプル、各列が特徴量) を中心化した後、SVDを使って主成分を抽出する例を示します。

using LinearAlgebra
using Statistics # mean関数を使うため

# 例としてデータ行列 X を定義 (各行が観測値、各列が特徴量)
# 10個のサンプル、3つの特徴量
X = [10.0 1.0 2.0;
     11.0 2.0 2.5;
     12.0 3.0 3.0;
     13.0 4.0 3.5;
     14.0 5.0 4.0;
     10.5 1.5 2.2;
     11.5 2.5 2.7;
     12.5 3.5 3.2;
     13.5 4.5 3.7;
     14.5 5.5 4.2]

println("元のデータ行列 X:")
println(X)

# データを中心化 (各列の平均を引く)
# PCAでは通常、データの平均を0にする
X_centered = X .- mean(X, dims=1)

println("\n中心化されたデータ行列 X_centered:")
println(X_centered)

# 中心化されたデータ行列のSVDを計算
F_pca = svd(X_centered)

# 左特異ベクトル U (主成分の方向)
# 特異値 S (各主成分の「重要度」)
# 右特異ベクトル V (ここでは直接使わないことが多いが、変換の基底)
U_pca = F_pca.U
S_pca = F_pca.S
V_pca = F_pca.V

println("\n特異値 (主成分の分散の大きさに対応):")
println(S_pca)

# 主成分の方向 (Principal Components directions) は V の列ベクトルに対応
# SVDの定義 A = U * S * V' で、A がデータ行列の場合、V の列ベクトルが主成分の方向
# データが (サンプル数 x 特徴量数) の場合、V は (特徴量数 x 特徴量数) の行列となり、
# その列ベクトルが主成分の方向を示す。
println("\n主成分の方向 (V の列ベクトル):")
println(V_pca)

# データを新しい主成分空間に射影 (次元削減)
# 例えば、上位2つの主成分のみを保持する場合
k = 2 # 保持する主成分の数

# 射影されたデータ (スコア)
# U行列の最初のk列と、特異値の最初のk個の積 (U * Diagonal(S)) を使って計算することもできますが、
# 実際には X_centered * V_pca[:, 1:k] で計算することが多いです。
# SVDによるPCAでは、U * Diagonal(S) が主成分スコアに相当します。
# あるいは、V_pcaの最初のk列を使ってデータを変換します。
# ここでは U * Diagonal(S) の形で計算します。
# S は既に降順に並んでいるので、最初の k 個を取るだけでよい。
projected_data = U_pca[:, 1:k] * Diagonal(S_pca[1:k])
# もしくは X_centered * V_pca[:, 1:k]
# projected_data = X_centered * V_pca[:, 1:k] # これは PCA のスコア (主成分の値) を与える
# 上記二つの結果は、符号が異なる場合がありますが、数学的には同等です。

println("\n上位 $(k) 個の主成分に射影されたデータ:")
println(projected_data)

# 各主成分が説明する分散の割合
explained_variance_ratio = (S_pca.^2) / sum(S_pca.^2)
println("\n各主成分が説明する分散の割合:")
println(explained_variance_ratio)

# 累積説明分散
cumulative_explained_variance = cumsum(explained_variance_ratio)
println("\n累積説明分散:")
println(cumulative_explained_variance)

解説

• U * Diagonal(S) あるいは X_centered * V の最初の k 列を取ることで、元のデータを k 次元の主成分空間に射影した結果（主成分スコア）が得られます。
• V 行列の列ベクトルは、主成分の方向（新しい基底）を表します。
• 特異値 S は、各主成分が説明する分散の大きさに比例します。大きいほど重要です。
• 中心化されたデータ行列に対してSVDを実行します。
• PCAでは、まずデータを行列として表し、各特徴量（列）の平均を引いて中心化します。

画像圧縮

SVDは、画像を低ランク近似することで圧縮するのに非常に効果的です。画像はピクセル値の行列として表現できるため、SVDを適用して、最も重要な情報（大きい特異値に対応する成分）のみを保持することで、画像を再構築します。

この例では、Images.jl パッケージを使用します。事前に Pkg.add("Images") と Pkg.add("TestImages") を実行してください。

using Images
using TestImages # 例の画像用
using LinearAlgebra
using Plots # 画像表示用 (必要であれば Pkg.add("Plots") も)

# グレースケール画像を読み込む
# "cameraman" はよく使われるテスト画像
img_gray = Gray.(testimage("cameraman"))
img_matrix = Float64.(channelview(img_gray)) # 行列に変換 (0.0-1.0の範囲)

println("元の画像のサイズ: $(size(img_matrix))")

# SVDを計算
F_img = svd(img_matrix)
U_img, S_img, V_img = F_img.U, F_img.S, F_img.V

# 異なるランクで画像を再構成し、圧縮効果を見る
ranks = [5, 20, 50, 100] # 使用する特異値の数 (ランク)

reconstructed_images = []

for k in ranks
    # 上位 k 個の特異値と対応する特異ベクトルを使用
    U_k = U_img[:, 1:k]
    S_k = Diagonal(S_img[1:k])
    V_k = V_img[:, 1:k]

    # 画像を再構成
    img_reconstructed = U_k * S_k * V_k'

    # 値が0.0-1.0の範囲に収まるようにクランプ (clamp01!)
    clamp01!(img_reconstructed)

    push!(reconstructed_images, Gray.(img_reconstructed))
    println("ランク $(k) で再構成された画像のサイズ: $(size(img_reconstructed))")
end

# 元の画像と再構成された画像を並べて表示
plot_list = [plot(img_gray, title="Original")]
for (i, img_rec) in enumerate(reconstructed_images)
    push!(plot_list, plot(img_rec, title="Rank $(ranks[i])"))
end

plot(plot_list..., layout=(1, length(plot_list)), size=(1200, 400), margin=5Plots.mm)

# 圧縮率の目安を計算 (元のデータ量に対する保存すべきデータ量の比率)
# 元の行列サイズ: m * n
# ランク k の場合: U (m*k) + S (k) + V (n*k) = k*(m+n+1)
m, n = size(img_matrix)
println("\n圧縮率の目安:")
for k in ranks
    original_data_points = m * n
    compressed_data_points = k * (m + n + 1) # Sはベクトルなので k 個
    compression_ratio = compressed_data_points / original_data_points * 100
    println("ランク $(k): $(round(compression_ratio, digits=2))%")
end

解説

• 圧縮率の目安は、必要な特異値と特異ベクトルの要素数から計算されます。
• Plots.jl を使って、元の画像と圧縮された画像を比較表示します。
• clamp01! は、再構成された画像のピクセル値が0.0から1.0の範囲に収まるように調整します。
• ranks で指定された異なる数の特異値を使って画像を再構成します。ランクが低いほど圧縮率は高いですが、画質は低下します。
• svd() でSVDを計算します。
• 画像をGray.(testimage("cameraman"))で読み込み、Float64.(channelview(img_gray))で数値行列に変換します。

SVDは、正方行列でない、または特異な（逆行列を持たない）行列の擬似逆行列（Moore-Penrose pseudoinverse）を計算するためにも使用できます。Juliaでは pinv() 関数がこの計算を行います。

using LinearAlgebra

# 例として特異な行列を定義
A_singular = [1.0 2.0;
              2.0 4.0] # ランク1の特異行列

println("特異な行列 A_singular:")
println(A_singular)

# Juliaの標準関数 `pinv()` を使って擬似逆行列を計算
A_pinv_builtin = pinv(A_singular)
println("\n`pinv()` で計算された擬似逆行列:")
println(A_pinv_builtin)

# SVDを使って手動で擬似逆行列を計算
F_pinv = svd(A_singular)
U_pinv, S_pinv, V_pinv = F_pinv.U, F_pinv.S, F_pinv.V

# 特異値の逆数を取る (0に近い特異値は逆数を0にする)
tolerance = eps(maximum(S_pinv)) * max(size(A_singular)...) # 小さな特異値を0とみなす閾値
S_inv = zeros(length(S_pinv))
for i in 1:length(S_pinv)
    if S_pinv[i] > tolerance
        S_inv[i] = 1.0 / S_pinv[i]
    end
end

# S_invを対角行列に変換 (転置して、S_invの対角成分はS_pinvの転置と同じ次元になる)
S_plus = zeros(size(A_singular, 2), size(A_singular, 1)) # n x m のゼロ行列
for i in 1:length(S_inv)
    S_plus[i, i] = S_inv[i]
end

# 擬似逆行列 A^+ = V * S^+ * U'
A_pinv_manual = V_pinv * S_plus * U_pinv'

println("\nSVDを使って手動で計算された擬似逆行列:")
println(A_pinv_manual)

# 2つの結果がほぼ同じであることを確認
println("\n手動計算とbuiltinの差のノルム:")
println(norm(A_pinv_manual - A_pinv_builtin))

# 擬似逆行列の性質を確認: A * A^+ * A ≈ A
println("\nA * A_pinv_manual * A:")
println(A_singular * A_pinv_manual * A_singular)
println("\nA * A_pinv_manual * A - A のノルム:")
println(norm(A_singular * A_pinv_manual * A_singular - A_singular))

# 最小二乗問題の解法
# Ax = b の最小二乗解は x = A^+ b で与えられる
b = [1.0, 3.0] # 解が存在しないベクトル
x_least_squares = A_pinv_manual * b
println("\nAx = b の最小二乗解 x:")
println(x_least_squares)
println("A * x_least_squares (bの近似):")
println(A_singular * x_least_squares)

• 擬似逆行列は、連立一次方程式 Ax=b の解が存在しない場合（過剰決定系など）に、最小二乗の意味で最適な解 x を見つけるために使用されます。
• S+ は S の対角成分の逆数を取ったもので、非常に小さい特異値は0として扱われます。この閾値処理は数値安定性のために重要です。eps(maximum(S_pinv)) * max(size(A_singular)...) のようなヒューリスティックな閾値がよく使われます。
• SVDによる擬似逆行列の計算式は $A^+ = V S^+ U^\*$ です。
• pinv(A) は、Juliaの組み込み関数で簡単に擬似逆行列を計算できます。

---

LinearAlgebra モジュールには、svd() が返す SVD オブジェクトから特定の情報だけを取得するための関数が用意されています。