numpyを使用しないで共役勾配法の実装を試みていますが途中の式の計算がうまくできず困っている.
numpyを使用しないで共役勾配法の解を算出するプログラムを作っています.
途中の a = float( np.dot(r0.T,r0) / np.dot(np.dot(p.T, A),p) )という計算をnpを使わないで計算したいのですがうまく変えずに困っております.
自分なりに内包表記やzipを使用して試行錯誤はしたのですが,通りません.
アドバイスいただけたら幸いです.
また,みにくいソースかもしれませんがご容赦ください.
発生している問題・エラーメッセージ
ValueError Traceback (most recent call last) <ipython-input-173-78082ced4a70> in <module> 1 #print(b.shape) ----> 2 ans = cgm(mat,vec,x0) 3 print(ans) <ipython-input-172-15436a5d72cf> in cgm(A, b, x_init) 14 15 for i in range(k < k+1): ---> 16 d1 = dot(list(map(list,zip(*r0))),r0) 17 print(d1) 18 d2 = dot(dot(list(map(list,zip(*p))),A),p) <ipython-input-170-9fb70407aafc> in dot(a, b) 5 for i,v in enumerate(a): 6 for j,u in enumerate(zip(*b)): ----> 7 tmp[i][j] = sum([x*y for x,y in zip(v,u)]) 8 9 return tmp <ipython-input-170-9fb70407aafc> in <listcomp>(.0) 5 for i,v in enumerate(a): 6 for j,u in enumerate(zip(*b)): ----> 7 tmp[i][j] = sum([x*y for x,y in zip(v,u)]) 8 9 return tmp ~/.pyenv/versions/3.6.5/lib/python3.6/site-packages/numpy/matrixlib/defmatrix.py in __mul__(self, other) 218 if isinstance(other, (N.ndarray, list, tuple)) : 219 # This promotes 1-D vectors to row vectors --> 220 return N.dot(self, asmatrix(other)) 221 if isscalar(other) or not hasattr(other, '__rmul__') : 222 return N.dot(self, other) ValueError: shapes (1,48) and (1,48) not aligned: 48 (dim 1) != 1 (dim 0)
該当のソースコード
Python
1import random 2from scipy.io import mminfo,mmread 3 4 5mat = mmread('bcsstm01.mtx').todense() 6print(type(mat)) 7print(mat.shape) 8 9random.seed(1) 10vec = [[random.randrange(10) for i in range(48)]] 11print(vec) 12print(type(vec)) 13 14x0 = [[0]*48] 15print(type(x0)) 16 17 18def dot(a,b): 19 20 tmp = [[0]*len(b[0]) for i in range(len(a))] 21 22 for i,v in enumerate(a): 23 for j,u in enumerate(zip(*b)): 24 tmp[i][j] = sum([x*y for x,y in zip(v,u)]) 25 26 return tmp 27 28 29def cgm(A, b, x_init): 30 x = x_init 31 32 d = dot(A,x)[0] 33 #print(d) 34 r0 = [e1 - e2 for e1,e2 in zip(b,d)] 35 #r0 = [e1 - e2 for e1,e2 in zip(b,(dot(A,x)[0]))] 36 #print(dot(A,x)) 37 #print(r0) 38 39 p = r0 40 k = 0 41 42 43 for i in range(k < k+1): 44 **d1 = dot(list(map(list,zip(*r0))),r0) 45 print(d1) 46 d2 = dot(dot(list(map(list,zip(*p))),A),p) 47 print(d2) 48 a = [e3 / e4 for e3,e4 in zip(d1,d2)]** 49 #a = dot(r0.T,r0) / dot(dot(p.T, A),p) 50 x = x + p*a 51 r1 = r0 - dot(A*a, p) 52 if np.linalg.norm(r1) < 1.0e-10: 53 k = k+1 54 return x 55 b = dot(r1.T, r1) / dot(r0.T, r0) 56 p = r1 + b * p 57 r0 = r1 58 k = k+1 59 return x 60 61ans = cgm(mat,vec,x0) 62print(ans)
試したこと
a = float( np.dot(r0.T,r0) / np.dot(np.dot(p.T, A),p) )という計算はnumpyを使用しているため下記のようにzipとリスト内包表記を使用したものに変更してみたがエラーは解決できませんでした.
**でくくってあるところの下のコメントa = dot(r0.T,r0) / dot(dot(p.T, A),p)というよな計算を実装したいです.
助言いただけたら幸いです.
補足情報(FW/ツールのバージョンなど)
Python version 3.5.6
jupyter notebook 使用
numpy sympyは使用しないで内積の計算を行いたいためnp.は使用できません.
回答2件
あなたの回答
tips
プレビュー