muladd(matrix, matrix, scalar) gives unexpected result

[stevengj]

This seems wrong:

julia> muladd(zeros(2,2), zeros(2,2), 3)
2×2 Matrix{Float64}:
 3.0  3.0
 3.0  3.0

It's not clear that muladd(X, Y, z) should even be defined when z is a scalar, since it is defined as X * Y + z and we stopped doing implicit vector + scalar broadcasting a long time ago.

This was added in https://github.com/JuliaLang/julia/pull/37065 by @mcabbot, but I don't see much discussion of the scalar case in that PR?

[mcabbott]

That PR changed the definition to X * Y .+ z, with the restriction that the answer not be larger than X * Y alone. The motivating case (as in the top post) wanted z::Vector. (That PR was a bit of a mess, sorry, as it was written under the belief that there no existing non-Number uses... was tightened in follow-up PRs.)

[LilithHafner]

For the evalpoly algorithm we'd want this to return [3 0; 0 3]. Subjectively, I think that [3 0; 0 3] makes more sense mathematically and [3 3; 3 3] is more useful for hand coded neural networks and other ML stuff.

Unfortunately, I don't think this is fixable in a non-breaking way. We could try nanosoldier, but I'm not optimistic.

[stevengj]

My preference would be for muladd(matrix, matrix, scalar) (or muladd(matrix, matrix, vector)) to simply throw an error, since we don't usually do implicit broadcasting in arithmetic operations.

But I agree that it's maybe too late to change in Julia 1.x.