2

In my program, I need the following matrix multiplication:

A = U * B * U^T

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.

However, Python doesn't say so.

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)

And B is indeed symmetric:

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

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 is not a symmetric matrix.

Besides, I checked np.all(U.T*U == (U.T*U).T) would be False.

Is this the reason that my A matrix is not symmetric? In other words, is this a numerical issue one cannot avoid?

In practice, how would one avoid this kind of issue and get a symmetric matrix A?

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?

ali_m
  • 71,714
  • 23
  • 223
  • 298
panc
  • 817
  • 2
  • 14
  • 30
  • `*` is not a matrix multiplication :) – cel Jul 14 '15 at 19:35
  • 3
    @cel: For the numpy `matrix` class, `*` *is* matrix multiplication. – Warren Weckesser Jul 14 '15 at 19:41
  • 3
    @cel I purposefully used `np.matrix` to covert all numpy arrays to matrix, so that I can simply use `*` for matrix multiplication. – panc Jul 14 '15 at 19:43
  • @WarrenWeckesser, ouch, you're right. – cel Jul 14 '15 at 19:44
  • 2
    @Pan Chao: This is normal 64 bit floating point error. What are you doing with `A` that requires it to be *exactly* symmetric? – Warren Weckesser Jul 14 '15 at 19:51
  • 3
    As a general rule, if you have stuff in your code which depends on `np.all(A == A.T)` being True and not simply on `np.isclose(A, A.T).all()` being True then you're asking for trouble. – DSM Jul 14 '15 at 19:56
  • @cel - http://docs.scipy.org/doc/numpy/reference/generated/numpy.matrix.html?highlight=numpy.matrix#numpy-matrix .... *".... It has certain special operators, such as `*` (matrix multiplication) and `**` (matrix power).*" – rayryeng Jul 14 '15 at 21:01
  • Use `np.allclose` to compare floating point arrays. – hpaulj Jul 14 '15 at 21:02

1 Answers1

6

You observed that So U, a 10*5 matrix, is indeed orthonormal except numerical rounding causes not exactly identity.

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

I use np.allclose when comparing float arrays.

I also prefer ndarray and np.dot over np.matrix because element by element multiplication is just as common as matrix multiplication.

If the rest of the code depends on A being symmtric, then your trick may be a good choice. It's not computationally expensive.

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
hpaulj
  • 221,503
  • 14
  • 230
  • 353
  • That is interesting! I wonder if `einsum` is accumulating intermediate results in quad-precision floats? – ali_m Jul 14 '15 at 21:41
  • 2
    The double `dot` is doing a sum of products twice. `einsum` is calculating all the products, and then doing the sum over the 2 dimensions. There's no intermediate sum to cause rounding errors. – hpaulj Jul 14 '15 at 22:15