invenia / PDMatsExtras.jl

Extra Positive (Semi-)Definite Matricies
MIT License
8 stars 6 forks source link

`string(::WoodburyPDMat)` has worse than cubic complexity #30

Closed tpgillam closed 2 years ago

tpgillam commented 2 years ago

Consider the following little example that converts a n * n woodbury matrix to a string. This operation should be quadratic in n, both in terms of time and space.

function f(n::Int64)
    x = WoodburyPDMat(ones(n, 5), Diagonal(ones(5)), Diagonal(ones(n)))
    string(x)
    return nothing
end

A little benchmarking...

julia> using BenchmarkTools
julia> @btime f(10)
  240.060 μs (1314 allocations: 511.56 KiB)

julia> @btime f(20)
  2.219 ms (5215 allocations: 6.01 MiB)

julia> @btime f(30)
  8.254 ms (11716 allocations: 28.16 MiB)

julia> @btime f(40)
  22.650 ms (20817 allocations: 85.41 MiB)

n: 10 -> 20 => 9x increase in time n: 20 -> 40 => 10x increase in time

... which is more than O(n^3).

Note that if we convert to a standard matrix first, things seem happier:

function g(n::Int64)
    x = WoodburyPDMat(ones(n, 5), Diagonal(ones(5)), Diagonal(ones(n)))
    string(Matrix(x))  # NB!
    return nothing
end
julia> @btime g(10)
  32.008 μs (318 allocations: 48.23 KiB)

julia> @btime g(20)
  124.851 μs (1219 allocations: 187.16 KiB)

julia> @btime g(30)
  275.677 μs (2720 allocations: 418.97 KiB)

julia> @btime g(40)
  489.239 μs (4821 allocations: 741.66 KiB)

n: 10 -> 20 => 4x increase in time n: 20 -> 40 => 4x increase in time

... which is nearly bang on O(n^2).

Needless to say, for moderately large n, this becomes most unpleasant:

julia> @btime g(400)
  49.998 ms (480029 allocations: 71.96 MiB)

julia> @btime f(400)
  136.658 s (2720023 allocations: 767.87 GiB)

I think this bad scaling is fundamentally due to the iteration protocol not being (efficiently) implemented, although I haven't looked into this in too much detail, beyond noting that Array(::WoodburyPDMat) also seems to suffer from something similar.