Creating one's own Matrix type

General Usage
1

I have a question related to 13:45 in Alan Edelman’s recent award presentation:

In the following two code snippets, is there a deeper difference other than MyMatrix being more descriptive, i.e. the type tells me the type of the elements and that it stores those elements as a vector? Is the first version more Julian? If so, why? Are there performance differences? The struct definition of Version 2 seems a little cleaner?!

Version 1

struct MyMatrix{T, V<:AbstractVector{T}} <: AbstractMatrix{T}
    v::V
end

Base.size(A::MyMatrix) = length(A.v), length(A.v)
Base.getindex(A::MyMatrix, i, j) = A.v[i]*(i==j) + A.v[i]*A.v[j]

julia> A = MyMatrix([1,2,3])

3×3 MyMatrix{Int64,Array{Int64,1}}:
 1  0  0
 0  2  0
 0  0  3

Version 2

struct MyNewMatrix{T} <: AbstractMatrix{T}
    v::AbstractVector{T}
end

Base.size(A::MyNewMatrix) = length(A.v), length(A.v)
Base.getindex(A::MyNewMatrix, i, j) = A.v[i]*(i==j) + A.v[i]*A.v[j]

julia>B = MyNewMatrix([1,2,3])

3×3 MyNewMatrix{Int64}:
 1  0  0
 0  2  0
 0  0  3
1 Like
2

Version 2 is not concrete. v is of abstract type at the moment of creation and (in principle) it will be much slower than the first one; You can always fix this by fixing the actual concrete type, i.e. declaring v::Vector{T}, but that makes the approach less flexible (you can’t use it with e.g. sparse vectors, before the’re converted to dense ones)

2 Likes
3

When you say

do you refer to creating matrices or do you refer to doing computations with them?

4

As a follow up:

julia> using BenchmarkTools
julia> n = 10000000
julia> A = @btime MyMatrix(rand(n))
julia> B = @btime MyNewMatrix(rand(n))

  20.005 ms (3 allocations: 76.29 MiB)
  19.737 ms (3 allocations: 76.29 MiB)
5

you’re benchmarking mostly the rand function:

julia> let v = rand(n)
          @btime MyMatrix($v);
          @btime MyNewMatrix($v);
       end;
  5.882 ns (1 allocation: 16 bytes)
  5.854 ns (1 allocation: 16 bytes)

however using those matrices has significant penalty leading to

julia> let v = rand(10_000)
          @btime sum(MyMatrix($v));
          @btime sum(MyNewMatrix($v));
       end;
  207.626 ms (7 allocations: 144 bytes)
  28.082 s (984670013 allocations: 14.67 GiB)

Edit: simple indexing is 70x slower

let v = rand(100)
          @btime MyMatrix($v)[10,40]
          @btime MyNewMatrix($v)[10,40]
       end
  2.320 ns (0 allocations: 0 bytes)
  152.165 ns (7 allocations: 112 bytes)
4 Likes
6

See Performance Tips · The Julia Language and the following section. The short answer is that you should declare the type in such as way that the type system knows exactly what it’s going to get from A.v.

3 Likes
7

Thanks for the pointer!