scipy 中的复杂矩阵,SVD 不能按预期工作?
SVD not working as expected with complex matrix in scipy?
我试图在 scipy 中使用 SVD 找到 M@x = 0
的近似非零解,其中 M 是一个复值 4x4 矩阵。
首先是一个玩具示例:
M = np.array([
[1,1,1,1],
[1,0,1,1],
[1,-1,0,0],
[0,0,0,1e-10]
])
U, s, Vh = scipy.linalg.svd(M)
print(s) # [2.57554368e+00 1.49380718e+00 3.67579714e-01 7.07106781e-11]
print(Vh[-1]) # [ 0.00000000e+00 2.77555756e-16 -7.07106781e-01 7.07106781e-01]
print(np.linalg.norm( M@Vh[-1] )) # 7.07106781193738e-11
所以在这种情况下,s中的最小(最后)值很小,对应的最后一列Vh[-1]是M@x=0
的近似解,其中M@Vh[-1]也很小,和 s[-1].
的阶数大致相同
现在的真实例子不是以同样的方式工作:
M = np.array([[ 1.68572560e-01-3.98053448e-02j, 5.61165939e-01-1.22638499e-01j,
3.39625823e-02-1.16216469e+00j, 2.65140034e-06-4.10296457e-06j],
[ 4.17991622e-01+1.33504182e-02j, -4.79190633e-01-2.08562169e-01j,
4.87429517e-01+3.68070222e-01j, -3.63710538e-05+6.43912577e-06j],
[-2.18353842e+06-4.20344071e+05j, -2.52806647e+06-2.08794519e+05j,
-2.01808847e+06-1.96246695e+06j, -5.77147300e-01-3.12598394e+00j],
[-3.03044160e+05-6.45842521e+04j, -6.85879183e+05+2.07045473e+05j,
6.14194217e+04-1.28864668e+04j, -7.08794838e+00+9.70230041e+00j]])
U, s, Vh = scipy.linalg.svd(M)
print(s) # [4.42615634e+06 5.70600901e+05 4.68468171e-01 5.21600592e-13]
print(Vh[-1]) # [-5.35883825e-05+0.00000000e+00j 3.74712739e-05-9.89288566e-06j 4.03111556e-06+7.59306578e-06j -8.20834667e-01+5.71165865e-01j]
print(np.linalg.norm( M@Vh[-1] )) # 35.950705194666476
这是怎么回事? s[-1]很小,所以M@x原则上应该有解,但是Vh[-1]看起来不像解。这是 M 和 Vh 是复数的问题吗?数值 stability/accuracy 问题?还有别的吗?
我真的很想知道 x 会给 M@x 多少,其数量级与 s[-1] 大致相同,请告诉我任何解决此问题的方法。
你忘了共轭转置
SVD给出的分解为np.allclose(M, U @ np.diag(s) @ Vh)
,如果s[-1]
小则表示U @ np.diag(s) ~ M @ np.inv(Vh) ~ M @ Vh.T.conj()
的最后一列。所以你可以找到使用
M @ Vh[-1].T.conj() # [-7.77136331e-14-3.74441041e-13j,
# 4.67810503e-14+3.45797987e-13j,
# -2.84217094e-14-1.06581410e-14j,
# 7.10542736e-15+3.10862447e-15j]
我试图在 scipy 中使用 SVD 找到 M@x = 0
的近似非零解,其中 M 是一个复值 4x4 矩阵。
首先是一个玩具示例:
M = np.array([
[1,1,1,1],
[1,0,1,1],
[1,-1,0,0],
[0,0,0,1e-10]
])
U, s, Vh = scipy.linalg.svd(M)
print(s) # [2.57554368e+00 1.49380718e+00 3.67579714e-01 7.07106781e-11]
print(Vh[-1]) # [ 0.00000000e+00 2.77555756e-16 -7.07106781e-01 7.07106781e-01]
print(np.linalg.norm( M@Vh[-1] )) # 7.07106781193738e-11
所以在这种情况下,s中的最小(最后)值很小,对应的最后一列Vh[-1]是M@x=0
的近似解,其中M@Vh[-1]也很小,和 s[-1].
现在的真实例子不是以同样的方式工作:
M = np.array([[ 1.68572560e-01-3.98053448e-02j, 5.61165939e-01-1.22638499e-01j,
3.39625823e-02-1.16216469e+00j, 2.65140034e-06-4.10296457e-06j],
[ 4.17991622e-01+1.33504182e-02j, -4.79190633e-01-2.08562169e-01j,
4.87429517e-01+3.68070222e-01j, -3.63710538e-05+6.43912577e-06j],
[-2.18353842e+06-4.20344071e+05j, -2.52806647e+06-2.08794519e+05j,
-2.01808847e+06-1.96246695e+06j, -5.77147300e-01-3.12598394e+00j],
[-3.03044160e+05-6.45842521e+04j, -6.85879183e+05+2.07045473e+05j,
6.14194217e+04-1.28864668e+04j, -7.08794838e+00+9.70230041e+00j]])
U, s, Vh = scipy.linalg.svd(M)
print(s) # [4.42615634e+06 5.70600901e+05 4.68468171e-01 5.21600592e-13]
print(Vh[-1]) # [-5.35883825e-05+0.00000000e+00j 3.74712739e-05-9.89288566e-06j 4.03111556e-06+7.59306578e-06j -8.20834667e-01+5.71165865e-01j]
print(np.linalg.norm( M@Vh[-1] )) # 35.950705194666476
这是怎么回事? s[-1]很小,所以M@x原则上应该有解,但是Vh[-1]看起来不像解。这是 M 和 Vh 是复数的问题吗?数值 stability/accuracy 问题?还有别的吗?
我真的很想知道 x 会给 M@x 多少,其数量级与 s[-1] 大致相同,请告诉我任何解决此问题的方法。
你忘了共轭转置
SVD给出的分解为np.allclose(M, U @ np.diag(s) @ Vh)
,如果s[-1]
小则表示U @ np.diag(s) ~ M @ np.inv(Vh) ~ M @ Vh.T.conj()
的最后一列。所以你可以找到使用
M @ Vh[-1].T.conj() # [-7.77136331e-14-3.74441041e-13j,
# 4.67810503e-14+3.45797987e-13j,
# -2.84217094e-14-1.06581410e-14j,
# 7.10542736e-15+3.10862447e-15j]