Skip to content

Speedup matrix multiplications #129

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 8 commits into
base: master
Choose a base branch
from
Open

Speedup matrix multiplications #129

Changes from 1 commit
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
2915481
fix tests
CarloLucibello Dec 4, 2020
2766105
faster matrix-fillmatrix multiplication
CarloLucibello Dec 4, 2020
e7dcb94
add hermitian to tests
CarloLucibello Dec 4, 2020
4f08ad6
cleanup
CarloLucibello Dec 4, 2020
5a999ba
bump version
CarloLucibello Dec 5, 2020
9ef396f
generalise to AbstractFill for consistency
dlfivefifty Dec 10, 2020
15a15f1
Add vector' tests
dlfivefifty Dec 10, 2020
8350786
Update runtests.jl
dlfivefifty Dec 10, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
84 changes: 70 additions & 14 deletions src/fillalgebra.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,3 +1,10 @@
const FillVector{F,A} = Fill{F,1,A}
const FillMatrix{F,A} = Fill{F,2,A}
const OnesVector{F,A} = Ones{F,1,A}
const OnesMatrix{F,A} = Ones{F,2,A}
const ZerosVector{F,A} = Zeros{F,1,A}
const ZerosMatrix{F,A} = Zeros{F,2,A}

## vec

vec(a::Ones{T}) where T = Ones{T}(length(a))
Expand Down Expand Up @@ -87,11 +94,22 @@ end
*(a::Zeros{<:Any,2}, b::Diagonal) = mult_zeros(a, b)
*(a::Diagonal, b::Zeros{<:Any,1}) = mult_zeros(a, b)
*(a::Diagonal, b::Zeros{<:Any,2}) = mult_zeros(a, b)
function *(a::Diagonal, b::AbstractFill{<:Any,2})

# Cannot unify following methods for Diagonal
# due to ambiguity with general array mult. with fill
function *(a::Diagonal, b::FillMatrix)
size(a,2) == size(b,1) || throw(DimensionMismatch("A has dimensions $(size(a)) but B has dimensions $(size(b))"))
a.diag .* b # use special broadcast
end
function *(a::FillMatrix, b::Diagonal)
size(a,2) == size(b,1) || throw(DimensionMismatch("A has dimensions $(size(a)) but B has dimensions $(size(b))"))
a .* permutedims(b.diag) # use special broadcast
end
function *(a::Diagonal, b::OnesMatrix)
size(a,2) == size(b,1) || throw(DimensionMismatch("A has dimensions $(size(a)) but B has dimensions $(size(b))"))
a.diag .* b # use special broadcast
end
function *(a::AbstractFill{<:Any,2}, b::Diagonal)
function *(a::OnesMatrix, b::Diagonal)
size(a,2) == size(b,1) || throw(DimensionMismatch("A has dimensions $(size(a)) but B has dimensions $(size(b))"))
a .* permutedims(b.diag) # use special broadcast
end
Expand All @@ -100,23 +118,61 @@ end
*(a::Transpose{T, <:StridedMatrix{T}}, b::Fill{T, 1}) where T = reshape(sum(parent(a); dims=1) .* b.value, size(parent(a), 2))
*(a::StridedMatrix{T}, b::Fill{T, 1}) where T = reshape(sum(a; dims=2) .* b.value, size(a, 1))

function *(a::Adjoint{T, <:StridedMatrix{T}}, b::Fill{T, 2}) where T
fB = similar(parent(a), size(b, 1), size(b, 2))
fill!(fB, b.value)
return a*fB
function *(x::AbstractMatrix, f::FillMatrix)
axes(x, 2) ≠ axes(f, 1) &&
throw(DimensionMismatch("Incompatible matrix multiplication dimensions"))
m = size(f, 2)
repeat(sum(x, dims=2) * f.value, 1, m)
end

function *(f::FillMatrix, x::AbstractMatrix)
axes(f, 2) ≠ axes(x, 1) &&
throw(DimensionMismatch("Incompatible matrix multiplication dimensions"))
m = size(f, 1)
repeat(sum(x, dims=1) * f.value, m, 1)
end

function *(a::Transpose{T, <:StridedMatrix{T}}, b::Fill{T, 2}) where T
fB = similar(parent(a), size(b, 1), size(b, 2))
fill!(fB, b.value)
return a*fB
function *(x::AbstractMatrix, f::OnesMatrix)
axes(x, 2) ≠ axes(f, 1) &&
throw(DimensionMismatch("Incompatible matrix multiplication dimensions"))
m = size(f, 2)
repeat(sum(x, dims=2) * one(eltype(f)), 1, m)
end

function *(a::StridedMatrix{T}, b::Fill{T, 2}) where T
fB = similar(a, size(b, 1), size(b, 2))
fill!(fB, b.value)
return a*fB
function *(f::OnesMatrix, x::AbstractMatrix)
axes(f, 2) ≠ axes(x, 1) &&
throw(DimensionMismatch("Incompatible matrix multiplication dimensions"))
m = size(f, 1)
repeat(sum(x, dims=1) * one(eltype(f)), m, 1)
end

*(x::FillMatrix, y::FillMatrix) = mult_fill(x, y)
*(x::FillMatrix, y::OnesMatrix) = mult_fill(x, y)
*(x::OnesMatrix, y::FillMatrix) = mult_fill(x, y)
*(x::OnesMatrix, y::OnesMatrix) = mult_fill(x, y)
*(x::ZerosMatrix, y::OnesMatrix) = mult_zeros(x, y)
*(x::ZerosMatrix, y::FillMatrix) = mult_zeros(x, y)
*(x::FillMatrix, y::ZerosMatrix) = mult_zeros(x, y)
*(x::OnesMatrix, y::ZerosMatrix) = mult_zeros(x, y)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think you are massively underestimating the number of ambiguities we need to take into account. For example we need to take into account AdjointAbsVec, TransposeAbsVec, AbstractTriangular, Adjoint, Transpose, .... A somewhat more exhaustive list is here:

https://github.com/JuliaMatrices/ArrayLayouts.jl/blob/9fbed3a8e49b26d4ccff5608635b8c986f702b7a/src/mul.jl#L158

It's also a bad idea to overload Ones and Fills seperatly since that just means more ambiguities downstream. The following is better:

*(a::AbstractMatrix, b::AbstractFill{<:Any,2}) = fill_mul(a, b)
fill_mul(a, b::Ones) = ...
fill_mul(a, b::Zeros) = ...

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

hopefully I addressed both concerns

Copy link
Contributor

@mcabbott mcabbott Dec 5, 2020 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there a chance that overloading something other than * might be simpler, and avoid all ambiguities? They are pretty thorny, as every other package that tries to play the same game introduces new possibilities...

Without special handling, it will end up after similar at some generic_mul!(C, A::Fill, B). Could one overload that instead?

Or maybe that's too late, since similar(::Filll)::Array, but is there somewhere in the chain between * and BLAS or whatever at which this could be caught?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This issue is precisely what ArrayLayouts.jl is designed to avoid, and would also give us the mul! variants. One option would be to make FillArrays.jl depend on ArrayLayouts.jl (it’s currently the other way around).

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Wait how does that work? Is this only if you call something special like ArrayLayouts.mul or does it stick its fingers into ordinary A * B?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We would do AbstractFill <: LayoutArray which then catches all the necessary * ovverrides

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK, got it. Overloading for one supertype stops all its children from stepping on each others' toes.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

although that's probably the right direction, It involves a lot of churn and the coordination of two packages, so I would merge this as it is for the time being

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes I agree.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

the coordination of two packages

FYI I wrote both packages so this isn't quite so bad. The only reason I haven't done this is FillArrays.jl is clean and simple, so needs a good reason to make it more complicated. I'm not sure we are there yet.


# function *(a::Adjoint{T, <:StridedMatrix{T}}, b::Fill{T, 2}) where T
# fB = similar(parent(a), size(b, 1), size(b, 2))
# fill!(fB, b.value)
# return a*fB
# end

# function *(a::Transpose{T, <:StridedMatrix{T}}, b::Fill{T, 2}) where T
# fB = similar(parent(a), size(b, 1), size(b, 2))
# fill!(fB, b.value)
# return a*fB
# end

# function *(a::StridedMatrix{T}, b::Fill{T, 2}) where T
# fB = similar(a, size(b, 1), size(b, 2))
# fill!(fB, b.value)
# return a*fB
# end

function _adjvec_mul_zeros(a::Adjoint{T}, b::Zeros{S, 1}) where {T, S}
la, lb = length(a), length(b)
if la ≠ lb
Expand Down