JuliaでLU分解するときの注意点

Yuto Miyatake
Jan 3, 2022 2 min read

JuliaでLU分解を利用して連立一次方程式を数値計算するとき，枢軸選択 (pivoting) に注意する必要があります．みるべき情報はここですが，備忘録も兼ねてまとめておきます．

julia> using LinearAlgebra

julia> A = [0 0 2; 1 0 1 ; 2 -1 1]
3×3 Array{Int64,2}:
 0   0  2
 1   0  1
 2  -1  1

julia> L, U = lu(A) # 素直にLU分解
LU{Float64,Array{Float64,2}}
L factor:
3×3 Array{Float64,2}:
 1.0  0.0  0.0
 0.5  1.0  0.0
 0.0  0.0  1.0
U factor:
3×3 Array{Float64,2}:
 2.0  -1.0  1.0
 0.0   0.5  0.5
 0.0   0.0  2.0

julia> L*U
3×3 Array{Float64,2}:
 2.0  -1.0  1.0
 1.0   0.0  1.0
 0.0   0.0  2.0

L*UがAと一致しません．そもそもAの (1,1) 成分が 0 ですから，通常のLU分解は初手から破綻してしまいます．では何が行われているかというと枢軸選択 (pivoting) が行われているわけです．要するに行の入れ替えが行われています．

julia> L, U, p = lu(A)
LU{Float64,Array{Float64,2}}
L factor:
3×3 Array{Float64,2}:
 1.0  0.0  0.0
 0.5  1.0  0.0
 0.0  0.0  1.0
U factor:
3×3 Array{Float64,2}:
 2.0  -1.0  1.0
 0.0   0.5  0.5
 0.0   0.0  2.0

julia> p # 置換を表すベクトル
3-element Array{Int64,1}:
 3
 2
 1

julia> b = [4;3;5] # 連立一次方程式の右辺ベクトル
3-element Array{Int64,1}:
 4
 3
 5

julia> A\b #　連立一次方程式の解
3-element Array{Float64,1}:
  1.0
 -1.0
  2.0

julia> x1 = U\(L\b) # pivotingを無視して計算すると...
3-element Array{Float64,1}:
  0.5
 -0.5
  2.5

julia> x2 = U\(L\b[p]) # pivotingを考慮して計算するにはこのようにすればよい
3-element Array{Float64,1}:
  1.0
 -1.0
  2.0

置換を表す行列も簡単に取り出せます．

julia> F = lu(A);

julia> F.P
3×3 Array{Float64,2}:
 0.0  0.0  1.0
 0.0  1.0  0.0
 1.0  0.0  0.0

julia> (F.P)*A - L*U
3×3 Array{Float64,2}:
 0.0  0.0  0.0
 0.0  0.0  0.0
 0.0  0.0  0.0

なお，枢軸選択を無視してLU分解することも可能です．

julia> lu(A,Val(false)) # A は上の3x3行列
ERROR: ZeroPivotException: factorization encountered one or more zero pivots. Consider switching to a pivoted LU factorization.

怒られます．(1,1) 成分が 0 なのですから当然です．しかし，枢軸選択を無視してもLU分解が動く例もあります．

julia> A = [1 1; 2 1] # この行列に対して
2×2 Array{Int64,2}:
 1  1
 2  1

julia> L,U　= lu(A) #　LU分解すると
LU{Float64,Array{Float64,2}}
L factor:
2×2 Array{Float64,2}:
 1.0  0.0
 0.5  1.0
U factor:
2×2 Array{Float64,2}:
 2.0  1.0
 0.0  0.5

julia> L*U # 枢軸選択が適用されますが
2×2 Array{Float64,2}:
 2.0  1.0
 1.0  1.0

julia> L,U = lu(A,Val(false)) # このようにすれば，枢軸選択は適用されません
LU{Float64,Array{Float64,2}}
L factor:
2×2 Array{Float64,2}:
 1.0  0.0
 2.0  1.0
U factor:
2×2 Array{Float64,2}:
 1.0   1.0
 0.0  -1.0

julia> L*U # 実際に L*U が A に一致します
2×2 Array{Float64,2}:
 1.0  1.0
 2.0  1.0

あえて枢軸選択を行わないメリットはあまりないとは思いますが，例えば，matlabの素朴なLU分解と比較するときなどには lu(A,Val(false)) は役立つかもしれません．

Julia 線形計算
Yuto Miyatake
Yuto Miyatake
Associate Professor

Applied Mathematics, Numerical Analysis, Computational Uncertain Quantification