行列の累乗

python

行列の累乗を計算する方法。
実行時間はgoogle colabのGPUなしで計測。

Blocks POJ No.3734 (蟻本P182)などで使用できる。

①対角化して累乗。数学的な感じ。
#input

import numpy as np
np.set_printoptions(precision = 0)

A = np.array([[2,1,0],[2,2,2],[0,1,2]])

def diagonalization(x):
  eig = np.linalg.eig(x)
  e = np.diag(eig[0])
  p = eig[1]
  return e,p

e = diagonalization(A)
D = e[0]
P = e[1]
P_inv = np.linalg.inv(P)

B = P@(D**3)@P_inv
print("Aの3乗：")
print(B)

import time
t = time.time()
C = P@(D**(10**18))@P_inv
print("Aの10**18乗：")
print(C)
print("処理時間：")
print(time.time()-t)

#output

Aの3乗
[[20. 16. 12.]
 [32. 32. 32.]
 [12. 16. 20.]]
Aの10**18乗
[[inf nan nan]
 [nan inf inf]
 [nan inf inf]]
処理時間
0.004506111145019531

→累乗が大きいとオーバーフローしている。


②numpyのライブラリを使う。
#input

import numpy as np
np.set_printoptions(precision = 0)

A = np.array([[2,1,0],[2,2,2],[0,1,2]])

B = np.linalg.matrix_power(A,3)
print("Aの3乗：")
print(B)

import time
t = time.time()
C = np.linalg.matrix_power(A,10**18)
print("Aの10**18乗：")
print(C)
print("処理時間：")
print(time.time()-t)

#output

Aの3乗：
[[20 16 12]
 [32 32 32]
 [12 16 20]]
Aの10**18乗：
[[0 0 0]
 [0 0 0]
 [0 0 0]]
処理時間：
0.003349781036376953

→累乗が大きいとオーバーフローする。
意外とタイムアウトはしない。スカラーの累乗とは計算方法が違うのだろうか。


③numpyのdot()を使って計算。
単純に繰り返しをおこなう。

#input

import numpy as np

def mat_power(A,num):
  result = A
  for j in range(num-1):
    result = result.dot(A)
  return result

B = mat_power(A,3)
print("Aの3乗：")
print(B)

import time
t = time.time()
C = mat_power(A,10**18)
print("Aの10**18乗：")
print(C)
print("処理時間：")
print(time.time()-t)

#output

Aの3乗：
[[20 16 12]
 [32 32 32]
 [12 16 20]]

→計算時間10分以上で途中で手動中断。
繰り返し部分がめちゃくちゃ時間がかかる。


④繰り返し二乗法を使って計算。
np.dot()を使用した。

#input

import numpy as np

A = np.array([[2,1,0],[2,2,2],[0,1,2]])

def linearpow(A,n):
  if n == 1:
    return A

  if n%2 == 1:
    B = np.dot(A,linearpow(A,n-1))
    return B
  
  B = linearpow(A,n//2)
  B = np.dot(B,B)
  return B


B = linearpow(A,3)
print("Aの3乗：")
print(B)

import time
t = time.time()
C = linearpow(A,10**18)
print("Aの10**18乗：")
print(C)
print("処理時間：")
print(time.time()-t)

#output

Aの3乗：
[[20 16 12]
 [32 32 32]
 [12 16 20]]
Aの10**18乗：
[[0 0 0]
 [0 0 0]
 [0 0 0]]
処理時間：
0.003695964813232422

→なんにしてもオーバーフローはする。速度はmatrix_power()と変わらない。ただ、この方法が一番わかりやすくて自作プログラムで使うなら間違いがなさそう。


⑤結果のmodを取りたい場合のやりかた。
繰り返し二乗法とnp.dot()を使った。

#input()

import numpy as np

A = np.array([[2,1,0],[2,2,2],[0,1,2]])

mod = 10

def linearmodpow(A,n,mod):
  if n == 1:
    return A%mod

  if n%2 == 1:
    B = np.dot(A,linearmodpow(A,n-1,mod))
    return B%mod
  
  B = linearmodpow(A,n//2,mod)
  B = np.dot(B,B)
  return B%mod

B = linearmodpow(A,3,mod)
print("Aの3乗のmod：")
print(B)

import time
t = time.time()
C = linearmodpow(A,10**100,mod)
print("Aの10**18乗のmod：")
print(C)
print("処理時間：")
print(time.time()-t)

#output

Aの3乗のmod：
[[0 6 2]
 [2 2 2]
 [2 6 0]]
Aの10**18乗のmod：
[[2 4 6]
 [8 8 8]
 [6 4 2]]
処理時間：
0.00898122787475586

→計算結果は合っているようす？
計算時間も速い。
オーバーフローせずに済んでいる。
numpy.linalg.matrix_power()は途中でmodを挟めないので、
解答としてmodを出すならこちらで用は足りそう。

numpyはCの型の制約を受ける。とのこと。
オーバーフローには注意。

参考：
対角化とmatrix_power()について：
assam-blog.com

np.dot()について：
amaharawaya.livedoor.blog

オーバーフローについて：
ikatakos.com