在我的程序中,我需要以下矩阵乘法:
A = U * B * U^T
其中 B
是一个 M * M
对称矩阵,U
是一个 N * M
矩阵,其中它的列是正交的。所以我希望 A
也是一个对称矩阵。
但是,Python 并没有这么说。
import numpy as np
import pymc.gp.incomplete_chol as pyichol
np.random.seed(10)
# Create symmetric matrix B
b = np.matrix(np.random.randn(20).reshape((5,4)))
B = b * b.T
np.all(B== B.T)
而 B 确实是对称的:
In[37]: np.all(B== B.T)
Out[37]: True
# Create U
m = np.matrix(np.random.random(100).reshape(10,10))
M = m * m.T
# M
U, s, V = np.linalg.svd(M)
U = U[:, :5]
U.T * U
In[41]: U.T * U
Out[41]:
matrix([[ 1.00000000e+00, 0.00000000e+00, -2.77555756e-17,
-1.04083409e-17, -1.38777878e-17],
[ 0.00000000e+00, 1.00000000e+00, -5.13478149e-16,
-7.11236625e-17, 1.11022302e-16],
[ -2.77555756e-17, -5.13478149e-16, 1.00000000e+00,
-1.21430643e-16, -2.77555756e-16],
[ -1.04083409e-17, -7.11236625e-17, -1.21430643e-16,
1.00000000e+00, -3.53883589e-16],
[ 0.00000000e+00, 9.02056208e-17, -2.63677968e-16,
-3.22658567e-16, 1.00000000e+00]])
因此 U
,一个 10*5 矩阵,确实是标准正交的,除了数字舍入导致不完全相同。
# Construct A
A = U * B * U.T
np.all(A == A.T)
In[38]: np.all(A == A.T)
Out[38]: False
A
不是对称矩阵。
此外,我检查过 np.all(U.T*U == (U.T*U).T)
会是 False
。
这就是我的 A
矩阵不对称的原因吗?换句话说,这是一个无法回避的数值问题吗?
在实践中,如何避免此类问题并获得对称矩阵A
?
我使用技巧 A = (A + A.T)/2
来强制它是对称的。这是解决此问题的好方法吗?
最佳答案
您观察到 所以 U,一个 10*5 矩阵,确实是标准正交的,除了数字舍入导致不完全相同。
同样的推理适用于 A
- 除了数字舍入外,它是对称的:
In [176]: A=np.dot(U,np.dot(B,U.T))
In [177]: np.allclose(A,A.T)
Out[177]: True
In [178]: A-A.T
Out[178]:
array([[ 0.00000000e+00, -2.22044605e-16, 1.38777878e-16,
5.55111512e-17, -2.49800181e-16, 0.00000000e+00,
0.00000000e+00, -1.11022302e-16, -1.11022302e-16,
0.00000000e+00],
...
[ 0.00000000e+00, 0.00000000e+00, 1.11022302e-16,
2.77555756e-17, -1.11022302e-16, 4.44089210e-16,
-2.22044605e-16, -2.22044605e-16, 0.00000000e+00,
0.00000000e+00]])
我在比较 float 组时使用 np.allclose
。
与 np.matrix
相比,我也更喜欢 ndarray
和 np.dot
,因为逐个元素的乘法与矩阵乘法一样常见。
如果其余代码依赖于 A
是对称的,那么您的技巧可能是个不错的选择。它的计算成本并不高。
出于某种原因,einsum
避免了数值问题:
In [189]: A1=np.einsum('ij,jk,lk',U,B,U)
In [190]: A1-A1.T
Out[190]:
array([[ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]])
In [193]: np.allclose(A,A1)
Out[193]: True
关于python - Numpy 矩阵乘法 U*B*U.T 导致非对称矩阵,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/31415709/