Running two subsequent reconstructions with system matrix caching causes memory issues

[markus-jehl]

When I run two successive reconstructions with system matrix caching (recomputing the system matrix in both), then the second reconstruction gets killed while computing the system matrix. This new behaviour was introduced with PR 1592 (#1592). Sanitisers are not flagging any memory leaks, so the issue is most likely caused by allocating huge chunks of contiguous memory, where the second time around the memory may already be fragmented.

If fragmentation is indeed the problem, the following may be solutions:

• do not insist on contiguous memory for system matrices (which are huge: in my case around 190 GB)
• find a way to not just free the allocated memory, but actively pass it back to the operating system (apparently "malloc_trim(0)" could work for linux)
• another option is apparently to request the memory directly from the operating system, rather than with "new". On linux this would be with mmap and on Windows with VirtualAlloc. Freeing (munmap) this would immediately return the memory to the operating system.

[KrisThielemans]

I'm not sure why this is happening. In fact, the matrix is not allocated as 1 huge chunk.

STIR/src/recon_buildblock/ProjMatrixByBin.cxx

Lines 213 to 214 in c6ee75a

	cache_collection[bin.view_num()][bin.segment_num()].insert(
	MapProjMatrixElemsForOneBin::value_type(cache_key(bin), probabilities));

Definition of cache_collection

STIR/src/include/stir/recon_buildblock/ProjMatrixByBin.h

Line 212 in c6ee75a

mutable VectorWithOffset<VectorWithOffset<MapProjMatrixElemsForOneBin>> cache_collection;

I don't think #1592 actually modified the behaviour of VectorWithOffset (it added extra capability), but I might be wrong.

I suppose we could try to use std::vector for cache_collection. I'd just mean using view_num() - min_view_num() etc as indices for it (and the locks, I guess). That would be a good move anyway (reducing our dependency on our own classes as much as possible). But I'm not sure if this will remove the problem really.

It could also be the opposite problem: too fragmented.

Still, I don't see how #1592 modified the projection matrix, so I'm at a loss...

[markus-jehl]

Just adding "malloc_trim(0);" between my two reconstructions solved the issue for me. So I think it is really that the freed memory is still associated with the process and then fragmented by something between the two reconstructions. By specifically handing back all memory to the operating system it is possible to re-allocate 192 GB the second time.

[markus-jehl]

I'll trace the changes in the PR a bit further down to understand where it impacted the caching. Do you see a reason for having the system matrix in contiguous memory? Otherwise I think it would be nice to turn it off there, since that cache is usually huge.

[KrisThielemans]

Do you see a reason for having the system matrix in contiguous memory?

but it isn't

[markus-jehl]

Do you see a reason for having the system matrix in contiguous memory?

but it isn't

Hmm... true. Now I got there, my theory doesn't make any sense :-/ So could it be that some other array allocations (e.g. projdata) have been building up over time, causing an issue when I try to allocate num_view * num_segment ProjMatrices?

[markus-jehl]

I think you may be right, that the memory is too fragmented. The cache is a 2D Array containing unordered_maps containing ProjMatrixElemsForOneBin, which are short vectors. Especially, since the unordered_maps are not pre-allocated and will therefore be grown and rehashed several times (this could also be what is slowing down the first reconstruction in my setup with STIR in Docker on WSL2). Is there a reason to use a hash map for the cache? Shouldn't we know beforehand how many ProjMatrixElemsForOneBin we need, since we know how many bins we have in the ProjData?

[KrisThielemans]

I guess the hash-map was a lazy way around not having to think about it. Also, if you do list-mode recon (and have sensitivity), you might not need all of them.

However, knowing how many you need might not solve the problem, as the elements are just ProjMatrixElemsForOneBin, which in itself is allocated.

STIR/src/include/stir/recon_buildblock/ProjMatrixElemsForOneBin.h

Line 192 in 840f911

std::vector<value_type> elements;

The question really is how "large" each LOR is, which is very hard to know, without doing calculations first (which is something that the SPECTUB projector does do).

So, I'm afraid I have no good solution. I suppose we could do a malloc_trim (within some '#ifdef) within the ProjMatrixByBin destructor, although you hope that wouldn't be necessary. Much more work would be to do "placement new` and take care of all this ourselves (but I'm not looking forward to that)

[markus-jehl]

I agree, there seems to be no perfect solution. I've seen many forums, though, that say that std::unordered_map is very slow, especially when not pre-allocated - replacing this might at least speed up the system matrix allocation in my setup with docker on WSL2, and maybe reduce the fragmentation a tiny bit. I'll see if there is a moderately painful way of doing this!

Never heard of "placemement new", but this could also be interesting! Otherwise malloc_trim in the destructor is definitely the easiest, but wouldn't solve my parallel problem of the slow initial allocation.

[KrisThielemans]

I'd be fine with replace unordered_map with vector<vector<>> or a 2d version (STIR or std). the default-constructed ProjMatrixElemsForOneBin should be quite small anyway.

If you think the growing a map reallocates/copies (not sure why it would, but I haven't looked into it. it appears it could do a rehash, but I'd guess that's relatively fast compared to the rest of the calculations), then we can very easily call reserve presumably.

As usual, some tests with valgrind/cachegrind seem in order (on a smaller system!). I don't have time for this for at least another month though.

I think we will have to release 6.3 with code as-is.