Use a native list of numpy arrays, then np.concatenate.
The native list will multiply (by ~1.125) in size when needed, so not too many reallocations will occur, moreover, it will only hold pointers to scattered (non contiguous in memory) np.arrays holding the actual data.
Calling concatenate only once will solve your problem.
Pseudocode
dataset = []
for f in glob.glob('*.png'):
x = read_as_numpyarray(f) # custom function; x is a matrix of shape (n, 1000)
dataset.append(x)
dataset_np = np.concatenate(dataset)
Notice vstack internally uses concatenate.
Edit to address the edited question:
Let's say the total size of data is 20 GB. When concatenating, the
system will have to still keep 20 GB (for each individual array) and
also allocate 20 GB for the new concatenated array, thus requiring 40
GB of RAM (double of the dataset). How to do this without requiring
the double of RAM? (Example: is there a solution if we only have 32 GB
of RAM?)
I would attack this problem first by doing the same as proposed in this current answer, in phases. dataset_np1 = np.concatenate(half1_of_data), dataset_np2 = np.concatenate(half2_of_data), will only need 150% RAM (not 200%). This can be extended recursively at the expense of speed until the limit at which this becomes the proposition in the question. I can only assume the likes of dask can handle this better, but haven't tested myself.
Just to clarify, after you have dataset_np1 you no longer need the list of all the sharded small arrays, and can free that. Only then you start loading the other half. Thus, you only ever need to hold an extra 50% of the data in memory.
Pseudocode:
def load_half(buffer: np.array, shard_path: str, shard_ind: int):
half_dataset = []
for f in glob.glob(f'{shard_path}/*.png'):
x = read_as_numpyarray(f) # custom function; x is a matrix of shape (n, 1000)
half_dataset.append(x)
half_dataset_np = np.concatenate(half_dataset) # see comment *
buffer[:buffer.shape[0] // 2 * (shard_ind + 1), ...] = half_dataset_np
half1_path = r"half1" # preprocess the shards to be found by glob or otherwise
half2_path = r"half2"
assert os.path.isdir(half1_path)
assert os.path.isdir(half2_path)
buffer = np.zeros(size_shape)
half1_np = load_half(half1_path, buffer, 0) # only 50% of data temporarily loaded, then freed [can be done manually if needed]
half2_np = load_half(half2_path, buffer, 1) # only 50% of data temporarily loaded, then freed
One could (easily, or not so easily) generalize this to quarters, eighths, or recursively any required fraction to reduce memory costs at the expense of speed, with the limit at infinity being the original proposition in the question.
- Important comment (see "see comment * in the code):
One might notice half_dataset_np = np.concatenate(half_dataset)
actually allocates 50% of the dataset, with the other 50% allocated
in shards, apparently saving us nothing. That is correct, and I could
not find a way to concat into a buffer. However, implementing this
recursively as suggested (and not shown in pseudocode) will save
memory, as a quarter will only use 2* 25% every time. This is just an
implementation detail, but I hope the gist is clear.
On a different note, another approach would state "what if the dataset is 1000GB"? then no numpy array will do. This is why databases exist, and they can be queried quite efficiently using tools. But again, this is somewhat a research question, and depends heavily on your specific needs. As a very uninformed hunch, I would check out dask.
Such libraries will obviously tackle problems like this one as a subset of what they do, and I recommend not implementing these things yourself, as the total time you will spend will much outweigh the time choosing and learning a library.
On another different note, I wonder if this really has to be such a huge array, and maybe a slightly different design or formulation of the problem could alleviate us from this technical issue altogether.
