読者です 読者をやめる 読者になる 読者になる

LU分解

はじめに

LU分解(エルユーぶんかい)とは、正方行列  A を下三角行列 $ L $ と上三角行列 $ U $ の積に分解すること。すなわち $ A = LU $ が成立するような $ L $ と $ U $ を求めることをいう。 (LU分解 - Wikipediaより)

この分解は

等に使われることでよく知られており、 効率的に計算することが出来るようになります。

C言語による最新アルゴリズム事典 (ソフトウェアテクノロジー)

C言語による最新アルゴリズム事典 (ソフトウェアテクノロジー)

アルゴリズム

LU分解にはGauss法を用います。

hades-netherworld-service.hatenablog.com

次のような係数行列が連立方程式等から与えられたとします。 $$ \begin{eqnarray} A &= \begin{pmatrix} 2 & 5 & 7 \\ 4 & 13 & 20 \\ 8 & 29 & 50 \end{pmatrix} \end{eqnarray} $$ このとき、Gauss法を適用すると、 $$ \begin{eqnarray} A' &= \begin{pmatrix} 2 & 5 & 7 \\ 0 & 3 & 6 \\ 0 & 0 & 4 \end{pmatrix} \end{eqnarray} $$ となります。
この適用では、行列の行をある一定の順序、ルールで2,4,3倍したものを他の行から引くことで、 上のように左下部分を0にする という処理を行っています。

どの0を作るために何倍したのかという情報を、0の代わりに$ A' $に書き込むと、 $$ \begin{eqnarray} A'' &= \begin{pmatrix} 2 & 5 & 7 \\ 2 & 3 & 6 \\ 4 & 3 & 4 \end{pmatrix} \end{eqnarray} $$ となり、コレ一つを持っておけば、 どのような右辺が連立方程式に与えられたとしてもgauss法を右辺部分のデータにだけ適応させて、解を求めていくことができます。 便利ですね。

さて、本題のLU分解ですが、 求められた$ A'' $を、左下部分$ L $(対角成分は1にする)と右上部分$ U $に分けると、最初の行列$ A $は$ L $と$ U $を書けたものに相当します。 (すごいですねー)

つまり、 $$ \begin{eqnarray} \begin{split} A &= L U \\ \begin{pmatrix} 2 & 5 & 7 \\ 4 & 13 & 20 \\ 8 & 29 & 50 \end{pmatrix} &= \begin{pmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 4 & 3 & 1 \end{pmatrix} \begin{pmatrix} 2 & 5 & 7 \\ 0 & 3 & 6 \\ 0 & 0 & 4 \end{pmatrix} \end{split} \end{eqnarray} $$ となります。 これがLU分解です。

問題点

たしかにここまでで、LU分解のアルゴリズム自体は終わりです。 しかし、Gauss法を係数行列に適用する際、対角成分が0に近づくと、不整合を起こしてしまうという問題が発生してしまいます。

この回避方法として、ピボット法というのがあり、今回ではこれも盛り込んでみたいと思います。

ピボット選択

ピボット選択では、$ k $行目の計算時に、$ | a_ik | (i=k, k+1, \dots, n)$が最大となる行を$ k $行目と入れ替える手法です。 例えば、 $$ \begin{eqnarray} \begin{cases} 2x_1 + 5x_2 + 7x_3 &= 23 \\ 4x_1 + 13x_2 + 20x_3 &= 58 \\ 8x_1 + 29x_2 + 50x_3 &= 132 \end{cases}
\end{eqnarray} $$ という連立方程式が与えられた時、
1行目には $ 8x_1 + 29x_2 + 50x_3 = 132 $ を選択します(第1要素の係数が最大)。
その結果、対角行列が0ないしはそれに近い値になりすぎてしまうことを回避することができます。

コードにしてみる

#!/usr/bin/env python
# -*- coding: utf-8 -*-
import numpy as np

class LU:
  def __init__(self, matrix):
    if not len(matrix) == len(matrix[0]):
      raise Exception('matrix must be square')
    # copy matrix as float value data
    self.matrix = matrix.astype(float)
    self.ip = np.array(range(0, len(self.matrix)))

    for i in range(0, len(self.matrix) - 1):
      self._pibot(i)
      index = self.ip[i]
      for j in range(i + 1, len(self.matrix)):
        t = self.matrix[self.ip[j]][i] / self.matrix[index][i]
        for k in range(i, len(self.matrix)):
          self.matrix[self.ip[j]][k] -= t * self.matrix[index][k]
        self.matrix[self.ip[j]][i] = t

  def _pibot(self, i):
    index = self.ip[i]
    maxIndex = i
    maxValue = abs(self.matrix[index][index])
    for j in range(i + 1, len(self.matrix)):
      roopIndex = self.ip[j]
      if abs(self.matrix[roopIndex][index]) > maxValue:
        maxIndex = j
        maxValue = abs(self.matrix[roopIndex][index])
    if not i == maxIndex:
      self.ip[maxIndex] = self.ip[i] - self.ip[maxIndex]
      self.ip[i] -= self.ip[maxIndex]
      self.ip[maxIndex] += self.ip[i]

  def getL(self):
    l = len(self.matrix)
    res = np.identity(l)
    for i in range(1, l):
      for j in range(0, i):
        res[i][j] = self.matrix[self.ip[i]][j]
    return res

  def getU(self):
    l = len(self.matrix)
    res = np.identity(l)
    for i in range(0, l):
      for j in range(i, l):
        res[i][j] = self.matrix[self.ip[i]][j]
    return res

このコードではピボット選択をする際に配列を直接入れ替えるのではなく、
ip変数を用意してどのようにピボットを行ったかが分かるようにしています。
こうしておくことで、例えば連立方程式を解きたいという時に、
ピボットして位置が変わっていても正確に対応付けることが出来るようになります。

動かしてみる

では実際に動かしてみましょう

qa = np.array([[2,5,7],[4,13,20],[8,29,50]])
lu = LU(qa)
print(lu.ip)
print(lu.getL(), '\n')
print(lu.getU())
>>> [2 0 1]
>>> [[ 1.          0.          0.        ]
>>>  [ 0.25        1.          0.        ]
>>>  [ 0.5         0.66666667  1.        ]]
>>>
>>> [[  8.          29.          50.        ]
>>>  [  0.          -2.25        -5.5       ]
>>>  [  0.           0.          -1.33333333]]

※LとUの計算時には、ピボットして行列を変換した後の形で値を返すようにしています。

こうして得られたLとUの積を見てみると、

[[  8.   29.            50.         ]
 [  2.   5.             7.          ]
 [  4.   12.9999999925  19.999999985]]
途中で割り切れない値が出てしまったため、若干ずれてしまっていますが、無限精度という優しい目で見てあげれば一致します。

## おわりに
今回はピボット選択に、行だけの入れ替えを行う`部分ピボット選択`を用いましたが、  
列も入れ替えの対象とする`完全ピボット選択`というものがあります。  
ただ、概ね部分ピボット選択で事足りるようです。