Make matrix multiplication work for more types #18218

blegat · 2016-08-24T14:00:35Z

Matrix multiplication does not work for arrays of elements types T and S such that the type TS = promote_op(*, T, S) does not satisfy TS == promote_op(+, TS, TS).

For example, in this package, I have the type Term which is one term of a polynomial and the type Polynomial which is the sum of several terms.
Suppose I am trying to compute A * x where A is a matrix of Int64 and x is a vector of Term.
Then TS will be Term so it will allocate a vector of Term to store the result. However when doing the sum of Term it will get Polynomial so what should have been allocated is a vector of Polynomial.

The fix is simple, allocate an array of type promote_op(+, TS, TS).

joehuchette · 2016-08-24T14:04:46Z

cc @mlubin, this may help with JuMP's internal matmul code.

andreasnoack · 2016-08-24T21:28:27Z

Couldn't this just be promote_op((t,s) -> t*s + t*s, T,S) instaed of a new function?

mlubin · 2016-08-24T21:30:08Z

@joehuchette, yep

pabloferz · 2016-08-25T01:19:53Z

Along the lines of @andreasnoack, you could define something like (probably with a better name) matprod(x, y) = x*y + x*y and just change * by matprod in the promote_ops

blegat · 2016-08-25T09:27:35Z

In fact, since 3ace05a and 6fd91b2 were merged, this change has become a lot more important. Indeed, before these commits, promote_op would just return Any for non-Number types and for a Number type TS, we have promote_op(+, TS, TS) == TS.
Therefore, since those commits have been merged into v0.5, I think it would be appropriate to merge this pull-request also into v0.5.

@andreasnoack @pabloferz yes we could also do something like that. The main difference is that in the solution of this pull request, if promote_op(*, .., ..) and/or promote_op(+,...,..) gets a specialized methods then it will be taken into account.
What is the advantage of having just one promote_op as you suggest ?
If this is the fact that I created a new promote variant, we can rename it. I just created a new function to avoid duplicated code and added promote_ as prefix because it does some kind of promotion.

martinholters · 2016-08-25T09:38:47Z

In the long run, I hope (specializations of) promote_op will disappear. For 0.5, if we backport this PR, it might be a good idea to honor specializations of promote_op for * and + (although making sure the called functions are inferable is better than specializing promote_op already).

pabloferz · 2016-08-25T09:38:54Z

I don't think that the plan is to eventually keep overloading promote_op. The work that was done, was to achieve the opposite (remove specializations). If for some reason promote_op gives you Any, instead of specializing it, you should check for type instabilities on your methods definitions.

What is the advantage of having just one promote_op as you suggest ?

Since promote_op relies on inference, and inference does not have a well defined behaviour, I don't think its safe to call promote_op on top of promote_op. In fact, the errors you get with the current proposal, I believe, are due to this.

I tried locally the alternative proposed solution and it works just fine.

blegat · 2016-08-25T10:16:47Z

@pabloferz I might not understand inference correctly but what scares me a bit is that matmul is not a generic addition of 2 elements of type typeof(t*s). It is t*s + t*s so the 2 elements are the same.
In my example of polynomial, since t*s and t*s are 2 terms of the same monomial, their sum is still a term (and not a polynomial !). In fact the sum of 2 terms can either be a Term of a Polynomial. If the inference was very smart it would return Union{Term, Polynomial}. Currently, both promote_matmul and promote_op(matmul, ..., ...) returns Any which is also very good :)

If you guarantee me that the inference will not use the fact that t*s and t*s are the same and only use their type then your solution is perfect :)

pabloferz · 2016-08-25T13:57:14Z

base/linalg/matmul.jl

@@ -424,7 +426,7 @@ end
 function generic_matmatmul{T,S}(tA, tB, A::AbstractVecOrMat{T}, B::AbstractMatrix{S})
    mA, nA = lapack_size(tA, A)
    mB, nB = lapack_size(tB, B)
-    C = similar(B, promote_op(*, arithtype(T), arithtype(S)), mA, nB)
+    C = similar(B, promote_op(matprod, arithtype(T), arithtype(S)), mA, nB)


I believe you missed two of these below (for matmul2x2 and matmul3x3)

Yes, you are right. Why aren't they using arithtype by the way ?

I guess they just were unintentionally forgot. But, now that you mention it, matprod also solves (more generally) the problem they where supposed to tackle, so in fact we can remove them and their definitions.

Probably by mistake. But IIUC, arithtype is only needed to convert Bool to Int when doing matrix multiplication, which the promote_op(matprod, ...) does even without the arithtype, so no worries here, but maybe a chance for cleanup.

blegat · 2016-08-25T16:30:17Z

I have removed arithtype with 2286b3e :)

pabloferz · 2016-08-25T16:35:31Z

test/linalg/matmul.jl

+  import Base.*, Base.+, Base.zero
+  immutable TypeA
+    x::Int
+  end


Not that it matters too much, but you can save some space by changing, e.g. TypeA by Int, so you only need two new types for the test.

blegat · 2016-08-25T22:52:50Z

@tkelman Fixed :)

martinholters · 2016-08-26T06:35:47Z

I have removed arithtype with 2286b3e :)

Looks like it's being used in DistributedArrays.jl, so better to deprecate than to remove.

blegat · 2016-08-26T08:30:24Z

@martinholters I have added a deprecation warning

KristofferC · 2016-08-26T08:52:00Z

base/deprecated.jl

@@ -778,6 +778,18 @@ end
 @deprecate cholfact!(A::Base.LinAlg.HermOrSym, uplo::Symbol, ::Type{Val{false}}) cholfact!(A, Val{false})
 @deprecate cholfact!(A::Base.LinAlg.HermOrSym, uplo::Symbol = :U) cholfact!(A)

+# #18218
+function arithtype(T)


Can't you use the @deprecate macro like used above?

Doesn't add @deprecate automatically add an export? We probably don't want that here.

I don't think I can because the macro is when a function should be used instead of another one. Here there is not replacement.

Makes sense. Thanks.

tkelman · 2016-09-23T01:20:30Z

base/deprecated.jl

+# #18218
+eval(Base.LinAlg, quote
+    function arithtype(T)
+        depwarn("arithtype is deprecated as it is no longer needed in Julia Base.",


what's the replacement?

Instead of using arithtype which replaces Bool by Int for matrix multiplications, we use an internal method matprod(x, y) = x*y + x*y which should cover the roll of arithtype while also being more general.

are packages using this? if so, the deprecation should tell them what to use instead

With this, I mean that there is no need now for arithtype or anything like it, which used to operate on types. Now its functionality (but not the method itself) is replaced by the definition and use of matprod.

If packages are using this together with promote_op then the new way of doing things would be promote_op(matprod, S, T) instead of promote_op(*, arithtype(S), arithtype(T)). If packages are using it some other way, they would have to define their own replacement, I guess.

this is not a helpful warning

Yeah, you're right. @blegat Can you specify something more useful? Maybe like

"arithtype is now deprecated. If you were using inside a promote_op call, use promote_op(LinAlg.matprod, Ts...) instead. Otherwise, if you need its functionality, consider defining it locally."

Of course, it doesn't have to be that if you find a better way to express it.

LinAlg.matprod right?

Yeah, that's the one (edited above).

Thanks for your help :) I have updated the warning message with your suggestion.

tkelman · 2016-09-23T10:06:09Z

base/deprecated.jl

+# #18218
+eval(Base.LinAlg, quote
+    function arithtype(T)
+        depwarn("arithtype is now deprecated. If you were using it inside a promote_op call, use promote_op(LinAlg.matprod, Ts...) instead. Otherwise, if you need its functionality, consider defining it locally.",


wrap this over multiple lines

Sorry about that, this is now fixed :)

andreasnoack · 2016-10-24T00:56:58Z

@blegat There seems to be a small conflict here. Can you please rebase. Then it should be ready to merge.

Currently it is assumed that the type of a sum of x::T and y::T is T but this may not be the case

blegat · 2016-10-24T01:49:22Z

@andreasnoack Thanks, that is now rebased :)

Same idea as JuliaLang/julia#18218

tkelman · 2016-10-27T01:04:07Z

cc @ajkeller34 it appears this broke Unitful.jl http://pkg.julialang.org/logs/Unitful_0.6.log

ajkeller34 · 2016-10-27T02:16:08Z

I just read that log quickly but it seems like all of those failures could be fixed by defining:

+{T<:Units}(::T, ::T) = T()

I guess that's not unreasonable: m+m may as well equal m. I'll try it out on master and see if any further errors happen upon implementing that. If things start getting hairy with promotion, I'll probably wait until 0.5.1 comes out with the fix for #18465 until working on this carefully.

ajkeller34 · 2016-10-27T06:08:49Z

@tkelman Fixed.

tkelman · 2016-10-27T06:26:56Z

I'm not sure having to define addition on types is really the best solution here?

ajkeller34 · 2016-10-27T07:43:30Z

Yes, ordinarily addition is only defined for quantities, not for units (not to be pedantic, but I've defined it for Units objects, not types). I haven't added docs or tests for this behavior because it is a temporary solution, which shouldn't be required in the future. There are people running Unitful who want it to work properly on 0.6, and this patch immediately restores functionality to valid operations with slim chance of causing problems for anyone.

I had planned on putting more effort/thought into promotion in Unitful once 0.5.1 is out. When inference is working well, it's my understanding that I shouldn't need any of the several promote_op definitions I have---as you know, promote_op is not even documented in the release-0.5 docs. I intend to clean up everything in one go and then require julia 0.5.1. I'll report back here or in a separate issue if I can't get this to work some other way.

* Make matrix multiplication work for more types Currently it is assumed that the type of a sum of x::T and y::T is T but this may not be the case * Remove arithtype in matmul and deprecate it

(cherry picked from commit 6ba63e5)

kshyatt added linear algebra Linear algebra types and dispatch Types, subtyping and method dispatch needs tests Unit tests are required for this change labels Aug 24, 2016

blegat force-pushed the matmul branch 3 times, most recently from 05e2791 to c5e38fa Compare August 25, 2016 13:43

pabloferz reviewed Aug 25, 2016
View reviewed changes

blegat force-pushed the matmul branch from c5e38fa to bacf5ba Compare August 25, 2016 14:31

pabloferz reviewed Aug 25, 2016
View reviewed changes

blegat force-pushed the matmul branch 2 times, most recently from 679d266 to 0157d0c Compare August 25, 2016 22:51

kshyatt removed the needs tests Unit tests are required for this change label Aug 25, 2016

blegat force-pushed the matmul branch from 0157d0c to f9dc357 Compare August 26, 2016 08:22

KristofferC reviewed Aug 26, 2016
View reviewed changes

blegat force-pushed the matmul branch 2 times, most recently from ff59690 to 81fbf1f Compare August 26, 2016 09:15

tkelman reviewed Sep 23, 2016

View reviewed changes

blegat force-pushed the matmul branch from c694678 to 39b794d Compare September 23, 2016 10:01

tkelman reviewed Sep 23, 2016

View reviewed changes

blegat force-pushed the matmul branch from 39b794d to dd36210 Compare September 23, 2016 10:42

andreasnoack approved these changes Oct 24, 2016

View reviewed changes

blegat added 2 commits October 23, 2016 21:44

Make matrix multiplication work for more types

4df067b

Currently it is assumed that the type of a sum of x::T and y::T is T but this may not be the case

Remove arithtype in matmul and deprecate it

bd7e894

blegat force-pushed the matmul branch from dd36210 to bd7e894 Compare October 24, 2016 01:46

andreasnoack merged commit 1efe487 into JuliaLang:master Oct 24, 2016

andyferris pushed a commit to JuliaArrays/StaticArrays.jl that referenced this pull request Oct 25, 2016

Improved output type inference for matrix products

a310bf5

Same idea as JuliaLang/julia#18218

blegat mentioned this pull request Mar 8, 2017

DimensionMismatch when performing vector-matrix multiplication with custom types jump-dev/JuMP.jl#988

Closed

Sacha0 added a commit to Sacha0/julia that referenced this pull request May 18, 2017

Add NEWS.md entry for arithtype deprecation (JuliaLang#18218).

5580725

Sacha0 added needs news A NEWS entry is required for this change deprecation This change introduces or involves a deprecation labels May 18, 2017

ararslan pushed a commit that referenced this pull request May 18, 2017

Add NEWS.md entry for arithtype deprecation (#18218). (#21942)

6ba63e5

Sacha0 removed the needs news A NEWS entry is required for this change label May 19, 2017

tkelman pushed a commit that referenced this pull request Jun 3, 2017

Add NEWS.md entry for arithtype deprecation (#18218). (#21942)

758c817

(cherry picked from commit 6ba63e5)

blegat mentioned this pull request Mar 6, 2018

Make LU factorization work for more types #26344

Merged

mlubin mentioned this pull request Sep 6, 2018

clean up some linear algebra primitives for JuMP eltypes jump-dev/JuMP.jl#1456

Merged

blegat mentioned this pull request Nov 7, 2018

Implement value for GenericQuadExpr. jump-dev/JuMP.jl#1588

Merged

blegat mentioned this pull request May 5, 2020

Upstream jump-dev/MutableArithmetics.jl#47

Open

8 tasks

blegat mentioned this pull request Jun 16, 2020

Fix Matrix on diagonal and triangular matrices for types without zero #36256

Open

Make matrix multiplication work for more types #18218

Make matrix multiplication work for more types #18218

Conversation

blegat commented Aug 24, 2016

joehuchette commented Aug 24, 2016

andreasnoack commented Aug 24, 2016

mlubin commented Aug 24, 2016

pabloferz commented Aug 25, 2016 • edited Loading

blegat commented Aug 25, 2016 • edited Loading

martinholters commented Aug 25, 2016 • edited Loading

pabloferz commented Aug 25, 2016 • edited Loading

blegat commented Aug 25, 2016

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

blegat commented Aug 25, 2016

pabloferz Aug 25, 2016 • edited Loading

Choose a reason for hiding this comment

blegat commented Aug 25, 2016

martinholters commented Aug 26, 2016

blegat commented Aug 26, 2016

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

pabloferz Sep 23, 2016 • edited Loading

Choose a reason for hiding this comment

tkelman Sep 23, 2016 • edited Loading

Choose a reason for hiding this comment

pabloferz Sep 23, 2016 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

andreasnoack commented Oct 24, 2016

blegat commented Oct 24, 2016

tkelman commented Oct 27, 2016

ajkeller34 commented Oct 27, 2016

ajkeller34 commented Oct 27, 2016

tkelman commented Oct 27, 2016

ajkeller34 commented Oct 27, 2016

pabloferz commented Aug 25, 2016 •

edited

Loading

blegat commented Aug 25, 2016 •

edited

Loading

martinholters commented Aug 25, 2016 •

edited

Loading

pabloferz commented Aug 25, 2016 •

edited

Loading

pabloferz Aug 25, 2016 •

edited

Loading

pabloferz Sep 23, 2016 •

edited

Loading

tkelman Sep 23, 2016 •

edited

Loading

pabloferz Sep 23, 2016 •

edited

Loading