つるかめ算メモ その2

数値的に連立方程式を解く関数をクラスにしてみた
なんかごちゃごちゃしてるので もう少し交通整理したほうが良いと思うけど実際のユースケースに応じて修正したほうが良いかも。

  • 違う問題をたくさん解くのか
  • 同じ係数の問題で入力が違うのを何度も解くのか
  • 連立方程式の項の数(係数の数)がどれぐらいか
    • 正規化した係数をキャッシュしたほうが良いのかメモリを節約しないといけないぐらい大きい問題か、

などなど
(これ、もともとは数学セミナーに掲載されてたCTスキャンのデータを数値的に解くのに使われてたアルゴリズムなので、結構大きい問題もありうる
あんまり大きい問題をPythonで解くのはどうかと思うけど)


import numpy as np
class Turukame(object):
def __init__(self,a,b,x=None,err=1e-3):
L=self.L=np.linalg.norm(a,axis=1)
self.a=a.copy()
self.b=b.copy()
self.o=np.zeros(a.shape,dtype=b.dtype)
self.err=err
if x:
self.x=x.copy()
else:
self.x=np.random.randn( b.shape[0])

for i in range(len(a)):
self.a[i]=a[i]/L[i]
self.b[i]=b[i]/L[i]
self.o[i]=self.a[i]*self.b[i]

assert (np.linalg.norm(a[i]))

def solve(self,x=None):
if x:
x=x.copy()
else:
x=self.x
for i in range(len(self.a)):
a=self.a[i]
b=self.b[i]
o=self.o[i]
x=x-np.dot(x-o,a)*a
self.x=x
return x
def valid(self,x=None,err=None):
if x is None:
x=self.x
e=self.sqerr(x)
if err is None:
err=self.err
if e