numpy矩阵乘法U * B * U.T导致非对称矩阵 [英] Numpy Matrix Multiplication U*B*U.T Results in Non-symmetric Matrix
问题描述
在我的程序中，我需要以下矩阵乘法:
In my program, I need the following matrix multiplication:
A = U * B * U^T
其中，B是M * M对称矩阵，而U是N * M矩阵，其列是正交的.因此，我希望A也是一个对称矩阵.
where B is an M * M symmetric matrix, and U is an N * M matrix where its columns are orthonormal. So I would expect A is also a symmetric matrix.
但是，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的矩阵，实际上确实是正交的，只是数值舍入导致不完全相同.
So U, a 10*5 matrix, is indeed orthonormal except numerical rounding causes not exactly identity.
# 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.
Besides, I checked np.all(U.T*U == (U.T*U).T) would be False.
这是我的A矩阵不对称的原因吗?换句话说，这是一个无法避免的数字问题吗?
Is this the reason that my A matrix is not symmetric? In other words, is this a numerical issue one cannot avoid?
在实践中，如何避免这种问题并获得对称矩阵A?
In practice, how would one avoid this kind of issue and get a symmetric matrix A?
我使用了技巧A = (A + A.T)/2强制将其对称.这是解决这个问题的好方法吗?
I used the trick A = (A + A.T)/2 to force it to be symmetric. Is this a good way to get around this problem?
推荐答案
您观察到So U, a 10*5 matrix, is indeed orthonormal except numerical rounding causes not exactly identity.
A同样的道理-除数字舍入外，它是对称的:
The same reasoning applies to A - it is symmetric except for numerical rounding:
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]])
在比较浮点数组时，我使用np.allclose.
I use np.allclose when comparing float arrays.
与np.matrix相比，我还更喜欢ndarray和np.dot，因为逐元素乘法与矩阵乘法一样普遍.
I also prefer ndarray and np.dot over np.matrix because element by element multiplication is just as common as matrix multiplication.
如果其余代码取决于A是否对称，那么您的把戏可能是一个不错的选择.它在计算上并不昂贵.
If the rest of the code depends on A being symmtric, then your trick may be a good choice. It's not computationally expensive.
出于某些原因einsum避免了数值问题:
For some reason einsum avoids the numerical issues:
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
这篇关于numpy矩阵乘法U * B * U.T导致非对称矩阵的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
