4

while reading the scipy lecture : http://scipy-lectures.github.io/intro/numpy/operations.html

there is an example: 

>>> a = np.triu(np.ones((3, 3)), 1)
>>> a
array([[ 0.,  1.,  1.],
       [ 0.,  0.,  1.],
       [ 0.,  0.,  0.]])
>>> a.T
array([[ 0.,  0.,  0.],
       [ 1.,  0.,  0.],
       [ 1.,  1.,  0.]])

and then it says: enter image description here

and I don't understand why. I ran an experiment, and it does make it symmetric. enter image description here

EDIT 1 :

The same results holds for random matrix: enter image description here

fast tooth
  • 2,317
  • 4
  • 25
  • 34

1 Answers1

5

I think I understand what the purpose of the warning was, though I don't know enough about Numpy internals to know why the issue isn't happening.

The point is that the transposed matrix that you have on the right side of the in-place addition operation is a shallow view of the matrix that is being modified (not a copy). If Numpy were to perform the operation by iterating over the values, some results would end up incorrect.

Consider this (crude) implementation of something like the += operator for 2d matrices:

def matrix_iadd(lhs, rhs):
    for i in range(lhs.shape[0]):
        for j in range(lhs.shape[1]):
            lhs[i,j] += rhs[i,j]
    return lhs

This will not work as expected if rhs is a view that overlaps with lhs (like a transpose). Here's an example:

>>> a = np.arange(9)
>>> a.shape=(3,3)
>>> a
array([[0, 1, 2],
       [3, 4, 5],
       [6, 7, 8]])
>>> a.T
array([[0, 3, 6],
       [1, 4, 7],
       [2, 5, 8]])
>>> matrix_iadd(a, a.T)
array([[ 0,  4,  8],
       [ 7,  8, 12],
       [14, 19, 16]])

The reason the results are not symmetric is that the right hand side was changing in the middle of the operation, as it was going on. When lhs[0,1] got updated (with rhs[0,1]'s 3 being added to the existing value of 1), rhs[1,0] changed with it. So, when lhs[1,0]'s turn came to be updated, rhs[1,0] held a modified value.

I suspect the warning is about this, though the effect seems not to apply to numpy's builtin operators in normal useage. Perhaps the guide and its warning were written for a previous version of Numpy that had less helpful behavior, or perhaps the bad behavior can still show up with some kinds of arrays even in the current version of Numpy. I'm not sure. I suspect it is good advice in general, to avoid in-place modifications using views, even if it is not guaranteed to break in some uses.

Blckknght
  • 100,903
  • 11
  • 120
  • 169
  • is "a = a + a.T" safe? is python going to compute and create an object for the RHS and bind it with the name 'a'? – fast tooth Jun 09 '14 at 15:41
  • As someone who actually *had* this problem and came to stackoverflow to solve it (http://stackoverflow.com/questions/17091285/semantics-of-generating-symmetric-matrices-in-numpy), I am confident in saying that it is an issue that is not just simply solved by recent updates to the numpy internal `__iadd__`. – aestrivex Jun 09 '14 at 15:41
  • 1
    @fast tooth, yes, `a=a+a.T` is safe. – aestrivex Jun 09 '14 at 15:43
  • it seems to me that when the matrix becomes large, the error starts to show. – fast tooth Jun 09 '14 at 15:44
  • 1
    If it's working is because numpy is doing [buffering](http://docs.scipy.org/doc/numpy/reference/ufuncs.html#use-of-internal-buffers). It is an implementation detail that you should not rely on: it will fail for arrays larger than the buffer size, and is not guaranteed to keep working on future releases. If @seberg is around, he can probably give the proper lecture of exactly when does it currently work. – Jaime Jun 09 '14 at 17:15
  • 1
    That's the correct answer: it's about buffering, and will fail for larger arrays. – Gael Varoquaux Oct 24 '15 at 22:20