How to avoid inserting unnecessary zeroes in a scipy sparse matrix
I ended up implementing the solution I proposed in my last paragraph.
For each 2nd diagonal, I construct a block sparse matrix without any padding, then I transform it into a matrix of the proper shape by a right-hand side product with a block matrix whose blocks are identity. I do need to store zeroes here to use bsr_matrix (I first gave a try to the method scipy.sparse.block_diag, but it was extremely slow), but less of them compared to my semi-padding solution: (4n+7)(n+1) vs (4n+8)(n+3); and they can be represented with 8 bits instead of 128. Execution time is increased by ~40% but I can live with that (and it's still a decrease of 20% compared to the first solution).
I might be missing something here, but for now I'm pretty satisfied with this solution.
EDIT
When trimming the zeroes of the RHS matrices before effecting the product, the execution time is reduced by 30% compared to the previously most efficient solution; all's well that ends well.
I'm toying with a problem modelled by a linear system, which can be written as a square block-tridiagonal matrix. These blocks are of size b = 4n+8, and the full matrix is of size Nb; N could be arbitrarily large (reasonably, of course) while n is kept rather small (typically less than 10).
The blocks themselves are sparse, the first diagonal being only identity matrices, and the second diagonals having only n+1 non-zero columns (so 3n+7 columns of zeroes) per block. These columns are contiguous, either zeroes then non-zeroes or the other way around.
Building all these blocks in memory results in a 3N-2 x b x b array that can be turned into a sparse matrix with scipy.sparse.bsr_matrix, then cast to CSR format and trimmed of the excess zeroes. It works nicely but I'd rather skip this temporary large and sparse array (for N = 1e4, n = 5 it's 5.6 zeros for every relevant entry!) altogether.

I had a look at scipy.sparse.dok_matrix, recommended for slicing and incremental building. Creating my entries fits in a tidy loop but the process is ~10 times longer than using bsr_matrix with my unnecessary dense array, which will be detrimental to the future use cases.
It doesn't seem like bsr_matrix can be used directly with scipy sparse matrices as input.
Using bsr_matrix without including the diagonal blocks, then adding a sparse eye greatly reduces the number of zeros (3.5 per relevant entry in my test configuration) and speeds up the process by a third compared to the original solution. Score!

Any clue on things that I could do to further reduce the original imprint of this matrix? The obvious goal being to give me more freedom with the choice of N.
EDIT
I managed to improved things a tad more by constructing the three block-diagonals separately. By doing so, I need less padding for my 2nd diagonals (n+3 instead of 3n+7; 1.3 zeroes per relevant entry), dividing my original blocks into two vertical blocks (one full of zeroes) and I only need it one diagonal at a time, cutting the memory cost in half on top of that. The main diagonal remains constructed with the eye method. The icing on the cake: a speed up of 25% compared to my 3rd bullet point, probably because separating the two 2nd diagonals saves some array reshaping operations needed before using bsr_matrix. Compared to the original method, for my (N, n) = (1e4, 5) test case it's ~20M zeroes saved when comparing matrices before trimming. At 128 bits each, it's a decent gain already!
The only possible improvement that I can picture now is building these diagonals separately, without any padding, then inserting columns of zeros (probably via  products with block-matrices of identities) and finally adding everything together.
I also read something about using a dict to update an empty dok_matrix, but in my case I think I would need to expand lists of indices, take their Cartesian product to construct the keys and each element of my blocks would need to be an individual value as one apparently cannot use slices as dictionary keys.