# Creating one's own Matrix type

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

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

When you say

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

``````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)
``````

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

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

Thanks for the pointer!