Default `Diagonal` in Julia doesn't support sparse subarray: - https://github.com/JuliaArrays/FillArrays.jl/issues/265#issuecomment-2285492169 - https://github.com/JuliaApproximation/ClassicalOrthogonalPolynomials.jl/issues/202
