Why is assignment operation still giving MemoryError for large arrays in Python?

Question

I have a large array K (29000 x 29000):

K= numpy.random.random((29000, 29000))

I want to apply the following operation on K:

output = K* (1.5 - 0.5 * K* K)

To try preventing 'MemoryError' , I am doing my computations as suggested on the answer from this thread.

However, when I try to do the assignment operation on the large array as follows, I still get the MemoryError:

K *= 1.5 - 0.5 * K * K

Any help welcome.

NOTE: this is not a duplicate post. There is a suggestion on this post using cython. But I am looking for alternative solutions which may not rely on Cython.

Use `dask.array`. chek it here http://dask.pydata.org/en/latest/array.html — Sahil Dahiya, Jan 05 '18 at 15:35
Possible duplicate of [How can I apply the assignment operator correctly in Python?](https://stackoverflow.com/questions/48115074/how-can-i-apply-the-assignment-operator-correctly-in-python) — Till Hoffmann, Jan 05 '18 at 15:40
@Till Hoffman not a duplicate.......this post is made separately to raise the issue of MemoryError to the wider audience. Your suggestion using cython is appreciated but the aim of this new post is to make clear my current issue of MemoryError as I am looking for alternative solutions which may not rely on Cython. — user121, Jan 05 '18 at 15:49
The expression `1.5 - 0.5 * K * K` still requires the creation of temporary arrays to hold intermediate results. Your array requires almost 7 gigabytes. How much RAM does your computer have? — Warren Weckesser, Jan 05 '18 at 16:00
@WarrenWeckesser 16GB of RAM...can the assignment operation be done incrementally for the given expression `1.5 - 0.5 * K * K` ? — user121, Jan 05 '18 at 16:08
Perhaps `temp = K*K; temp *= -0.5; temp += 1.5; K *= temp; del temp`. This should avoid ever having to have 3 arrays in memory. — Steven Rumbalski, Jan 05 '18 at 16:13
unknown121, @StevenRumbalski just gave the sequence that I was going to suggest, but it still requires one more array of the same size as K. My other suggestion is to work in batches, which is what user `if....` just suggested in an answer. — Warren Weckesser, Jan 05 '18 at 16:25

score 5 · Accepted Answer · 2018-01-05T17:59:43.890

5

You can do assignment in blocks, say, of 1000 rows. The additional array this creates will be 1/29 of the size of your array, and having a for loop running 29 times shouldn't be much of a speed problem. Typical memory/speed tradeoff.

block = 1000          # the size of row blocks to use 
K = np.random.random((29000, 29000))
for i in range(int(np.ceil(K.shape[0] / block))):
    K[i*block:(i+1)*block, :] *= 1.5 - 0.5 * K[i*block:(i+1)*block, :]**2

Since there was some concern about the performance on smaller matrices, here is a test for those:

block = 1000
K = np.arange(9).astype(np.float).reshape((3, 3))
print(1.5 * K - 0.5 * K**3)
for i in range(int(np.ceil(K.shape[0] / block))):
    K[i*block:(i+1)*block_size, :] *= 1.5 - 0.5 * K[i*block:(i+1)*block_size, :]**2
print(K)

This prints

[[   0.    1.   -1.]
 [  -9.  -26.  -55.]
 [ -99. -161. -244.]]

twice.

edited Jan 05 '18 at 17:59

answered Jan 05 '18 at 16:20

is this superior to the answer by @StevenRumbalski in the comments above? – user121 Jan 05 '18 at 16:23
Test both and you'll find out. I don't have your memory/CPU configuration. – Jan 05 '18 at 16:23
is the suggested code robust to all array sizes (e.g., 30213 x 30213)? – user121 Jan 05 '18 at 16:30
Yes, because I round up with `ceil`, and because NumPy handles slicing that goes beyond the index bounds. – Jan 05 '18 at 16:38
hmmm.....if a 3x3 array is used as input, I do not seem to be getting a consistent output as the solution from @StevenRumbalski – user121 Jan 05 '18 at 16:42
1

I'd say this is better solution than @StevenRumbalski's, because in his approach he's iterating through all cells many times, where as here only as many as needed are iterated (basic algebra is almost free and even in CPU you can do many computations in parallel). Also, consider using `np.power(K, 2)` instead of `K**2`. – Dawid Laszuk Jan 05 '18 at 16:55
indeed, but how can it handle an array size of 3x3 for instance.....as I mention, it does not seem to produce the correct output. – user121 Jan 05 '18 at 17:03
@unknown121 I added a test with 3 by 3 to my answer. It works because that's how NumPy slicing works: when you ask for the first 1000 entries of an array with 3 entries, you get those 3 entries. – Jan 05 '18 at 18:00
@if.... great. For a large array (29000 x29000), if I wanted to apply the following: `C = linalg.cholesky(K, lower=True)`, I get MemoryError. Do you have any tip on how to handle this case? – user121 Jan 05 '18 at 18:07
1

@unknown121 This is a different question. Try casting the array to single precision (np.float32) if that is enough precision for your purpose. There is also [`memmap`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.memmap.html) that could be used to hold the array on disk inside of memory, freeing the memory for temporary arrays creating in the computation. – Jan 05 '18 at 18:47
@if....I posted the question [here](https://stackoverflow.com/questions/48119493/how-to-overcome-memoryerror-for-scipy-function). Would you be able to demonstrate how `memmap` can be utilized? – user121 Jan 05 '18 at 18:59

Why is assignment operation still giving MemoryError for large arrays in Python?

1 Answers1

Linked