作者alen332l (alen3321)
看板Python
標題Re: [問題] 關於反矩陣inv計算爆值是正常的嗎?
時間Thu Jun 23 11:02:36 2016
※ 引述《alen332l (alen3321)》之銘言：
: ※ 引述《ray0215 (❺❺❻❻)》之銘言：
: : 小弟我最近需要計算相關係數(correlation coefficient matri)矩陣的反矩陣
: : 例如
: : Martrix array([[ 1. , 0.92531576, 0.81998247, 0.92586894],
: : [ 0.92531576, 1. , 0.70974912, 0.72358554],
: : [ 0.81998247, 0.70974912, 1. , 0.72669914],
: : [ 0.92586894, 0.72358554, 0.72669914, 1. ]])
: : 這樣用numpy.linalg.inv算出來的值
: : 會計算出
: : array([[ 6.51500296e+15, -3.18743561e+15, -7.89299760e+14,
: : -3.15207310e+15],
: : [ -3.18743561e+15, 1.55943840e+15, 3.86161323e+14,
: : 1.54213745e+15],
: : [ -7.89299760e+14, 3.86161323e+14, 9.56245323e+13,
: : 3.81877116e+14],
: : [ -3.15207310e+15, 1.54213745e+15, 3.81877116e+14,
: : 1.52502845e+15]])
: : 每個都十五次方
: : 我想應該不是正確解答吧@@
: : 請問要怎麼解決
: : 先謝謝版上的大大了
: 你推文說要計算R^(-1)-I 的值
: R = correlation matrix
: I = idendity matrix
: 處理矩陣時
: 有一個原則
: 能不要反矩陣 就不要算 因為算反矩陣很昂貴 又不穩定
: 你原本題目應該是：
: R ( X + I ) = I ,其中X為你要求的矩陣.
: 那麼，這樣處理比較穩定:
: R X = (I - R)
: └已知┘ └已知┘
: X 切成column表達:
: X = [x1 │x2 │x3...] ,其中x1為向量,x2,x3...亦同
: (I-R)定義成 B 好了
: 也切成column:
: (I-R) =def= B = [b1│b2│b3...]
: 原問題變成
: R ‧ x1 = b1
: (nxn) (nx1) (nx1) →維度驗算
: R ‧ x2 = b2
: R ‧ x3 = b3
: 每個問題用 BLAS 算完後
: 再把 x1,x2,x3...拼回來
: 印象中是這樣,幾年前處理的印象...
: p.s.
: 1.回文如有錯請高手指正 已經有一段時間沒處理matrix了 謝謝
: 2.我沒實際算...如果算出來是合理值 那就是你原本用inv時會有誤差造成
: 行列式小時，算inv的誤差會累積
: 通常基礎線性代數不會講
: 數學系的純數課程也不會講
我把我剛才的回文寫成code
如下: 但答案還是有10的7次方就是了
----------------------------------------------------------------
import numpy as np
R =np.array([[ 1. , 0.92531576, 0.81998247, 0.92586894],
[ 0.92531576, 1. , 0.70974912, 0.72358554],
[ 0.81998247, 0.70974912, 1. , 0.72669914],
[ 0.92586894, 0.72358554, 0.72669914, 1. ]])
I = np.array([[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]])
B = I - R
print("R",R)
print("B",B)
X = np.linalg.solve(R, B)
print("X",X)
-----------------------------------------------
答案就是 X
--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 36.228.190.212
※ 文章網址: https://webptt.com/m.aspx?n=bbs/Python/M.1466650959.A.B2B.html
※ 編輯: alen332l (36.228.190.212), 06/23/2016 18:47:27