M.Hiroi's Home Page
http://www.geocities.jp/m_hiroi/

Python3 Programming

お気楽 NumPy プログラミング超入門

[ Home | Light | Python3 ]

●数で遊ぼう

今回は簡単な数理パズルを出題します。Python でも簡単にプログラムを作ることができますが、できるだけ NumPy を使ってみてください。

●多角数

点を多角形の形に並べたとき、その総数を「多角数 (polygonal number)」といいます。詳しい説明は拙作のページ Puzzle DE Programming 多角数 をお読みください。

  1. 三角数と四角数を求めるプログラムを作ってください。
  2. 三角錐数と四角錐数を求めるプログラムを作ってください。
  3. p 角数とを p 角錐数求めるプログラムを作ってください。
  4. 三角数かつ四角数である数を「平方三角数 (square triangular number)」といいます。1,500,000 以下の平方三角数をすべて求めてください。

解答

●フィボナッチ数

フィボナッチ (fibonacci) 数はイタリアの数学者レオナルド・フィボナッチにちなんで名付けられた数です。詳しい説明は拙作のページ Puzzle DE Programming フィボナッチ数 をお読みください。

  1. n 以下のフィボナッチ数列を求めるプログラムを作ってください。
  2. 300,000,000 以下で最も大きいフィボナッチ数 Fn を求めてください。
  3. 初項 F0 から n 番目までのフィボナッチ数の和 F0 + F1 + F2 + ... + Fn を考えます。300,000,000 以下のフィボナッチ数の総和を求めてください。
  4. 300,000,000 以下のフィボナッチ数で、偶数になる項の総和を求めてください。
  5. 300,000,000 を 1 つ以上の連続しない相異なるフィボナッチ数の和として表してください。たとえば、7 = 1 + 1 + 2 + 3 と表すことができますが、1 を 2 回使っていることと、1, 2, 3 は連続したフィボナッチ数 (F2, F3, F4) なので条件を満たしていません。7 = 2 + 5 とすると条件を満たします。

解答

●素数

「素数 (prime number)」は正の約数が 1 と自分自身しかない自然数のことです。1 は素数に含めません。詳しい説明は拙作のページ Puzzle DE Programming 素数 をお読みください。

  1. 3,000,000 以下の素数の個数とその最大値を求めてください。
  2. 3,000,001 を素因数分解してください。素因数分解とは、素数でない整数 (合成数) を素数の積の形に書き表すことです。たとえば、12 は 22 * 3 と素因数分解することができます。
  3. 3,000,000 以下のフィボナッチ素数 (フィボナッチ数で素数) の個数とその最大値を求めてください。
  4. 3,000,000 以下の素数で、差の最大値とその素数の組を求めてください。
  5. 差が 2 である素数の組を「双子素数 (twin prime)」といいます。3,000,000 以下の双子素数の個数とその最大値を求めてください。

解答

●約数

約数 (divisor) は整数 n を割り切る整数のことです。詳しい説明は拙作のページ Puzzle DE Programming 約数 をお読みください。

  1. 自然数 n の約数の個数を求めるプログラムを作ってください。
  2. 自然数 n の約数の和を求めるプログラムを作ってください。
  3. 自然数 n の約数をすべて求めるプログラムを作ってください。

解答

●ハミング数

7 以上の素数で割り切れない正の整数を「ハミング数 (Hamming Numbers)」といいます。ハミング数は素因子が 2, 3, 5 しかない自然数のことで、素因数分解したとき 2i * 3j * 5k (i, j, k >= 0) の形式になります。たとえば、100 以下のハミング数は次のようになります。

1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36, 40, 45, 48, 50, 
54, 60, 64, 72, 75, 80, 81, 90, 96, 100

N 以下のハミング数を求めるプログラムを作ってください。

解答

●魔方陣

3 行 3 列の魔方陣を解くプログラムを作ってください。

 ┌─┬─┬─┐   式
 │A│B│C│   A + B + C = N, A + E + I = N
 ├─┼─┼─┤   D + E + F = N. C + E + G = N
 │D│E│F│   G + H + I = N
 ├─┼─┼─┤   A + D + G = N
 │G│H│I│   B + E + H = N
 └─┴─┴─┘   C + F + I = N

        図:魔方陣

上図の A から I の場所に 1 から 9 までの数字をひとつずつ配置します。縦横斜めの合計が等しくなるように数字を配置してください。

解答


●多角数の解答

リスト : 多角数の解答

import numpy as np

# 三角数と四角数
def triangular_number(n):
    return np.arange(1, n + 1, dtype=np.int64).cumsum()

def square_number(n):
    xs = np.arange(1, n + 1, dtype=np.int64)
    return xs ** 2

# 三角錐数と四角錐数
def triangular_pyramidal(n):
    return triangular_number(n).cumsum()

def square_pyramidal(n):
    return square_number(n).cumsum()

# 多角数と多角錐数
def polygonal_number(p, n):
    xs = np.arange(1, n + 1, dtype=np.int64)
    return ((p - 2) * xs ** 2 - (p - 4) * xs) // 2

def polygonal_pyramidal(p, n):
    return polygonal_number(p, n).cumsum()


# 平方三角数
def square_triangular_number(n):
    r = np.roots(np.poly1d([1, 1, -2 * n]))
    x = np.floor(r[1])
    y = np.floor(np.sqrt(n))
    return np.intersect1d(triangular_number(x), square_number(y))

print(triangular_number(20))
print(square_number(20))
print(triangular_pyramidal(20))
print(square_pyramidal(20))
print(polygonal_number(5, 20))
print(polygonal_pyramidal(5, 20))
print(polygonal_number(6, 20))
print(polygonal_pyramidal(6, 20))
print(square_triangular_number(1500000))
[  1   3   6  10  15  21  28  36  45  55  66  78  91 105 120 136 153 171 190 210]  # 三角数
[  1   4   9  16  25  36  49  64  81 100 121 144 169 196 225 256 289 324 361 400]  # 四角数
[   1    4   10   20   35   56   84  120  165  220  286  364  455  560   680       # 三角錐数
  816  969 1140 1330 1540]
[   1    5   14   30   55   91  140  204  285  385  506  650  819 1015             # 四角錐数
 1240 1496 1785 2109 2470 2870]
[  1   5  12  22  35  51  70  92 117 145 176 210 247 287 330 376 425 477 532 590]  # 五角数
[   1    6   18   40   75  126  196  288  405  550  726  936 1183 1470             # 五角錐数
 1800 2176 2601 3078 3610 4200]
[  1   6  15  28  45  66  91 120 153 190 231 276 325 378 435 496 561 630 703 780]  # 六角数
[   1    7   22   50   95  161  252  372  525  715  946 1222 1547 1925             # 六角錐数
 2360 2856 3417 4047 4750 5530]
[      1      36    1225   41616 1413721]                                          # 平方三角数

●フィボナッチ数の解答

リスト : フィボナッチ数の解答

# フィボナッチ数
def fibo(n):
    if n == 0 or n == 1:
        return n
    else:
        f = np.matrix([[1, 1], [1, 0]], 'int64')
        return (f ** (n - 1))[0, 0]

# n 以下のフィボナッチ数列
def fibo_seq(n):
    fn = np.frompyfunc(fibo, 1, 1)
    phi = (1 + np.sqrt(5)) / 2
    x = np.floor(np.log(n * np.sqrt(5)) / np.log(phi))
    return fn(np.arange(x + 1, dtype='int64'))

# n 以下で最大のフィボナッチ数
def fibo_max(n):
    return fibo_seq(n)[-1]

# n 以下のフィボナッチ数の総和
def fibo_sum(n):
    return np.sum(fibo_seq(n))

# n 以下のフィボナッチ数の偶数の総和
def fibo_even_sum(n):
    xs = fibo_seq(n)
    return np.sum(xs[xs % 2 == 0])

# 整数を 1 つ以上の連続しない相異なるフィボナッチ数の和で表す
# ゼッケンドルフの表現
def solver(n):
    xs = []
    ys = fibo_seq(n + 1)
    while n > 0:
        m = ys[-1]
        xs.append(m)
        n -= m
        ys = ys[ys <= n]
    return xs

print(fibo_seq(30000))
print(fibo_sum(300000000))
print(fibo_even_sum(300000000))
print(solver(300000000))
print(solver(267914296))
print(solver(267914296 - 1))
[0 1 1 2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181 6765 10946 17711 28657]
701408732
350704366
[267914296, 24157817, 5702887, 2178309, 46368, 233, 89, 1]
[267914296]
[165580141, 63245986, 24157817, 9227465, 3524578, 1346269, 514229, 196418, 75025, 28657, 
 10946, 4181, 1597, 610, 233, 89, 34, 13, 5, 2]

●素数の解答

リスト : 素数の解答

# 素数
def sieve(n):
    ps = np.ones((n - 3) // 2 + 1, 'u1')
    x = 3
    while x * x <= n:
        y = (x - 3) // 2
        if ps[y]: ps[y + x::x] = 0
        x += 2
    xs = np.where(ps == 1)[0]
    return np.insert(xs * 2 + 3, 0, 2)

# 素因数分解
def factorization(n):
    xs = []
    ys = []
    for p in sieve(np.int(np.sqrt(n))):
        if n < p * p: break
        c = 0
        while n % p == 0:
            c += 1
            n //= p
        if c > 0:
            xs.append(p)
            ys.append(c)
    if n > 1:
        xs.append(n)
        ys.append(1)
    return np.array(xs), np.array(ys, dtype=np.int64)

# フィボナッチ素数
def fibo_prime(n):
    xs = fibo_seq(n)
    ys = sieve(n)
    return np.intersect1d(xs, ys)

# 素数の個数と最大値
xs = sieve(3000000)
print(len(xs), np.max(xs))

# 素因数分解
print(factorization(3000001))

# フィボナッチ素数
a = fibo_prime(3000001)
print(a)
print(len(a), a[-1])

# 素数の差
ys = np.diff(xs)
m = np.argmax(ys)
print(np.max(ys), xs[m], xs[m + 1])
print(np.histogram(ys, np.arange(0, 151, 10)))

# 双子素数
zs = xs[np.where(ys == 2)[0]]
print(len(zs), np.max(zs))
216816 2999999
(array([ 853, 3517], dtype=int64), array([1, 1], dtype=int64))
[2 3 5 13 89 233 1597 28657 514229]
9 514229
148 2010733 2010881
(array([91347, 74257, 28769, 13778,  5391,  1855,   853,   360,   121,
          50,    21,     6,     4,     2,     1], dtype=int64), array([  0,  10,  20,  30,  40,  50,
        60,  70,  80,  90, 100, 110, 120, 130, 140, 150]))
20932 2999831

●約数の解答

リスト : 約数の解答

# 約数の個数
def divisor_num(n):
    xs = factorization(n)[1] + 1
    return np.prod(xs)

# 約数の合計値
def divisor_sum(n):
    xs, ys = factorization(n)
    a = 1
    for x, y in zip(xs, ys):
        a *= np.sum(x ** np.arange(y + 1))
    return a

# 約数
def divisor(n):
    xs, ys = factorization(n)
    zs = xs[0] ** np.arange(ys[0] + 1)
    for x, y in zip(xs[1:], ys[1:]):
        zs = np.outer(zs, x ** np.arange(y + 1)).flatten()
    zs.sort()
    return zs

print(divisor_num(24))
print(divisor_num(12345678))
print(divisor_num(123456789))
print(divisor_num(1234567890))
print(divisor_num(1111111111))

print(divisor_sum(24))
print(divisor_sum(12345678))
print(divisor_sum(123456789))
print(divisor_sum(1234567890))
print(divisor_sum(1111111111))

print(divisor(24))
print(divisor(12345678))
print(divisor(123456789))
print(divisor(1234567890))
print(divisor(1111111111))
8
24
12
48
16

60
27319968
178422816
3211610688
1246404096

[ 1  2  3  4  6  8 12 24]

[       1        2        3        6        9       18       47       94
      141      282      423      846    14593    29186    43779    87558
   131337   262674   685871  1371742  2057613  4115226  6172839 12345678]

[        1         3         9      3607      3803     10821     11409
     32463     34227  13717421  41152263 123456789]

[         1          2          3          5          6          9
         10         15         18         30         45         90
       3607       3803       7214       7606      10821      11409
      18035      19015      21642      22818      32463      34227
      36070      38030      54105      57045      64926      68454
     108210     114090     162315     171135     324630     342270
   13717421   27434842   41152263   68587105   82304526  123456789
  137174210  205761315  246913578  411522630  617283945 1234567890]

[         1         11         41        271        451       2981
       9091      11111     100001     122221     372731    2463661
    4100041   27100271  101010101 1111111111]

●ハミング数の解答

ハミング数を求める一番簡単な方法は 1 から n までの数列を生成し、その要素が 2, 3, 5 だけで割り切れるか試すことです。プログラムはとても簡単ですが、引数 n の値が大きくなると時間がかかるようになります。n に比べてハミング数の個数は少ないようなので、式 2i * 3j * 5k (i, j, k >= 0) を使ってハミング数を生成したほうがよさそうです。この式の計算は NumPy の meshgrid() を使うと簡単です。プログラムと結果を以下に示します。

リスト : ハミング数の解答

def hamming(n):
    def init(m):
        xs = [1]
        k = m
        while k <= n:
            xs.append(k)
            k *= m
        return np.array(xs, dtype=np.float64)
    x, y, z = np.meshgrid(init(2), init(3), init(5))
    a = (x * y * z).flat
    return np.sort(a[a <= n])

print(hamming(100))
for x in [1000, 10000, 100000, 1000000, 10000000]:
    print(x, len(hamming(x)))
[  1.   2.   3.   4.   5.   6.   8.   9.  10.  12.  15.  16.  18.  20.
  24.  25.  27.  30.  32.  36.  40.  45.  48.  50.  54.  60.  64.  72.
  75.  80.  81.  90.  96. 100.]
1000 86
10000 175
100000 313
1000000 507
10000000 768

●魔方陣の解答

下記のプログラムは単純な生成検定法なので、実行時間はちょっと遅いです。

リスト : 魔方陣 (3 * 3 盤)

import itertools

def magic_square():
    for xs in itertools.permutations(range(1, 9 + 1)):
        a = np.array(xs).reshape((3, 3))
        n = a.trace()
        if a[:, ::-1].trace() != n: continue
        b = a.sum(axis=0)
        if not np.all(b == n): continue
        b = a.sum(axis=1)
        if not np.all(b == n): continue
        print(a)

magic_square()
[[2 7 6]
 [9 5 1]
 [4 3 8]]
[[2 9 4]
 [7 5 3]
 [6 1 8]]
[[4 3 8]
 [9 5 1]
 [2 7 6]]
[[4 9 2]
 [3 5 7]
 [8 1 6]]
[[6 1 8]
 [7 5 3]
 [2 9 4]]
[[6 7 2]
 [1 5 9]
 [8 3 4]]
[[8 1 6]
 [3 5 7]
 [4 9 2]]
[[8 3 4]
 [1 5 9]
 [6 7 2]]

●連立一次方程式の解法

今回は連立一次方程式を解く基本的なアルゴリズムを取り上げます。NumPy には連立一次方程式を解く関数 linalg.solve() が用意されているので、実用的にはそちらを使えばいいのですが、Python (NumPy) とアルゴリズムのお勉強ということで、あえてプログラムを作ってみましょう。

●ガウスの消去法

ガウスの消去法 (Gaussian elimination) は人が連立方程式を解くときの方法とほとんど同じです。次の図を見てください。

a1 * x + a2 * y + a3 * z = d1       a1 * x + a2  * y + a3  * z = d1        a1 * x + a2  * y + a3   * z = d1
b1 * x + b2 * y + b3 * z = d2  ==>           b2' * y + b3' * z = d2'  ==>           b2' * y + b3'  * z = d2'
c1 * x + c2 * y + c3 * z = d3                c2' * y + c3' * z = d3'                          c3'' * z = d3''

                                          図 : 前進消去

1 番目の式を b1 / a1 倍して 2 番目の式から引き算すると、x の項を消去することができます。同様の方法で 3 番目の式の x の項を消去することができます。次に、2 番目の式を c2' / b2' 倍して 3 番目の式から引き算すると y の項を消去することができます。この処理を「前進消去」といいます。

前進消去を行うと、最後の式には変数が z しかありません。z の値は d3'' / c3'' に決定することができます。次はひとつ前の式に戻り、 変数 z に d3'' / c3'' を代入すると変数 y の値を求めることができます。これを繰り返して、先頭の式に戻ると変数 x の値を求めることができます。この処理を「後退代入」といいます。

連立一次方程式は行列を使って Ax = b と表すことができます。A を係数行列といい、A と b を連結した行列を拡大係数行列といいます。

[[a1, a2, a3],    [[a1, a2, a3, d1],
 [b1, b2, b3],     [b1, b2, b3, d2],
 [c1, c2, c3]]     [c1, c2, c3, d3]]

   係数行列         拡大係数行列

ガウスの消去法は、拡大係数行列を使うと簡単にプログラムを作ることができます。Python でプログラムすると、次のようになります。

リスト : ガウスの消去法

def gauss(xs, ys):
    # 拡大係数行列の生成
    n = len(xs)
    zs = [x + [y] for x, y in zip(xs, ys)]
    # 前進消去
    for i in range(n - 1):
        for j in range(i + 1, n):
            temp = zs[j][i] / zs[i][i]
            # 係数を 0 にする計算は省略可
            # その場合は range(i + 1, n + 1) とする
            for k in range(i, n + 1):
                zs[j][k] -= temp * zs[i][k]
        # for debug
        print(zs)
    # 後退代入
    # 解は zs の n 列目にセットする
    for i in range(n - 1, -1, -1):
        for j in range(i + 1, n): zs[i][n] -= zs[i][j] * zs[j][n]
        zs[i][n] /= zs[i][i]
    return [z[-1] for z in zs]

#  x +  y = 100
# 2x + 4y = 272
print(gauss([[1, 1], [2, 4]], [100, 272]))

#  x +  y +  z = 10
# 2x + 4y + 6z = 38
# 2x      + 4z = 14
print(gauss([[1, 1, 1], [2, 4, 6], [2, 0, 4]], [10, 38, 14]))

#   a +  b +   c + d = -5
#  -a +  b -   c + d = -7
#  8a + 4b +  2c + d = -31
# -8a + 4b + -2c + d = -35
print(gauss([[1, 1, 1, 1], [-1, 1, -1, 1], [8, 4, 2, 1], [-8, 4, -2, 1]], [-5, -7, -31, -35]))

#   a -  b +  c - d + e = 1
# 12a - 6b + 2c         = 0 
#   a +  b +  c + d + e = 8
# 12a + 6b + 2c         = 0
#  4a + 3b + 2c + d     = 1
print(gauss([[1,-1,1,-1,1], [12,-6,2,0,0], [1,1,1,1,1], [12,6,2,0,0], [4,3,2,1,0]], [1,0,8,0,1]))

関数 gauss() の引数 xs が係数行列、ys が右辺の定数を格納したリストです。まず最初に、拡大係数行列を生成して変数 zs にセットします。あとはアルゴリズムをそのままプログラムしただけなので、特に難しいところはないと思います。実行結果は次のようになります。なお、拡大行列の表示は手作業で整形しています。

[[1,   1,  100],
 [0.0, 2.0, 72.0]]
[64.0, 36.0]

[[1,    1,   1,   10],
 [0.0,  2.0, 4.0, 18.0],
 [0.0, -2.0, 2.0, -6.0]]
[[1,   1,   1,   10],
 [0.0, 2.0, 4.0, 18.0],
 [0.0, 0.0, 6.0, 12.0]]
[3.0, 5.0, 2.0]

[[1,    1,    1,    1,    -5],
 [0.0,  2.0,  0.0,  2.0, -12.0],
 [0.0, -4.0, -6.0, -7.0,   9.0],
 [0.0, 12.0,  6.0,  9.0, -75.0]]
[[1,   1,    1,    1,    -5],
 [0.0, 2.0,  0.0,  2.0, -12.0],
 [0.0, 0.0, -6.0, -3.0, -15.0],
 [0.0, 0.0,  6.0, -3.0,  -3.0]]
[[1,   1,    1,    1,    -5],
 [0.0, 2.0,  0.0,  2.0, -12.0],
 [0.0, 0.0, -6.0, -3.0, -15.0],
 [0.0, 0.0,  0.0, -6.0, -18.0]]
[0.0, -9.0, 1.0, 3.0]

[[1, -1, 1, -1, 1, 1],
 [0.0, 6.0, -10.0, 12.0, -12.0, -12.0],
 [0.0, 2.0, 0.0, 2.0, 0.0, 7.0],
 [0.0, 18.0, -10.0, 12.0, -12.0, -12.0],
 [0.0, 7.0, -2.0, 5.0, -4.0, -3.0]]
[[1, -1, 1, -1, 1, 1],
 [0.0, 6.0, -10.0, 12.0, -12.0, -12.0],
 [0.0, 0.0, 3.333333333333333, -2.0, 4.0, 11.0],
 [0.0, 0.0, 20.0, -24.0, 24.0, 24.0],
 [0.0, 0.0, 9.666666666666668, -9.0, 10.0, 11.0]]
[[1, -1, 1, -1, 1, 1],
 [0.0, 6.0, -10.0, 12.0, -12.0, -12.0],
 [0.0, 0.0, 3.333333333333333, -2.0, 4.0, 11.0],
 [0.0, 0.0, 0.0, -11.999999999999998, -3.552713678800501e-15, -42.000000000000014],
 [0.0, 0.0, 0.0, -3.1999999999999984, -1.6000000000000032, -20.90000000000001]]
[[1, -1, 1, -1, 1, 1],
 [0.0, 6.0, -10.0, 12.0, -12.0, -12.0],
 [0.0, 0.0, 3.333333333333333, -2.0, 4.0, 11.0],
 [0.0, 0.0, 0.0, -11.999999999999998, -3.552713678800501e-15, -42.000000000000014],
 [0.0, 0.0, 0.0, -4.440892098500626e-16, -1.6000000000000023, -9.70000000000001]]
[0.3124999999999991, 0.0, -1.874999999999997, 3.5, 6.062499999999997]

NumPy を使ったプログラムと実行結果を以下に示します。

リスト : ガウスの消去法 (NumPy バージョン)

import numpy as np

def gauss_np(xs, ys):
    # 拡大係数行列の生成
    n = len(xs)
    zs = np.c_[xs.astype(np.float_), ys.astype(np.float_)]
    # 前進消去
    for i in range(n - 1):
        for j in range(i + 1, n):
            temp = zs[j, i] / zs[i, i]
            zs[j, i:] -= temp * zs[i, i:]
        # for debug
        print(zs)
    # 後退代入
    for i in range(n - 1, -1, -1):
        zs[i, n] -= zs[i, i+1:n] @ zs[i+1:, n]
        zs[i, n] /= zs[i, i]
    return zs[:, n]

print(gauss_np(np.array([[1, 1], [2, 4]]), np.array([100, 272])))
print(gauss_np(np.array([[1, 1, 1], [2, 4, 6], [2, 0, 4]]), np.array([10, 38, 14])))
print(gauss_np(np.array([[1, 1, 1, 1], [-1, 1, -1, 1], [8, 4, 2, 1], [-8, 4, -2, 1]]), np.array([-5, -7, -31, -35])))
print(gauss_np(np.array([[1,-1,1,-1,1], [12,-6,2,0,0], [1,1,1,1,1], [12,6,2,0,0], [4,3,2,1,0]]), np.array([1,0,8,0,1])))
[[   1.    1.  100.]
 [   0.    2.   72.]]
[ 64.  36.]
[[  1.   1.   1.  10.]
 [  0.   2.   4.  18.]
 [  0.  -2.   2.  -6.]]
[[  1.   1.   1.  10.]
 [  0.   2.   4.  18.]
 [  0.   0.   6.  12.]]
[ 3.  5.  2.]
[[  1.   1.   1.   1.  -5.]
 [  0.   2.   0.   2. -12.]
 [  0.  -4.  -6.  -7.   9.]
 [  0.  12.   6.   9. -75.]]
[[  1.   1.   1.   1.  -5.]
 [  0.   2.   0.   2. -12.]
 [  0.   0.  -6.  -3. -15.]
 [  0.   0.   6.  -3.  -3.]]
[[  1.   1.   1.   1.  -5.]
 [  0.   2.   0.   2. -12.]
 [  0.   0.  -6.  -3. -15.]
 [  0.   0.   0.  -6. -18.]]
[ 0. -9.  1.  3.]
[[  1.  -1.   1.  -1.   1.   1.]
 [  0.   6. -10.  12. -12. -12.]
 [  0.   2.   0.   2.   0.   7.]
 [  0.  18. -10.  12. -12. -12.]
 [  0.   7.  -2.   5.  -4.  -3.]]
[[  1.          -1.           1.          -1.           1.           1.        ]
 [  0.           6.         -10.          12.         -12.         -12.        ]
 [  0.           0.           3.33333333  -2.           4.          11.        ]
 [  0.           0.          20.         -24.          24.          24.        ]
 [  0.           0.           9.66666667  -9.          10.          11.        ]]
[[  1.00000000e+00  -1.00000000e+00   1.00000000e+00  -1.00000000e+00
    1.00000000e+00   1.00000000e+00]
 [  0.00000000e+00   6.00000000e+00  -1.00000000e+01   1.20000000e+01
   -1.20000000e+01  -1.20000000e+01]
 [  0.00000000e+00   0.00000000e+00   3.33333333e+00  -2.00000000e+00
    4.00000000e+00   1.10000000e+01]
 [  0.00000000e+00   0.00000000e+00   0.00000000e+00  -1.20000000e+01
   -3.55271368e-15  -4.20000000e+01]
 [  0.00000000e+00   0.00000000e+00   0.00000000e+00  -3.20000000e+00
   -1.60000000e+00  -2.09000000e+01]]
[[  1.00000000e+00  -1.00000000e+00   1.00000000e+00  -1.00000000e+00
    1.00000000e+00   1.00000000e+00]
 [  0.00000000e+00   6.00000000e+00  -1.00000000e+01   1.20000000e+01
   -1.20000000e+01  -1.20000000e+01]
 [  0.00000000e+00   0.00000000e+00   3.33333333e+00  -2.00000000e+00
    4.00000000e+00   1.10000000e+01]
 [  0.00000000e+00   0.00000000e+00   0.00000000e+00  -1.20000000e+01
   -3.55271368e-15  -4.20000000e+01]
 [  0.00000000e+00   0.00000000e+00   0.00000000e+00  -4.44089210e-16
   -1.60000000e+00  -9.70000000e+00]]
[ 0.3125  0.     -1.875   3.5     6.0625]

●ガウス・ジョルダン法

ガウスの消去法は、係数行列の部分を上三角行列に変換することで、連立一次方程式を解きました。ここで、係数行列を単位行列に変換できれば、後退代入することなく解を求めることができます。これを「ガウス・ジョルダン法 (Gauss - Jordan elimination)」といいます。簡単な動作例を下図に示します。

次の連立方程式をガウス・ジョルダン法で解く

  x +  y +  z = 10
 2x + 4y + 6z = 38
 2x      + 4z = 14

拡大係数行列を zs とする

(1) zs[0, 0] を 1 にする
[[  1.   1.   1.  10.],    # zs[0, 0:] /= 1
 [  2.   4.   6.  38.],
 [  2.   0.   4.  14.]]

(2) zs[1, 0], zs[2, 0] を 0 にする
[[  1.   1.   1.  10.]     
 [  0.   2.   4.  18.]     # zs[1, 0:] -= zs[0, 0:] * 2
 [  0.  -2.   2.  -6.]]    # zs[2, 0:] -= zs[0, 0:] * 2

(3) zs[1, 1] を 1 にする
[[  1.   1.   1.  10.]
 [  0.   1.   2.   9.]     # zs[1, 1:] /= 2
 [  0.  -2.   2.  -6.]]

(4) zs[0, 1], zs[2, 1] を 0 にする
[[  1.   0.  -1.   1.]     # zs[0, 1:] -= zs[1, 1:] * 1
 [  0.   1.   2.   9.]
 [  0.   0.   6.  12.]]    # zs[2, 1:] -= zs[1, 1:] * -2

(5) zs[2, 2] を 1 にする
[[  1.   0.  -1.   1.]     
 [  0.   1.   2.   9.]
 [  0.   0.   1.   2.]]    # zs[2, 2:] /= 6

(6) zs[2, 0], zs[2, 1] を 0 にする
[[  1.   0.   0.   3.]     # zs[0, 2:] -= zs[2, 2:] * -1
 [  0.   1.   0.   5.]     # zs[1, 2:] -= zs[2, 2:] * 2
 [  0.   0.   1.   2.]]

このように、ガウス・ジョルダン法の基本的な考え方は実にシンプルで、プログラムも簡単になるのですが、速度はガウスの消去法よりも少し遅くなるといわれています。以下にプログラムと実行例を示します。

リスト : ガウス・ジョルダン法

import numpy as np

def gauss_jordan(xs, ys):
    n = len(xs)
    zs = np.c_[xs.astype(np.float_), ys.astype(np.float_)]
    for i in range(n):
        # zs[i, i] を 1 にする
        temp = zs[i, i]
        zs[i, i:] /= temp
        # zs[i] 以外の行の i 番目の係数を 0 にする
        for j in range(n):
            if i != j:
                temp = zs[j, i]
                zs[j, i:] -= zs[i, i:] * temp
        # for debug
        print(zs)
    return zs[:, n]

print(gauss_jordan(np.array([[1, 1], [2, 4]]), np.array([100, 272])))
print(gauss_jordan(np.array([[1, 1, 1], [2, 4, 6], [2, 0, 4]]), np.array([10, 38, 14])))
print(gauss_jordan(np.array([[1, 1, 1, 1], [-1, 1, -1, 1], [8, 4, 2, 1], [-8, 4, -2, 1]]), np.array([-5, -7, -31, -35])))
print(gauss_jordan(np.array([[1,-1,1,-1,1], [12,-6,2,0,0], [1,1,1,1,1], [12,6,2,0,0], [4,3,2,1,0]]), np.array([1,0,8,0,1])))
[[   1.    1.  100.]
 [   0.    2.   72.]]
[[  1.   0.  64.]
 [  0.   1.  36.]]
[ 64.  36.]

[[  1.   1.   1.  10.]
 [  0.   2.   4.  18.]
 [  0.  -2.   2.  -6.]]
[[  1.   0.  -1.   1.]
 [  0.   1.   2.   9.]
 [  0.   0.   6.  12.]]
[[ 1.  0.  0.  3.]
 [ 0.  1.  0.  5.]
 [ 0.  0.  1.  2.]]
[ 3.  5.  2.]

[[  1.   1.   1.   1.  -5.]
 [  0.   2.   0.   2. -12.]
 [  0.  -4.  -6.  -7.   9.]
 [  0.  12.   6.   9. -75.]]
[[  1.   0.   1.   0.   1.]
 [  0.   1.   0.   1.  -6.]
 [  0.   0.  -6.  -3. -15.]
 [  0.   0.   6.  -3.  -3.]]
[[  1.    0.    0.   -0.5  -1.5]
 [  0.    1.    0.    1.   -6. ]
 [ -0.   -0.    1.    0.5   2.5]
 [  0.    0.    0.   -6.  -18. ]]
[[ 1.  0.  0.  0.  0.]
 [ 0.  1.  0.  0. -9.]
 [ 0.  0.  1.  0.  1.]
 [-0. -0. -0.  1.  3.]]
[ 0. -9.  1.  3.]

[[  1.  -1.   1.  -1.   1.   1.]
 [  0.   6. -10.  12. -12. -12.]
 [  0.   2.   0.   2.   0.   7.]
 [  0.  18. -10.  12. -12. -12.]
 [  0.   7.  -2.   5.  -4.  -3.]]
[[  1.           0.          -0.66666667   1.          -1.          -1.        ]
 [  0.           1.          -1.66666667   2.          -2.          -2.        ]
 [  0.           0.           3.33333333  -2.           4.          11.        ]
 [  0.           0.          20.         -24.          24.          24.        ]
 [  0.           0.           9.66666667  -9.          10.          11.        ]]
[[  1.    0.    0.    0.6  -0.2   1.2]
 [  0.    1.    0.    1.    0.    3.5]
 [  0.    0.    1.   -0.6   1.2   3.3]
 [  0.    0.    0.  -12.    0.  -42. ]
 [  0.    0.    0.   -3.2  -1.6 -20.9]]
[[ 1.   0.   0.   0.  -0.2 -0.9]
 [ 0.   1.   0.   0.   0.   0. ]
 [ 0.   0.   1.   0.   1.2  5.4]
 [-0.  -0.  -0.   1.  -0.   3.5]
 [ 0.   0.   0.   0.  -1.6 -9.7]]
[[ 1.      0.      0.      0.      0.      0.3125]
 [ 0.      1.      0.      0.      0.      0.    ]
 [ 0.      0.      1.      0.      0.     -1.875 ]
 [-0.     -0.     -0.      1.      0.      3.5   ]
 [-0.     -0.     -0.     -0.      1.      6.0625]]
[ 0.3125  0.     -1.875   3.5     6.0625]

●ピボット選択

ところで、ガウスの消去法とガウス・ジョルダン法は、拡大係数行列の対角成分 zs[i, i] が 0 になると計算できなくなります。たとえば、次の連立方程式は解くことができますが、zs[0, 0] が 0 になっているため、プログラムの実行中にエラー (divide by zero) が発生します。

      2y + 4z = 14   [[0.  2.  4. 14.]
  x +  y +  z = 10    [1.  1.  1. 10.]
 4x + 2y + 6z = 38    [4.  2.  6. 38.]]

    連立方程式        拡大係数行列 zs

このとき、拡大係数行列の行を交換すると連立一次方程式が解ける場合があります。また、zs[i, i] が 0 に等しくなくても 0 に近い値だと解の誤差が大きくなります。そこで、zs[i, i] の絶対値がなるべく大きくなるように行を交換します。これを「部分ピボット選択」といいます。なお、参考文献 1 によると、『行だけでなく列も交換する完全ピボット選択という方法もあるが、通常は部分ピボット選択で十分である』 とのことです。

部分ピボット選択のプログラムは NumPy を使うと簡単です。次のリストを見てください。

リスト : 部分ピボット選択

def select_pivot(xs, i):
    k = np.abs(xs[i:,i]).argmax() + i
    if k != i:
        temp = xs[i].copy()
        xs[i] = xs[k]
        xs[k] = temp

# ガウスの消去法
def gauss1(xs, ys):
    # 拡大係数行列の生成
    n = len(xs)
    zs = np.c_[xs.astype(np.float_), ys.astype(np.float_)]
    # 前進消去
    for i in range(n - 1):
        # ピボット選択
        select_pivot(zs, i)
        for j in range(i + 1, n):
            temp = zs[j, i] / zs[i, i]
            zs[j, i:] -= temp * zs[i, i:]
        # for debug
        print(zs)
    # 後退代入
    for i in range(n - 1, -1, -1):
        zs[i, n] -= zs[i, i+1:n] @ zs[i+1:, n]
        zs[i, n] /= zs[i, i]
    return zs[:, n]

# ガウス・ジョルダン法
def gauss_jordan1(xs, ys):
    # 拡大係数行列の生成
    n = len(xs)
    zs = np.c_[xs.astype(np.float_), ys.astype(np.float_)]
    for i in range(n):
        # ピボット選択
        select_pivot(zs, i)
        temp = zs[i, i]
        zs[i, i:] /= temp
        for j in range(n):
            if i != j:
                temp = zs[j, i]
                zs[j, i:] -= zs[i, i:] * temp
        # for debug
        print(zs)
    return zs[:, n]

print(gauss1(np.array([[0, 2, 4], [1, 1, 1], [4, 2, 6]]), np.array([14, 10, 38])))
print(gauss_jordan1(np.array([[0, 2, 4], [1, 1, 1], [4, 2, 6]]), np.array([14, 10, 38])))

ピボット選択は関数 select_pivot() で行います。np.abs() で xs[i:, i] の絶対値を求め、その中から最大値の位置を argmax() で求めます。実際の位置はそれに i を加算したものになります。これを変数 k にセットします。k と i が異なる値であれば i 行と k 行と交換します。Python のリストであれば xs[i], xs[k] = xs[k], xs[i] で交換できるのですが、NumPy の配列は交換できません。NumPy のスライス操作はビューを生成するだけなので、copy() を使って行の値をコピーする必要があります。ご注意くださいませ。

あとは select_pivot() を呼び出すだけです。それでは実際に試してみましょう。

[[  4.    2.    6.   38. ]
 [  0.    0.5  -0.5   0.5]
 [  0.    2.    4.   14. ]]
[[  4.    2.    6.   38. ]
 [  0.    2.    4.   14. ]
 [  0.    0.   -1.5  -3. ]]
[ 5.  3.  2.]

[[  1.    0.5   1.5   9.5]
 [  0.    0.5  -0.5   0.5]
 [  0.    2.    4.   14. ]]
[[ 1.   0.   0.5  6. ]
 [ 0.   1.   2.   7. ]
 [ 0.   0.  -1.5 -3. ]]
[[ 1.  0.  0.  5.]
 [ 0.  1.  0.  3.]
 [ 0.  0.  1.  2.]]
[ 5.  3.  2.]

正しく解を求めることができました。

●逆行列

連立方程式 Ax = b を解く場合、A の逆行列 A-1 がわかれば、x = A-1b で解を求めることができます。実をいうと、逆行列はガウス・ジョルダン法を使って簡単に求めることができるのです。次の図を見てください。

逆行列の定義 : A @ A-1 = I

[[a11, a12, a13],     [[b11, b12, b13],     [[1, 0, 0],
 [a21, a22, a23],  @   [b21, b22, b23],  =   [0, 1, 0],
 [a31, a32, a33]]      [b31, b32, b33]]      [0, 0, 1]]

    係数行列               逆行列            単位行列 I

上記の定義により逆行列 A-1 は次の方程式を解けば求めることができる

A @ [b11, b21, b31] = [1, 0, 0]
A @ [b12, b22, b32] = [0, 1, 0]
A @ [b13, b23, b33] = [0, 0, 1]

これは次の拡大係数行列にガウス・ジョルダン法を適用すること同じ

[[a11, a12, a13, 1, 0, 0],
 [a21, a22, a23, 0, 1, 0],
 [a31, a32, a33, 0, 0, 1]]

拡大係数行列の右半分に逆行列が求まる

プログラムも簡単です。次のリストを見てください。

リスト : ガウス・ジョルダン法で逆行列を求める

def matinv_gj(xs):
    # 拡大係数行列の生成
    n = len(xs)
    zs = np.c_[xs.astype(np.float_), np.eye(n)]
    for i in range(n):
        # ピボットの選択
        select_pivot(zs, i)
        temp = zs[i, i]
        zs[i, i:] /= temp
        for j in range(n):
            if i != j:
                temp = zs[j, i]
                zs[j, i:] -= zs[i, i:] * temp
    return zs[:, n:]

# 簡単なテスト
a1 = np.array([[1, 1], [2, 4]])
a2 = np.array([[1, 1, 1], [2, 4, 6], [2, 0, 4]])
a3 = np.array([[1, 1, 1, 1], [-1, 1, -1, 1], [8, 4, 2, 1], [-8, 4, -2, 1]])
a4 = np.array([[1, -1, 1, -1, 1], [12, -6, 2, 0, 0], [1, 1, 1, 1, 1], [12, 6, 2, 0, 0], [4, 3, 2, 1, 0]])
a5 = np.array([[0, 2, 4], [1, 1, 1], [4, 2, 6]])
for a in [a1, a2, a3, a4, a5]:
    b = matinv_gj(a)
    print(b)
    print(a @ b)

拡大係数行列を作るとき、右半分に単位行列 eye(n) を連結します。あとはガウス・ジョルダン法で解を求め、最後に拡大係数行列の右半分を返します。それでは実行してみましょう。

[[ 2.  -0.5]
 [-1.   0.5]]
[[ 1.  0.]
 [ 0.  1.]]

[[ 1.33333333 -0.33333333  0.16666667]
 [ 0.33333333  0.16666667 -0.33333333]
 [-0.66666667  0.16666667  0.16666667]]
[[  1.00000000e+00  -2.77555756e-17   0.00000000e+00]
 [  0.00000000e+00   1.00000000e+00   0.00000000e+00]
 [  4.44089210e-16  -1.11022302e-16   1.00000000e+00]]

[[-0.16666667  0.16666667  0.08333333 -0.08333333]
 [-0.16666667 -0.16666667  0.16666667  0.16666667]
 [ 0.66666667 -0.66666667 -0.08333333  0.08333333]
 [ 0.66666667  0.66666667 -0.16666667 -0.16666667]]
[[  1.00000000e+00   0.00000000e+00   0.00000000e+00  -2.77555756e-17]
 [  0.00000000e+00   1.00000000e+00   0.00000000e+00   2.77555756e-17]
 [ -1.11022302e-16   0.00000000e+00   1.00000000e+00  -1.11022302e-16]
 [  0.00000000e+00  -1.11022302e-16  -8.32667268e-17   1.00000000e+00]]

[[ -6.25000000e-02   4.16666667e-02   6.25000000e-02   8.33333333e-02  -1.25000000e-01]
 [  0.00000000e+00  -8.33333333e-02   0.00000000e+00   8.33333333e-02   0.00000000e+00]
 [  3.75000000e-01  -1.38777878e-17  -3.75000000e-01  -2.50000000e-01   7.50000000e-01]
 [ -5.00000000e-01   8.33333333e-02   5.00000000e-01  -8.33333333e-02   2.22044605e-16]
 [  1.87500000e-01  -4.16666667e-02   8.12500000e-01   1.66666667e-01  -6.25000000e-01]]
[[  1.00000000e+00  -6.93889390e-18   2.22044605e-16  -2.77555756e-17  -3.33066907e-16]
 [  1.11022302e-16   1.00000000e+00   4.44089210e-16   1.11022302e-16   2.22044605e-16]
 [  5.55111512e-17   1.38777878e-17   1.00000000e+00   0.00000000e+00   1.11022302e-16]
 [  1.11022302e-16  -2.77555756e-17   4.44089210e-16   1.00000000e+00   2.22044605e-16]
 [  0.00000000e+00  -2.77555756e-17   4.44089210e-16   5.55111512e-17   1.00000000e+00]]

[[-0.33333333  0.33333333  0.16666667]
 [ 0.16666667  1.33333333 -0.33333333]
 [ 0.16666667 -0.66666667  0.16666667]]
[[  1.00000000e+00   0.00000000e+00   0.00000000e+00]
 [  2.77555756e-17   1.00000000e+00   2.77555756e-17]
 [  1.11022302e-16   0.00000000e+00   1.00000000e+00]]

行列と逆行列の積を計算すると、浮動小数点数の計算による誤差がありますが、単位行列になることがわかります。

ところで、実際に連立一次方程式を解くとき、逆行列を使用することはほとんどありません。わざわざ逆行列を求めるよりも、ガウスの消去法を使ったほうがより速く解を求めることができるからです。同じ係数行列の連立方程式を何度も解く場合でも、逆行列を求めるのではなく、次に説明する LU 分解を使うことが多いようです。

●LU分解

連立一次方程式 Ax = b を解くとき、あらかじめ係数行列 A を下三角行列 L (lower) と上三角行列 U (Upper) に分解しておくと、簡単に解くことができます。Ax = b は LUx = b と表すことができるので、Ux = y とおくと Ly = b を解くことで y を求めることができます。次に、求めた y を使って Ux = y を解けば、Ax = b の解 x を求めることができます。簡単な例を示しましょう。

  連立方程式

  x +  y +  z = 10
 2x + 4y + 6z = 38
 2x      + 4z = 14

 [[1.  1.  1.]     [[1.  0.  0.]    [[1.  1.  1.]
  [2.  4.  6.]   =  [2.  1.  0.]  @  [0.  2.  4.] 
  [2.  0.  4.]]     [2. -1.  1.]]    [0.  0.  6.]]

               (1) LU 分解
  
 [[1.  0.  0.  10.]    y1 = 10
  [2.  1.  0.  38.]    y2 = 38 - 2 * 10 = 18 
  [2. -1.  1.  14.]]   y3 = 14 - 2 * 10 + 18 = 12

               (2) 前進代入

 [[1.  1.  1.  10.]   x1 = 10 - x2 - x3 = 3
  [0.  2.  4.  18.]   x2 = (18 - 4 * x3) / 2 = 5
  [0.  0.  6.  12.]]  x3 = 12 / 6 = 2

               (3) 後退代入

LU 分解はガウスの消去法と同じアルゴリズムで求めることができます。U は前進消去で得た行列で、L の要素は前進消去するとき係数を 0 にするための倍率になります。実際に L @ U を計算すると元の行列に戻ります。Ly = b を解く場合、L は下三角行列なので前進代入で y を求めることができます。次に Ux = y を解く場合、U は上三角行列なので後退代入で x を求めることができます。

ガウスの消去法の計算量は O(N3) ですが、前進代入と後進代入の計算量は O(N2) です。同じ係数行列を何度も使う場合、あらかじめ係数行列を LU 分解しておけば、ガウスの消去法よりも効率的に連立方程式を解くことができるわけです。

実際には A を L と U の二つに分けるのではなく、U の 0 の部分に L の要素をセットした行列 A' を作ることにします。プログラムは次のようになります。

リスト : LU 分解

def lu(xs):
    n = len(xs)
    zs = xs.astype(np.float_)
    for i in range(n):
        for j in range(i + 1, n):
            temp = zs[j, i] / zs[i, i]
            zs[j, i+1:] -= temp * zs[i, i+1:]
            zs[j, i] = temp
    return zs

def lu_solver(xs, ys):
    n = len(xs)
    zs = ys.astype(np.float_)
    # 前進代入
    for i in range(1, n):
        zs[i] -= xs[i, :i] @ zs[:i]
    # 後退代入
    for i in range(n - 1, -1, -1):
        zs[i] -= xs[i, i+1:] @ zs[i+1:]
        zs[i] /= xs[i, i]
    return zs

# 簡単なテスト
a1 = np.array([[1, 1], [2, 4]])
a2 = np.array([[1, 1, 1], [2, 4, 6], [2, 0, 4]])
a3 = np.array([[1, 1, 1, 1], [-1, 1, -1, 1], [8, 4, 2, 1], [-8, 4, -2, 1]])
a4 = np.array([[1, -1, 1, -1, 1], [12, -6, 2, 0, 0], [1, 1, 1, 1, 1], [12, 6, 2, 0, 0], [4, 3, 2, 1, 0]])
b1 = np.array([100, 272])
b2 = np.array([10, 38, 14])
b3 = np.array([-5, -7, -31, -35])
b4 = np.array([1, 0, 8, 0, 1])
for a, b in [(a1, b1), (a2, b2), (a3, b3), (a4, b4)]:
    xs = lu(a)
    print(xs)
    print(lu_solver(xs, b))

関数 lu() は引数の係数行列 xs を LU 分解します。アルゴリズムはガウスの消去法と同じですが、係数を 0 にするための倍率 temp を zs[j, i] にセットするところが異なります。

関数 lu_solver() は LU 分解した行列 xs と右辺の定数を格納したベクトル ys を受け取り、その解を求めます。前進代入は後退代入とは逆に前から順番に値を求めていきます。行列 L の対角成分 L[i, i] は 1 なので、割り算する必要はありません。後退代入はガウスの消去法のプログラムと同じです。

それでは実行してみましょう。

[[ 1.  1.]
 [ 2.  2.]]
[ 64.  36.]

[[ 1.  1.  1.]
 [ 2.  2.  4.]
 [ 2. -1.  6.]]
[ 3.  5.  2.]

[[ 1.  1.  1.  1.]
 [-1.  2.  0.  2.]
 [ 8. -2. -6. -3.]
 [-8.  6. -1. -6.]]
[ 0. -9.  1.  3.]

[[  1.00000000e+00  -1.00000000e+00   1.00000000e+00  -1.00000000e+00   1.00000000e+00]
 [  1.20000000e+01   6.00000000e+00  -1.00000000e+01   1.20000000e+01  -1.20000000e+01]
 [  1.00000000e+00   3.33333333e-01   3.33333333e+00  -2.00000000e+00   4.00000000e+00]
 [  1.20000000e+01   3.00000000e+00   6.00000000e+00  -1.20000000e+01  -3.55271368e-15]
 [  4.00000000e+00   1.16666667e+00   2.90000000e+00   2.66666667e-01  -1.60000000e+00]]
[ 0.3125  0.     -1.875   3.5     6.0625]

正常に動作していますね。

●ピボット選択を使った LU 分解

LU 分解にピボット選択を適用することも簡単です。次のリストを見てください。

リスト : LU 分解 (ピボット選択バージョン)

def select_pivot(xs, idx, i):
    k = np.abs(xs[i:,i]).argmax() + i
    if k != i:
        temp = xs[i].copy()
        xs[i] = xs[k]
        xs[k] = temp
        idx[i], idx[k] = idx[k], idx[i]

def lu_pv(xs):
    n = len(xs)
    zs = xs.astype(np.float_)
    idx = list(range(n))
    for i in range(n):
        select_pivot(zs, idx, i)
        for j in range(i + 1, n):
            temp = zs[j, i] / zs[i, i]
            zs[j, i+1:] -= temp * zs[i, i+1:]
            zs[j, i] = temp
    return zs, idx

def lu_solver_pv(xs, idx, ys):
    return lu_solver(xs, ys[idx])

select_pivot() の引数 idx に交換した行の情報をセットします。係数行列の大きさを n とすると、idx の初期値は 0 から n - 1 までの整数を格納したリストになります。i 行と k 行を交換した場合は、idx[i] と idx[k] の値も交換します。関数 lu_pv() はピボット選択しながら LU 分解を行い、L と U を格納した行列 zs とリスト idx を返します。関数 lu_solver_pv() は、lu_pv() が返した行列 xs とインデックスリスト idx を受け取り、引数 ys のベクトルの要素を idx に合わせて並べ替えてから lu_solver() を呼び出すだけです。

それでは実際に試してみましょう。以下に簡単なテストと実行結果を示します。

リスト :  簡単なテスト

a1 = np.array([[1, 1], [2, 4]])
a2 = np.array([[1, 1, 1], [2, 4, 6], [2, 0, 4]])
a3 = np.array([[1, 1, 1, 1], [-1, 1, -1, 1], [8, 4, 2, 1], [-8, 4, -2, 1]])
a4 = np.array([[1, -1, 1, -1, 1], [12, -6, 2, 0, 0], [1, 1, 1, 1, 1], [12, 6, 2, 0, 0], [4, 3, 2, 1, 0]])
a5 = np.array([[0, 2, 4], [1, 1, 1], [4, 2, 6]])
b1 = np.array([100, 272])
b2 = np.array([10, 38, 14])
b3 = np.array([-5, -7, -31, -35])
b4 = np.array([1, 0, 8, 0, 1])
b5 = np.array([14, 10, 38])
for a, b in [(a1, b1), (a2, b2), (a3, b3), (a4, b4), (a5, b5)]:
    xs, idx = lu_pv(a)
    print(xs)
    print(idx)
    print(lu_solver_pv(xs, idx, b))
[[ 2.   4. ]
 [ 0.5 -1. ]]
[1, 0]
[ 64.  36.]

[[ 2.    4.    6.  ]
 [ 1.   -4.   -2.  ]
 [ 0.5   0.25 -1.5 ]]
[1, 2, 0]
[ 3.  5.  2.]

[[ 8.      4.      2.      1.    ]
 [-1.      8.      0.      2.    ]
 [ 0.125   0.0625  0.75    0.75  ]
 [-0.125   0.1875 -1.      1.5   ]]
[2, 3, 0, 1]
[ 0. -9.  1.  3.]

[[ 12.          -6.           2.           0.           0.        ]
 [  1.          12.           0.           0.           0.        ]
 [  0.33333333   0.41666667   1.33333333   1.           0.        ]
 [  0.08333333  -0.04166667   0.625       -1.625        1.        ]
 [  0.08333333   0.125        0.625       -0.23076923   1.23076923]]
[1, 3, 4, 0, 2]
[ 0.3125  0.     -1.875   3.5     6.0625]

[[ 4.    2.    6.  ]
 [ 0.    2.    4.  ]
 [ 0.25  0.25 -1.5 ]]
[2, 0, 1]
[ 5.  3.  2.]

正常に動作していますね。

●参考文献・URL

  1. 奥村晴彦,『C言語による最新アルゴリズム事典』, 技術評論社, 1991
  2. ガウスの消去法 - Wikipedia

Copyright (C) 2018 Makoto Hiroi
All rights reserved.

[ Home | Light | Python3 ]