[Numpy-discussion] Problem with numpy's array reference assignment?

[Bernhard Spinnler]

I have problems to get a piece of code to work with a new numpy/scipy version. The code essentially sets up a matrix Ryy and a vector Rya and solves the system of linear equations Ryy*c = Rya for c. Then it checks whether the resulting vector c satisfies the equation: Ryy*c must be equal to Rya.

While this code runs fine on
- python-2.7.5.p1, numpy-1.7.0, scipy-0.12.0.p0
- python-2.7.3, numpy-1.7.1, scipy-0.12.0

it fails on
- python 2.7.2, numpy 1.9.0.dev-fde3dee, scipy 0.14.0.dev-4938da3
with error:

	AssertionError: 
	Arrays are not almost equal to 6 decimals

	(mismatch 100.0%)
	 x: array([ 9.+0.j,  8.+0.j,  7.+0.j])
	 y: array([ 7.+0.j,  6.+0.j,  5.+0.j])

The piece of code is:
------------------------------------
import numpy as np
from numpy.testing import assert_array_almost_equal

lag0 = 5
N = 3
D = 2
corr_ya = np.array([0,1,2,3,4,5,6,7,8,9+0j])
corr_yy = np.array([0,1,2,3,4,5,4,3,2,1])

Rya = corr_ya[lag0+D:lag0+D-N:-1]
Ryy = np.zeros((N, N), complex)
for i in np.arange(N):
    Ryy[i,:] = corr_yy[lag0-i:lag0-i+N]
    
c = np.linalg.solve(Ryy, Rya)

# check result
Ryy_x_c = np.dot(Ryy, c)
print "Ryy*c =", repr(Ryy_x_c)
print "Rya =", repr(Rya)
assert_array_almost_equal(Ryy_x_c, Rya)
------------------------------------

I guess that it has something to do with numpy assigning arrays by reference and not copying them, since the problem goes away when vector Rya is copied before passing it to solve, i.e. doing Rya_copy = copy(Rya) and passing Rya_copy to solve.

The problem also does not occur when the initial array corr_ya is an integer array, e.g. by deleting the "+0j" from the last element of corr_ya above. (Normally, my initial arrays are complex. I just used artificially simplified arrays above to show the problem.) I could imagine that this is due to numpy copying the integers into a new complex array instead of creating a reference that is passed to solve.

Could it be that a bug has been introduced in recent numpy/scipy version? Or am I misunderstanding something?

Thanks,
	Bernhard