generic matmul depends on zero

[StefanKarpinski]

Currently generic matmul uses zero(eltype(a)) to fill the initial output array with zero values. When defining a new numeric type, this can be annoying (you may not have defined a zero method yet), but there are also real-world case where this causes matmul to fail, such as when you have an Any array that just happens to contain numbers that you want to multiply:

julia> A = convert(Array{Any}, rand(3,4))
3x4 Array{Any,2}:
 0.994922  0.18993   0.826601  0.82501
 0.10378   0.937165  0.753685  0.321605
 0.269249  0.678925  0.749532  0.305874

julia> A*rand(4)
ERROR: MethodError: `zero` has no method matching zero(::Type{Any})
 in generic_matvecmul! at linalg/matmul.jl:319
 in * at linalg/matmul.jl:74

On the one hand, you really want this to be a Float64 matrix for performance. But on the other hand, this really ought to work in any case.

[jiahao]

Similar in spirit: JuliaLang/julia#8319 JuliaLang/julia#9276

[timholy]

See if you like a potential fix in JuliaLang/julia#10622.

[kshyatt]

@andreasnoack did this ever get resolved?

[tkelman]

Still needed. Bit of a confusing gotcha in every demo that tries to show generic linear algebra capabilities.

[pabloferz]

At least the OP example works now on master.

[andreasnoack]

This doesn't error anymore.

[tkelman]

It does with any user defined type though. Whatever workaround was added for Any is apparently not always applicable?

[tkelman]

And where is it tested?

[andreasnoack]

I added a test in JuliaLang/julia@351e8af but I guess that doesn't cover the demo example you tried? Please share the example if that is the case.

[tkelman]

julia> immutable DNum{T}
           v::T
           e::T
       end
julia> Base.:+(x::DNum, y::DNum) = DNum(x.v + y.v, x.e + y.e)
julia> Base.:*(x::DNum, y::DNum) = DNum(x.v * y.v, x.v * y.e + y.v * x.e)
julia> DNum.(rand(1:10, 6, 4), rand(1:10, 6, 4)) * DNum.(rand(1:10, 4, 3), rand(1:10, 4, 3))
ERROR: MethodError: no method matching zero(::DNum{Int64})
Closest candidates are:
  zero(::Type{Base.LibGit2.GitHash}) at libgit2\oid.jl:88
  zero(::Type{Base.Pkg.Resolve.VersionWeights.VWPreBuildItem}) at pkg/resolve/versionweight.jl:80
  zero(::Type{Base.Pkg.Resolve.VersionWeights.VWPreBuild}) at pkg/resolve/versionweight.jl:120
  ...
Stacktrace:
 [1] _generic_matmatmul!(::Array{DNum{Int64},2}, ::Char, ::Char, ::Array{DNum{Int64},2}, ::Array{DNum{Int64},2}) at .\linalg\matmul.jl:506
 [2] generic_matmatmul!(::Array{DNum{Int64},2}, ::Char, ::Char, ::Array{DNum{Int64},2}, ::Array{DNum{Int64},2}) at .\linalg\matmul.jl:478
 [3] *(::Array{DNum{Int64},2}, ::Array{DNum{Int64},2}) at .\linalg\matmul.jl:141

[StefanKarpinski]

Comment out the promote_rule and convert definitions in examples/ModInts.jl and you get this:

julia> A = map(ModInt{13}, rand(0:12, 5, 5))
5×5 Array{ModInts.ModInt{13},2}:
  4  5   4  12  6
  3  1  12  11  6
 11  9  12   5  0
 12  6   2   3  7
  7  7   8  10  3

julia> A*A
ERROR: MethodError: Cannot `convert` an object of type Int64 to an object of type ModInts.ModInt{13}
This may have arisen from a call to the constructor ModInts.ModInt{13}(...),
since type constructors fall back to convert methods.
Stacktrace:
 [1] _generic_matmatmul!(::Array{ModInts.ModInt{13},2}, ::Char, ::Char, ::Array{ModInts.ModInt{13},2}, ::Array{ModInts.ModInt{13},2}) at ./linalg/matmul.jl:506
 [2] generic_matmatmul!(::Array{ModInts.ModInt{13},2}, ::Char, ::Char, ::Array{ModInts.ModInt{13},2}, ::Array{ModInts.ModInt{13},2}) at ./linalg/matmul.jl:478
 [3] *(::Array{ModInts.ModInt{13},2}, ::Array{ModInts.ModInt{13},2}) at ./linalg/matmul.jl:141

[andreasnoack]

I've looked a bit at the code and think that we should require that appropriate promotion and convert methods are defined for new arithmetic types if they are subtypes of something with some fallbacks defined like Number and AbstractMatrix. If you don't want your subtype to be subtype of anything, then I think we should require that you define a zero(::T) but not necessarily a zero(::Type{T}) method. Otherwise, we need to change the way we initialize the accumulators and I'm not sure how to do that in general. E.g. how would you compute the zeros of the result of a m x 0 times 0 x n multiplication.

[StefanKarpinski]

Right, fair enough.

[tkelman]

do you need to call zero if all dimensions are nonzero?

[StefanKarpinski]

Probably not, but it seems like a gotcha corner case – probably better to depend on zero up front so that people implement enough functionality that matmul works, full stop.

[andreasnoack]

do you need to call zero if all dimensions are nonzero?

You'd need to define conversion from Int to your type or zero. We should probably write a minimal custom-type-linear-algebra example into the documentation.