Умножение матрицы Numpy U*B*UT приводит к несимметричной матрице

Question

В моей программе мне нужно следующее умножение матриц:

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 заставить его быть симметричным. Это хороший способ обойти эту проблему?

Answer 1

Вы заметили, что So U, a 10*5 matrix, is indeed orthonormal except numerical rounding causes not exactly identity.

То же самое относится и к 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]])

я использую np.allclose при сравнении массивов поплавков.

Я тоже предпочитаю ndarray а также np.dot над np.matrix потому что поэлементное умножение так же часто, как матричное умножение.

Если остальная часть кода зависит от 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