No meu programa, preciso da seguinte multiplicação de matrizes:

A = U * B * U^T

OndeB é umM * M matriz simétrica eU é umN * M matriz onde suas colunas são ortonormais. Então eu esperariaA também é uma matriz simétrica.

No entanto, o Python não diz isso.

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)

E B é de fato simétrico:

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]])

assimU, uma matriz 10 * 5, é realmente ortonormal, exceto que o arredondamento numérico não causa exatamente a identidade.

# Construct A
A = U * B * U.T
np.all(A == A.T)

In[38]: np.all(A == A.T)
Out[38]: False

A não é uma matriz simétrica.

Além disso, eu verifiqueinp.all(U.T*U == (U.T*U).T) seriaFalse.

É por isso que meuA matriz não é simétrica? Em outras palavras, esse é um problema numérico que não se pode evitar?

Na prática, como evitar esse tipo de problema e obter uma matriz simétricaA?

Eu usei o truqueA = (A + A.T)/2 para forçá-lo a ser simétrico. Essa é uma boa maneira de contornar esse problema?