MWE: Performance of views

[MikaelSlevinsky]

I'd really like to start using this library seriously, but certain elements (at least of my understanding) leave me perplexed. I can write two simple mul!-like functions that are faster than the built-in mul! for a diagonal matrix, but 1) simply recreating the blocks is fastest but should be unnecessary and 2) views, which do not recreate the banded-block-banded matrix, are slower and allocate. Here's the example:

EDIT: ok, the fastest way is to view the data blocks of the banded-block-banded matrix

julia> using BlockBandedMatrices

julia> function my_special_mul!(y, D, x)
           @assert size(D, 1) == size(D, 2) == length(x) == length(y)
           n = blocklength(x)
           for j in 1:n
               xb = x[Block(j)]
               DB = D[Block(j, j)]
               for i in 1:length(xb)
                   y.blocks[j][i] = DB[i, i]*xb[i]
               end
           end
           return y
       end
my_special_mul! (generic function with 1 method)

julia> function my_special_mul_with_a_view!(y, D, x)
           @assert size(D, 1) == size(D, 2) == length(x) == length(y)
           n = blocklength(x)
           for j in 1:n
               xb = x.blocks[j]
               DB = view(D, Block(j, j))
               for i in 1:length(xb)
                   y.blocks[j][i] = DB[i, i]*xb[i]
               end
           end
           return y
       end
my_special_mul_with_a_view! (generic function with 1 method)

julia> function my_special_mul_through_data!(y, D, x)
           @assert size(D, 1) == size(D, 2) == length(x) == length(y)
           n = blocklength(x)
           @inbounds for j in 1:n
               xb = view(x, Block(j))
               DB = view(D.data, Block(j))
               yb = view(y, Block(j))
               for i in 1:length(xb)
                   yb[i] = DB[i]*xb[i]
               end
           end
           return y
       end
my_special_mul_through_data! (generic function with 1 method)

julia> for n in 20:20:200
           D = BandedBlockBandedMatrix{Float64}(I, 1:n, 1:n, (0, 0), (0, 0))
           x = BlockVector(randn(sum(1:n)), 1:n)
           y = zero(x)
           @time my_special_mul_through_data!(y, D, x);
           @time my_special_mul!(y, D, x);
           @time my_special_mul_with_a_view!(y, D, x);
           @time mul!(y, D, x);
       end
  0.056795 seconds (38.60 k allocations: 2.499 MiB, 99.94% compilation time)
  0.032453 seconds (16.87 k allocations: 1.132 MiB, 99.91% compilation time)
  0.039557 seconds (39.35 k allocations: 2.721 MiB, 99.92% compilation time)
  0.000117 seconds (22 allocations: 1.344 KiB)
  0.000008 seconds
  0.000018 seconds (80 allocations: 17.656 KiB)
  0.000028 seconds (820 allocations: 51.250 KiB)
  0.000195 seconds (42 allocations: 2.594 KiB)
  0.000010 seconds
  0.000029 seconds (120 allocations: 36.375 KiB)
  0.000061 seconds (1.83 k allocations: 114.375 KiB)
  0.000403 seconds (62 allocations: 3.844 KiB)
  0.000015 seconds
  0.000041 seconds (160 allocations: 61.438 KiB)
  0.000093 seconds (3.24 k allocations: 202.500 KiB)
  0.000736 seconds (82 allocations: 5.094 KiB)
  0.000025 seconds
  0.000124 seconds (200 allocations: 93.312 KiB)
  0.000211 seconds (5.05 k allocations: 315.625 KiB)
  0.000701 seconds (102 allocations: 6.344 KiB)
  0.000031 seconds
  0.000078 seconds (240 allocations: 131.750 KiB)
  0.000186 seconds (7.26 k allocations: 453.750 KiB)
  0.000683 seconds (122 allocations: 7.594 KiB)
  0.000041 seconds
  0.000113 seconds (280 allocations: 176.125 KiB)
  0.000239 seconds (9.87 k allocations: 616.875 KiB)
  0.000671 seconds (142 allocations: 8.844 KiB)
  0.000052 seconds
  0.000126 seconds (320 allocations: 227.281 KiB)
  0.000312 seconds (12.88 k allocations: 805.000 KiB)
  0.000514 seconds (162 allocations: 10.094 KiB)
  0.000066 seconds
  0.000158 seconds (360 allocations: 284.938 KiB)
  0.000383 seconds (16.29 k allocations: 1018.125 KiB)
  0.001663 seconds (182 allocations: 11.344 KiB)
  0.000081 seconds
  0.000180 seconds (400 allocations: 349.719 KiB)
  0.000613 seconds (20.10 k allocations: 1.227 MiB)
  0.001204 seconds (202 allocations: 12.594 KiB)

julia> 

[dlfivefifty]

Here was my first attempt by stoping axes from allocating but it didn't change much:

julia> for n in 20:20:200
           ax = BlockArrays.BlockedOneTo(ArrayLayouts.RangeCumsum(Base.OneTo(n)))
           D = BandedBlockBandedMatrix{Float64}(I, (ax,ax), (0, 0), (0, 0))
           x = BlockVector(randn(sum(1:n)), (ax,))
           y = zero(x)
           @time my_special_mul_through_data!(y, D, x);
           @time my_special_mul!(y, D, x);
           @time my_special_mul_with_a_view!(y, D, x);
           @time mul!(y, D, x);
       end
  0.000004 seconds
  0.000019 seconds (80 allocations: 5.938 KiB)
  0.000013 seconds (210 allocations: 13.125 KiB)
  0.000074 seconds (22 allocations: 1.312 KiB)
  0.000005 seconds
  0.000011 seconds (160 allocations: 18.219 KiB)
  0.000054 seconds (820 allocations: 51.250 KiB)
  0.000041 seconds (42 allocations: 2.562 KiB)
  0.000007 seconds
  0.000041 seconds (240 allocations: 37.188 KiB)
  0.000059 seconds (1.83 k allocations: 114.375 KiB)
  0.000050 seconds (62 allocations: 3.812 KiB)
  0.000012 seconds
  0.000034 seconds (320 allocations: 62.438 KiB)
  0.000148 seconds (3.24 k allocations: 202.500 KiB)
  0.000078 seconds (82 allocations: 5.062 KiB)
  0.000017 seconds
  0.000036 seconds (400 allocations: 94.500 KiB)
  0.000158 seconds (5.05 k allocations: 315.625 KiB)
  0.000082 seconds (102 allocations: 6.312 KiB)
  0.000025 seconds
  0.000044 seconds (480 allocations: 133.469 KiB)
  0.000179 seconds (7.26 k allocations: 453.750 KiB)
  0.000108 seconds (122 allocations: 7.562 KiB)
  0.000032 seconds
  0.000093 seconds (560 allocations: 178.156 KiB)
  0.000209 seconds (9.87 k allocations: 616.875 KiB)
  0.000116 seconds (142 allocations: 8.812 KiB)
  0.000042 seconds
  0.000098 seconds (640 allocations: 229.531 KiB)
  0.000410 seconds (12.88 k allocations: 805.000 KiB)
  0.000141 seconds (162 allocations: 10.062 KiB)
  0.000052 seconds
  0.000110 seconds (720 allocations: 287.469 KiB)
  0.000339 seconds (16.29 k allocations: 1018.125 KiB)
  0.000162 seconds (182 allocations: 11.312 KiB)
  0.000065 seconds
  0.000112 seconds (800 allocations: 351.938 KiB)
  0.000429 seconds (20.10 k allocations: 1.227 MiB)
  0.000220 seconds (202 allocations: 12.562 KiB)

EDIT: using your updated version

[dlfivefifty]

Oh wait it seems to have made mul! much faster...

[dlfivefifty]

Oh my "fix" doesn't actually do anything since it changes the lazy axes to a Vector.... so I need to fix that

[dlfivefifty]

OK I fixed that, and the improved version is now both slower and uses more memory 😂

  0.000010 seconds (40 allocations: 1.250 KiB)
  0.000028 seconds (120 allocations: 7.188 KiB)
  0.000022 seconds (630 allocations: 26.250 KiB)
  0.000094 seconds (62 allocations: 2.562 KiB)
  0.000006 seconds (80 allocations: 2.500 KiB)
  0.000016 seconds (240 allocations: 20.719 KiB)
  0.000067 seconds (2.46 k allocations: 102.500 KiB)
  0.000029 seconds (122 allocations: 5.062 KiB)
  0.000009 seconds (120 allocations: 3.750 KiB)
  0.000033 seconds (360 allocations: 40.938 KiB)
  0.000178 seconds (5.49 k allocations: 228.750 KiB)
  0.000085 seconds (182 allocations: 7.562 KiB)
  0.000014 seconds (160 allocations: 5.000 KiB)
  0.000038 seconds (480 allocations: 67.438 KiB)
  0.000207 seconds (9.72 k allocations: 405.000 KiB)
  0.000101 seconds (242 allocations: 10.062 KiB)
  0.000020 seconds (200 allocations: 6.250 KiB)
  0.000050 seconds (600 allocations: 100.750 KiB)
  0.000515 seconds (15.15 k allocations: 631.250 KiB)
  0.000108 seconds (302 allocations: 12.562 KiB)
  0.000029 seconds (240 allocations: 7.500 KiB)
  0.000067 seconds (720 allocations: 140.969 KiB)
  0.000384 seconds (21.78 k allocations: 907.500 KiB)
  0.000119 seconds (362 allocations: 15.062 KiB)
  0.000036 seconds (280 allocations: 8.750 KiB)
  0.000063 seconds (840 allocations: 186.906 KiB)
  0.000581 seconds (29.61 k allocations: 1.205 MiB)
  0.000682 seconds (422 allocations: 17.562 KiB)
  0.000047 seconds (320 allocations: 10.000 KiB)
  0.000134 seconds (960 allocations: 239.531 KiB)
  0.000870 seconds (38.64 k allocations: 1.572 MiB)
  0.000208 seconds (482 allocations: 20.062 KiB)
  0.000055 seconds (360 allocations: 11.250 KiB)
  0.000146 seconds (1.08 k allocations: 298.719 KiB)
  0.001011 seconds (48.87 k allocations: 1.989 MiB)
  0.000321 seconds (542 allocations: 22.562 KiB)
  0.000072 seconds (400 allocations: 12.500 KiB)
  0.001355 seconds (1.20 k allocations: 364.438 KiB)
  0.001470 seconds (60.30 k allocations: 2.454 MiB)
  0.000292 seconds (602 allocations: 25.062 KiB)

[dlfivefifty]

There seems to be issues with views allocating. If we change view to return a banded matrix it will fix the allocations:

function my_special_mul_through_data2!(y, D, x)
    @assert size(D, 1) == size(D, 2) == length(x) == length(y)
    n = blocklength(x)
    @inbounds for j in 1:n
        xb = view(x, Block(j))
        DB = _BandedMatrix(view(D.data, Block(j)), j, 0, 0 )
        yb = view(y, Block(j))
        mul!(yb, DB, xb)
    end
    return y
end

for n in 20:20:200
    D = BandedBlockBandedMatrix{Float64}(I, 1:n, 1:n, (0, 0), (0, 0))
    x = BlockVector(randn(sum(1:n)), 1:n)
    y = zero(x)
    @time my_special_mul_through_data!(y, D, x);
    @time my_special_mul!(y, D, x);
    @time my_special_mul_with_a_view!(y, D, x);
    @time mul!(y, D, x); 
    @time my_special_mul_through_data2!(y, D, x);
end

0.000005 seconds
  0.000012 seconds (80 allocations: 5.938 KiB)
  0.000025 seconds (210 allocations: 13.125 KiB)
  0.000070 seconds (22 allocations: 1.312 KiB)
  0.066451 seconds (45.65 k allocations: 2.180 MiB, 99.91% compilation time)
  0.000005 seconds
  0.000013 seconds (160 allocations: 18.219 KiB)
  0.000025 seconds (820 allocations: 51.250 KiB)
  0.000032 seconds (42 allocations: 2.562 KiB)
  0.000031 seconds
  0.000007 seconds
  0.000028 seconds (240 allocations: 37.188 KiB)
  0.000052 seconds (1.83 k allocations: 114.375 KiB)
  0.000044 seconds (62 allocations: 3.812 KiB)
  0.000055 seconds
  0.000012 seconds
  0.000027 seconds (320 allocations: 62.438 KiB)
  0.000108 seconds (3.24 k allocations: 202.500 KiB)
  0.000082 seconds (82 allocations: 5.062 KiB)
  0.000076 seconds
  0.000093 seconds
  0.000046 seconds (400 allocations: 94.500 KiB)
  0.000186 seconds (5.05 k allocations: 315.625 KiB)
  0.000085 seconds (102 allocations: 6.312 KiB)
  0.000324 seconds
  0.000024 seconds
  0.000051 seconds (480 allocations: 133.469 KiB)
  0.000276 seconds (7.26 k allocations: 453.750 KiB)
  0.000167 seconds (122 allocations: 7.562 KiB)
  0.000167 seconds
  0.000032 seconds
  0.000093 seconds (560 allocations: 178.156 KiB)
  0.000259 seconds (9.87 k allocations: 616.875 KiB)
  0.000151 seconds (142 allocations: 8.812 KiB)
  0.000133 seconds
  0.000045 seconds
  0.000096 seconds (640 allocations: 229.531 KiB)
  0.004092 seconds (12.88 k allocations: 805.000 KiB, 90.75% gc time)
  0.000364 seconds (162 allocations: 10.062 KiB)
  0.000167 seconds
  0.000056 seconds
  0.000108 seconds (720 allocations: 287.469 KiB)
  0.000558 seconds (16.29 k allocations: 1018.125 KiB)
  0.000535 seconds (182 allocations: 11.312 KiB)
  0.000258 seconds
  0.000106 seconds
  0.000136 seconds (800 allocations: 351.938 KiB)
  0.000651 seconds (20.10 k allocations: 1.227 MiB)
  0.000501 seconds (202 allocations: 12.562 KiB)
  0.000213 seconds

I think to match the speed we might need to special case diagonal banded matrix multiplication?

[dlfivefifty]

Actually I think the speed is because OpenBLAS sucks for banded matrices. Try AppleAccelerate/MKL.

[MikaelSlevinsky]

So basically if I write a banded_view, it can reuse the same data and use mul! From banded matrices

[MikaelSlevinsky]

I guess a view doesn't even need to return a SubArray. From the docstring example:

julia> typeof(view(2:5, 2:3))
UnitRange{Int64}

[dlfivefifty]

That’s correct and BlockArray already just returns a Matrix from view

I made a PR to change the view to a BandedMatrix