Reinterpreting a matrix view?

General Usage
1

Hi all. I have a use case where I make a large memory buffer for a matrix, but I will regularly want to take a slightly smaller sub-array of it. Views are great, particularly for version 1.5+, but I’m wondering what happens when I pass my view of a matrix into LAPACK because it isn’t just a simple contiguous chunk of memory anymore. As an example, consider this:

_A = randn(36, 36)
A  = _A'_A
B = view(A, 1:24, 1:24)
_C = randn(24, 24)
C = _C'_C
cholesky!(B) 
cholesky!(C) # won't this be faster because of the memory layout?

I’m having a little trouble figuring out the smart way to benchmark the in-place Cholesky, so apologies if there’s some easy way to see that the factorizations for B and C cost the same. But I would think they don’t, particularly at non-toy sizes.

With this issue in mind, I’m wondering if there’s some way to do the following thing: take the first 24*24 entries of the buffer for A and then treat that as a native 24x24 matrix, disregarding the structure about A that the view operation maintains. Something like reinterpret(Matrix{Float64}, (24, 24), view(A, 1:24*24)) or something, which I recognize is not anywhere near a valid function call. Is there some tidy way to do that?

Am I overthinking this? Is this not a problem? I’d appreciate any thoughts you have to share. Thanks in advance for reading.

2

I think the impact is minimal. views keep track of the fact that within your submatrix, the entries within a column are contiguous (the important part for speed), even if the end on one column is not adjacent to the start of the next (the relatively unimportant part). You need to be slightly more clever about design choices if you’re slicing a higher dimensional array.

4 Likes
3

Use

using BenchmarkTools
@btime cholesky!(C) setup=(C=copy(A)) evals=1

I think there will be no miracle there. You will probably have to try using slices (that copy, and then copy back) and views and see what is faster at the end for your particular case. If the elements of the slice are not spread on the original matrix I guess the views will be the best choice.

Edit: To evaluate the time of the view option you will need to do:

@btime cholesky!(B) setup=(C=copy(A); B=view(C,1:24,1:24)) evals=1
1 Like
4

LAPACK can still handle it as long as your submatrix is contiguous within each column, i.e. unit-stride ranges in the first dimension.

2 Likes
5

Thanks to all of you for the helpful tips and guidance! I think I’ll just not play with this micro-optimization after seeing your comments.