Describe the workflow you want to enable

As discussed in #27359, make_sparse_spd_matrix actually returns a dense numpy array (with many zero values).

I think it should be possible to rewrite this code to compose operations on sparse arrays/matrices from the start instead of allocating large dense square matrices.

I have not tried myself, but I think it should be doable and it should be much more memory efficient.