Juliaのコレスキー分解：potrf!()とcholesky()の性能比較と最適化

JuliaのLinearAlgebra.LAPACK.potrf!()関数は、LAPACK（Linear Algebra PACKage）のルーチンをJuliaから呼び出すことで、行列のコレスキー分解（Cholesky decomposition）を行う関数です。特に、この関数は正方行列の下三角行列または上三角行列への分解をインプレースで行います。つまり、元の行列そのものが分解された結果で上書きされます。

もう少し詳しく説明します。

• 返り値
  potrf!()は、分解が成功したかどうかを示す整数値を返します。通常、0は成功を意味し、0以外の値はエラーを示します。 しかし、分解後の行列そのものがAに格納されているため、通常は返り値よりもAの中身を確認します。

• ! (感嘆符)
  Juliaの慣習で、関数名に!が付いている場合、その関数が引数を変更（破壊的変更）することを意味します。potrf!()は、入力行列Aを分解後の行列で上書きするため、!が付いています。

• potrf!()の引数
  

  • UPLO (Character): 'U' (Upper) または 'L' (Lower) を指定します。
    • 'U': 入力行列Aの上三角部分を使ってコレスキー分解を行い、上三角行列Uを返します（A = U^HU）。この場合、Aの下三角部分は変更されません。
    • 'L': 入力行列Aの下三角部分を使ってコレスキー分解を行い、下三角行列Lを返します（A = LL^H）。この場合、Aの上三角部分は変更されません。
  • A (Matrix): コレスキー分解を行う正方行列です。この行列は関数によって上書きされます。 つまり、分解後のLまたはUがこの場所に格納されます。

• コレスキー分解とは
  正定値エルミート行列Aを、下三角行列Lとその転置（または共役転置）L^T（またはL^H）の積、つまり A = LL^T (または A = LL^H) に分解することを指します。 potrf!()は、この分解を計算します。

例

using LinearAlgebra

A = [4.0 1.0 1.0; 1.0 5.0 2.0; 1.0 2.0 6.0]  # 正定値エルミート行列

# 下三角行列への分解
L = potrf!('L', A) # AがLで上書きされる
println("L = \n", L)
println("A = \n", A) # Aの内容が変化していることを確認

# 上三角行列への分解 (元のAを復元してから実行)
A = [4.0 1.0 1.0; 1.0 5.0 2.0; 1.0 2.0 6.0]
U = potrf!('U', A) # AがUで上書きされる
println("U = \n", U)
println("A = \n", A) # Aの内容が変化していることを確認

---

よくあるエラー

1. ArgumentError: matrix must be square
   入力行列Aが正方行列でない場合に発生します。potrf!()は正方行列に対してのみ定義されています。

2. PosDefException
   入力行列Aが正定値行列でない場合に発生します。コレスキー分解は正定値行列に対してのみ可能です。半正定値行列や不定値行列に対しては失敗します。

3. DimensionMismatch
   UPLOに指定した文字('U'または'L')と、行列Aの形状が矛盾する場合に発生することがあります。例えば、'U'を指定したのに、Aの下三角部分に分解に必要な情報がない場合などです。

4. LAPACKのエラーコード
   potrf!()は、内部でLAPACKのルーチンを呼び出しています。LAPACKのルーチンがエラーを返した場合、Julia側で捕捉され、エラーメッセージとして表示されます。これらのエラーコードは、LAPACKのマニュアルを参照することで、詳細な原因を調べることができます。

トラブルシューティング

---

基本的なコレスキー分解 (下三角行列)

using LinearAlgebra

A = [4.0 1.0 1.0; 1.0 5.0 2.0; 1.0 2.0 6.0]  # 正定値対称行列
L = potrf!('L', A) # 下三角行列 L を計算 (A は L で上書きされる)

println("L (下三角行列):\n", L)
println("A (上書きされた行列):\n", A) # A の下三角部分が L で上書きされている

# L * L' が元の A に等しいことを確認 (誤差を考慮)
println("L * L' (再構成):\n", L * L')
println("元の A:\n", [4.0 1.0 1.0; 1.0 5.0 2.0; 1.0 2.0 6.0])

# 元の A を保持したい場合はコピーを作成
A_original = [4.0 1.0 1.0; 1.0 5.0 2.0; 1.0 2.0 6.0]
A_copy = copy(A_original)
L_copy = potrf!('L', A_copy)
println("L_copy:\n", L_copy)
println("A_original:\n", A_original) # A_original は変更されない
println("A_copy:\n", A_copy) # A_copy は L_copy で上書きされる

基本的なコレスキー分解 (上三角行列)

using LinearAlgebra

A = [4.0 1.0 1.0; 1.0 5.0 2.0; 1.0 2.0 6.0]  # 正定値対称行列
U = potrf!('U', A) # 上三角行列 U を計算 (A は U で上書きされる)

println("U (上三角行列):\n", U)
println("A (上書きされた行列):\n", A) # A の上三角部分が U で上書きされている

# U' * U が元の A に等しいことを確認 (誤差を考慮)
println("U' * U (再構成):\n", U' * U)
println("元の A:\n", [4.0 1.0 1.0; 1.0 5.0 2.0; 1.0 2.0 6.0])

エラー処理

using LinearAlgebra

A = [1.0 2.0; 3.0 4.0; 5.0 6.0] # 正方行列ではない
try
    potrf!('L', A)
catch e
    println("エラーが発生しました: ", e)
end

A = [1.0 2.0; 2.0 1.0] # 正定値ではない
try
    potrf!('L', A)
catch e
    println("エラーが発生しました: ", e)
end

A = [4.0 1.0; 1.0 4.0] # 正定値
try
  L = potrf!('L', A)
  println("コレスキー分解成功:\n", L)
catch e
  println("エラーが発生しました: ", e)
end

コレスキー分解を用いた線形方程式の求解

using LinearAlgebra

A = [4.0 1.0 1.0; 1.0 5.0 2.0; 1.0 2.0 6.0] # 正定値対称行列
b = [1.0, 2.0, 3.0] # 右辺ベクトル

L = potrf!('L', copy(A)) # コレスキー分解 (A のコピーを使用)
y = L \ b # 前進代入
x = L' \ y # 後退代入

println("解 x:\n", x)

# A * x が b に等しいことを確認 (誤差を考慮)
println("A * x:\n", A * x)
println("b:\n", b)

using LinearAlgebra

A = [4.0 1.0 1.0; 1.0 5.0 2.0; 1.0 2.0 6.0] # 正定値対称行列

L = potrf!('L', copy(A)) # コレスキー分解 (A のコピーを使用)
Linv = inv(L)
Ainv = Linv' * Linv

println("A の逆行列:\n", Ainv)

# A * Ainv が単位行列に近いことを確認 (誤差を考慮)
println("A * Ainv:\n", A * Ainv)

---

cholesky()関数

Juliaの標準ライブラリLinearAlgebraには、cholesky()関数が用意されています。これはpotrf!()と同様にコレスキー分解を行う関数ですが、いくつかの点で異なります。

• 使いやすさ
  cholesky()はpotrf!()よりも高レベルな関数であり、より簡単に使うことができます。
• 戻り値
  cholesky()はCholesky型のオブジェクトを返します。このオブジェクトには、分解された行列（LまたはU）や、分解の種類（下三角または上三角）などの情報が含まれています。
• インプレース操作
  cholesky()はデフォルトではインプレース操作を行いません。つまり、元の行列を変更せずに分解結果を返します。必要であれば、cholesky!(A)のように!をつけてインプレース操作を行うことも可能です。

using LinearAlgebra

A = [4.0 1.0 1.0; 1.0 5.0 2.0; 1.0 2.0 6.0]

# cholesky() でコレスキー分解
C = cholesky(A)

# 下三角行列 L を取得
L = C.L
println("L:\n", L)

# 上三角行列 U を取得
U = C.U
println("U:\n", U)

# 分解の種類を確認
println("分解の種類: ", C.UPLO)

# インプレース操作
A_copy = copy(A)
C_inplace = cholesky!(A_copy)
println("A_copy (上書き):\n", A_copy)

# 線形方程式の求解
b = [1.0, 2.0, 3.0]
x = C \ b  # C は Cholesky 型のオブジェクトなので、\ 演算子で直接解を求められる
println("解 x:\n", x)


# 逆行列の計算
Ainv = inv(C) # C は Cholesky 型のオブジェクトなので、inv() で直接逆行列を計算できる
println("Ainv:\n", Ainv)

lu()関数

コレスキー分解は正定値行列に特化した分解方法ですが、正定値でない行列に対してはLU分解を用いることができます。lu()関数は、行列を下三角行列Lと上三角行列Uに分解します。

using LinearAlgebra

A = [1.0 2.0; 3.0 4.0] # 正定値行列でなくてもOK
F = lu(A)

L = F.L
U = F.U
P = F.P # 置換行列 (Permutation matrix)

println("L:\n", L)
println("U:\n", U)
println("P:\n", P)

# A == P * L * U
println("P * L * U:\n", P * L * U)
println("A:\n", A)

JuliaのLinearAlgebraパッケージには、他にも様々な行列分解関数が用意されています。例えば、QR分解 (qr())、SVD分解 (svd()) などがあります。これらの分解は、コレスキー分解が適用できない場合に利用することができます。

外部パッケージ

Juliaのパッケージエコシステムには、高度な行列計算や特化した分解方法を提供するパッケージが存在します。例えば、GenericLinearAlgebra.jl などがあります。これらのパッケージを利用することで、より複雑な行列計算を行うことができます。