Implement MutableArithmetics

[blegat]

MutableArithmetics is designed to be the API for mutable types in Julia. It allows both algorithms to be written in a generic way for any type while still exploiting mutability of any type implementing MutableArithmetics.
It's a bit similar to the mutability API implemented in AbstractAlgebra so I would like to try to merge them.
This PR is a first attempt at this. I implemented backward compatibility which shows what's the difference between the two API's.
In MutableArithmetics, we distinguish operate_to! and operate! which are what the user should use and that may or may not mutate and mutable_operate_to!/mutable_operate! which are what the types should implement and that should always mutate. It's also clearer that there should be not alias and we differentiate the cases with the _to suffix. It seems to me that the aliasing can be confusing. For instance,

GroupsCore.jl/src/group_elements.jl

Lines 263 to 267 in 10a138a

	function comm!(out::GEl, g::GEl, h::GEl) where {GEl <: GroupElement}
	# TODO: can we make comm! with 3 arguments without allocation??
	out = conj!(out, g, h)
	return div_left!(out, out, g)
	end

seems to be incorrect in case out === g since g is used in the last line.
However, this is the case here

GroupsCore.jl/src/group_elements.jl

Line 144 in 10a138a

res = comm!(similar(g), g, h)

As we need to use if-else everytime with === we might as well differentiate with the _to! suffix ?
Let me know what you think. I can do a similar PR to AbstractAlgebra if you think this goes in the right direction.

[kalmarek]

this comm! is actually correct as conj!, or mul will dealias elements if necessary.
In the test suite we check that none of the arguments (except the first) are mutated by the ! operation.

I'd be happy to accept this!
The plan is that this interface makes its way into AA (or AA will start depending on it). However if such change is rejected by AbstractAlgebra, I'm not sure how feasible would this be to have it here, but not there?

I'm alsow afraid that with such pull against AA you're facing uphill battle there... even implementing some kind of iteration on AA matrix elements took over a year :) Anyway, the discussion must start somewhere, so I'm summoning @fingolfin, @thofma, @rforquet, @fieker

[fieker]

On Wed, Mar 24, 2021 at 04:23:04PM -0700, Marek Kaluba wrote: this comm! is actually correct as `conj!`, or `mul` will dealias elements if necessary. In the test suite we check that both g and h are not mutated by the mutating operation. I'd be happy to accept this! The plan is that this interface makes its way into AA (or AA will start depending on it). However if such change is rejected by AbstractAlgebra, I'm not sure how feasible would this be to have it here, but not there? I'm alsow afraid that with such pull against AA you're facing uphill battle there... even implementing some kind of iteration on AA matrix elements took over a year :) Anyway, the discussion must start somewhere, so I'm summoning @fingolfin, @thofma, @rforquet, @fieker

What is the advantage of using this over the current (in AA and upstream) approach to mutability? In particular, given that we'd need to rewrite the entire generic to use this? A different point is, but that is easier to address: mutablity is not everything: delayed reduction is also important to all (most) quotient rings: e.g. given vectors over Z/nZ (n big) the "best" way to compute a scalar product is to compute the lifted scalar product followed by one reduction in the end, so this would have to be incorporated as well. I agree fully that, if this had been available 4 years ago, we'd possibly use it. But what is the advantage now of rewriting all foundations?

…

-- You are receiving this because you were mentioned. Reply to this email directly or view it on GitHub: #21 (comment)

[blegat]

this comm! is actually correct as conj!, or mul will dealias elements if necessary.
In the test suite we check that none of the arguments (except the first) are mutated by the ! operation.

IIUC, conj!(g, g, h) is allowed to modified g (but not forced to). In the fallback of conj! implemented in GroupsCore, it allocates a new element: inv(h) and modifies and returns this. However, some type of group element could have a faster implementation for this and be able to implement conj!(g, g, h) without any allocation by modifying g.
In that case,

GroupsCore.jl/src/group_elements.jl

Line 266 in 10a138a

return div_left!(out, out, g)

will not do what we expect because g has changed.

What is the advantage of using this over the current (in AA and upstream) approach to mutability?

I think there is two sides of this question:

1. What is the advantage of the MA mutable API over the AA mutable API
2. What is the advantage of having a package defining a mutable API outside of AA

I'd be happy to discuss 1. We could still change MA API if there are advantages in the AA mutable API but it's probably best to discuss 2. first.
Currently, if you define + and * for a mutable type then all of LinearAlgebra and SparseArrays will work but be very slow and allocate a lot. This is also the case for some functions in Base such as sum and prod.
In JuMP, it required a significant amount of work to dispatch matrix multiplication, sum, etc... to alternative implementations that would use the mutable API of JuMP.
But then for functions defined outside of Base and stdlib (say in some package X), if we wanted to make it work efficiently with JuMP, we should either make JuMP depend on X and rewrite it for JuMP types or make X depend on JuMP.
The first solution seems unmanageable and the second solution does not seem right: why would package X add such a large dependency just for having the function work better for JuMP types ?

So there was two downsides in having this mutable API specific to JuMP:

1. All the work of rewriting Base and stdlib was only used for JuMP types while it could be also useful for polynomials, etc...
2. Packages wanting to use the mutable API needed to have JuMP as a dependency which is a bit too much.

This was the motivation for creating MA and starting from JuMP v0.21, the transition was done and JuMP now use MA for its mutable API. This same line of reasoning seems to apply to AA with the exception that the mutable API of AA seems to be closer to an generic API than the one of JuMP before MA.

So if AA implements MA, then if the user calls MA.operate(*, A::Matrix, b::Vector) where A and/or b contains any ring element, this will use an efficient implementation instead of the default one, same thing if A is a sparse matrix.
In addition, if the ring element is a subtype of MA.AbstractMutable then A * b will be dispatched to MA.operate(*, A::Matrix, b::Vector) automatically.
If we have a significant amount of types implementing MA, it will provide sufficient motivation for packages to start relying on MA. For example, GenericLinearAlgebra or RowEchelon would directly use the MA API and no need for the user to use MA.operate or for MA to re-dispatch in case one of the arguments is a subtype of MA.AbstractMutable.
In the long term, we could even imagine having some part of MA be moved to Base so that LinearAlgebra and SparseArrays can directly be written on top of MA and we don't have to reimplement and re-dispatch in MA.

In particular, given that we'd need to rewrite the entire generic to use this?

Writing this PR gave me a good idea of the small differences between the AA and MA mutable APIs so I could make a PR to AA that implements MA while still keeping the AA API for compatibility for existing code using this API (like I did for this PR).

A different point is, but that is easier to address: mutablity is not everything: delayed reduction is also important to all (most) quotient rings: e.g. given vectors over Z/nZ (n big) the "best" way to compute a scalar product is to compute the lifted scalar product followed by one reduction in the end, so this would have to be incorporated as well.

That's interesting and I agree that this should also be part of the API. We have a similar need for MathOptInterface where you would want to canonicalize the functions only once at the end. How is that handled in AA ?

I agree fully that, if this had been available 4 years ago, we'd possibly use it. But what is the advantage now of rewriting all foundations?

I totally agree. I would say that it's better late than even later :) I have submitted a talk on MA to JuliaCon 2021 to try to advertise the idea that packages should use this API so that their code is also efficient for types in AA, GroupsCore, JuMP, MultivariatePolynomials, BigInt, ...

[thofma]

Base.LinearAlgebra is carefully constructed to work only in a numeric bubble (hello zero(typeof(x)) and det([1 0; 0 1])) , so adding more maintenance on our side for zero gain is not very appealing. There is also no interest to change this in Base and since code in the wild assumes that this is the interface, it will be broken by default for types in GroupsCore or AA too.

[fieker]

Apart from everything else, how do you handle recursinve mutable structures? E.g polynomials over BigRational, or better matrices of polynomials over BigRational: here the polynomial (as a matrix entry) might be mutable but the coefficients might not (shared between others)? By trait, all of matrix, polynomial and BigRational are mutable, but for any given object this is more tricky (and we don't have a (good) solution there as well. Consider creating a matrix - which will produce lots of zero entries. Each entry is a polynomial - maybe even the identical one. Take a zero polynomial and add x^n, then the sum needs to have many zeroes (for a dense poly). Will they be identical or different objects? MA.op!(+, poly_x, poly_y, poly_z): it is allowed to call MA.op! on the coefficients? ======================= We cannot use the Julia linear algebra, mostly, at all, zero(Type) can not be made to work for all our types. Putting types into numbers also makes no sense at all for finite fields, they are not real or complex. Similarly, matrices with zero rows/ columns do not have the information of their ring... So, out of neccessity, we build our own lin. alg. for our own matrix type(s). So: what is the benefit to us of also supporting MA? If there are intersting real applications, we can think about it, but at this point, there is no advantage for AA (or upstream), instead we have a new dependency and a new interface. It is an intersting interface, but too late for us ...

…

On Thu, Apr 01, 2021 at 02:24:53AM -0700, Benoît Legat wrote: > What is the advantage of using this over the current (in AA and upstream) approach to mutability? I think there is two sides of this question: 1. What is the advantage of the MA mutable API over the AA mutable API 2. What is the advantage of having a package defining a mutable API outside of AA I'd be happy to discuss 1. We could still change MA API if there are advantages in the AA mutable API but it's probably best to discuss 2. first. Currently, if you define `+` and `*` for a mutable type then all of LinearAlgebra and SparseArrays will work but be very slow and allocate a lot. This is also the case for some functions in Base such as `sum` and `prod`. In JuMP, it required a significant amount of work to dispatch matrix multiplication, `sum`, etc... to alternative implementations that would use the mutable API of JuMP. But then for functions defined outside of Base and stdlib (say in some package X), if we wanted to make it work efficiently with JuMP, we should either make JuMP depend on X and rewrite it for JuMP types or make X depend on JuMP. The first solution seems unmanageable and the second solution does not seem right: why would package X add such a large dependency just for having the function work better for JuMP types ? So there was two downsides in having this mutable API specific to JuMP: 1) All the work of rewriting Base and stdlib was only used for JuMP types while it could be also useful for polynomials, etc... 2) Packages wanting to use the mutable API needed to have JuMP as a dependency which is a bit too much. This was the motivation for creating MA and starting from JuMP v0.21, the transition was done and JuMP now use MA for its mutable API. This same line of reasoning seems to apply to AA with the exception that the mutable API of AA seems to be closer to an generic API than the one of JuMP before MA. So if AA implements MA, then if the user calls `MA.operate(*, A::Matrix, b::Vector)` where `A` and/or `b` contains any ring element, this will use an efficient implementation instead of the default one, same thing if `A` is a sparse matrix. In addition, if the ring element is a subtype of `MA.AbstractMutable` then `A * b` will be dispatched to `MA.operate(*, A::Matrix, b::Vector)` automatically. If we have a significant amount of types implementing MA, it will provide sufficient motivation for packages to start relying on MA. For example, [GenericLinearAlgebra](https://github.com/JuliaLinearAlgebra/GenericLinearAlgebra.jl) or [RowEchelon](https://github.com/blegat/RowEchelon.jl) would directly use the MA API and no need for the user to use `MA.operate` or for MA to re-dispatch in case one of the arguments is a subtype of `MA.AbstractMutable`. In the long term, we could even imagine having some part of MA be moved to Base so that LinearAlgebra and SparseArrays can directly be written on top of MA and we don't have to reimplement and re-dispatch in MA. > In particular, given that we'd need to rewrite the entire generic to use this? Writing this PR gave me a good idea of the small differences between the AA and MA mutable APIs so I could make a PR to AA that implements MA while still keeping the AA API for compatibility for existing code using this API (like I did for this PR). > A different point is, but that is easier to address: mutablity is not everything: delayed reduction is also important to all (most) quotient rings: e.g. given vectors over Z/nZ (n big) the "best" way to compute a scalar product is to compute the lifted scalar product followed by one reduction in the end, so this would have to be incorporated as well. That's interesting and I agree that this should also be part of the API. We have a similar need for MathOptInterface where you would want to canonicalize the functions only once at the end. How is that handled in AA ? > I agree fully that, if this had been available 4 years ago, we'd possibly use it. But what is the advantage now of rewriting all foundations? I totally agree. I would say that it's better late than even later :) I have submitted a talk on MA to JuliaCon 2021 to try to advertise the idea that packages should use this API so that their code is also efficient for types in AA, GroupsCore, JuMP, MultivariatePolynomials, BigInt, ... -- You are receiving this because you were mentioned. Reply to this email directly or view it on GitHub: #21 (comment)

[blegat]

Base.LinearAlgebra is carefully constructed to work only in a numeric bubble

If I understand you correctly, this is exactly the kind of issues listed in jump-dev/MutableArithmetics.jl#47. We have the same issues with MultivariatePolynomials and JuMP types and MutableArithmetics also solve some of these issues by re-dispatching the calls (for subtypes of AbstractMutable) to implementations that in addition to exploiting mutability, are more cautious with this, e.g. thanks to MA.promote_operation.

so adding more maintenance on our side for zero gain is not very appealing

This would be handled by MutableArithmetics so there will not be more maintenance for AA. Why would this be zero-gain ? The speed gain for matrix-vector multiplication or sum is not only a fixed speedup, it even changes the complexity in some cases.

There is also no interest to change this in Base and since code in the wild assumes that this is the interface, it will be broken by default for types in GroupsCore or AA too.

We are still far from it but as attested by jump-dev/MutableArithmetics.jl#47, we've already been able to merge a few PR's to make Base and stdlib more aware that more exotic types exist.

here the polynomial (as a matrix entry) might be mutable but the coefficients might not (shared between others)?

This is indeed tricky and we had exactly this use case with MultivariatePolynomials/DynamicPolynomials/TypedPolynomials.
MA.operate!(::typeof(op) redirects to op if it's not mutable and to mutable_operate! if it's mutable so if the polynomial uses MA.operate! on its coefficients, it will work both with mutable coefficients and non-mutable coefficients.
I have worked with 4 nested levels of mutability, e.g. array of polynomials with coefficients being MathOptInterface functions with BigInt coefficients, you can exploit the mutability at all levels with MA.operate!.
An example with array of polynomials of BigInt is here: JuliaMath/Polynomials.jl#331. You can see that you can reach an allocation-free operation, you can even detect that the inner level of BigInt multiplication needs a buffer and create it once before the operation on the arrays.

Take a zero polynomial and add x^n, then the sum needs to have many zeroes (for a dense poly). Will they be identical or different objects?

That's a tricky one. In MA, the convention is that: you to be sure that an object can be mutated, it should be the result of MA.copy_if_mutable or MA.operate (which again, works recursively). For instance, you may have x === +x but modifying MA.operate(+, x) through the MA API is guaranteed to have no effect on x.

MA.op!(+, poly_x, poly_y, poly_z): it is allowed to call MA.op! on the coefficients?

If poly_x modifies the the coefficients through any of the MA operations it should ensure that MA.copy_if_mutable and MA.operate creates independent copies for their returned object. So yes poly_x is allowed to call MA.operate! on its coefficients. But you may choose not to modify some fields through the MA API and then you don't need to copy them in MA.copy_if_mutable and MA.operate.

We cannot use the Julia linear algebra, mostly, at all, zero(Type) can not be made to work for all our types.
Putting types into numbers also makes no sense at all for finite fields, they are not real or complex.
Similarly, matrices with zero rows/ columns do not have the information of their ring...

I completely understand, I had the same issues with JuMP and MultivariatePolynomials as detailed in jump-dev/MutableArithmetics.jl#47.
We can go around some of them, sometimes by making PR's to Base (there, having more use-cases help of course).
Matrices with zero elements is indeed difficult but it might still be useful to have it work for matrices with a nonzero number of zero elements.

[blegat]

I understand that there is a large ecosystem around AA with a lot of code using it's mutable API and a lot of types implementing it. Hence making a breaking change would be a lot of work and you need a very good incentive to do it.
So to clarify, I am not suggesting to do anything breaking, we could keep the two APIs. For instance, we would add

function AbstractAlgebra.inv!(out, g)
    if out === g
        return MA.operate!(inv, g)
    else
        return MA.operate_to!(out, inv, g)
    end
end
function AbstractAlgebra.mul!(out, g, h)
    if out === g
        return MA.operate!(*, g, h)
    elseif out === h
        return MA.operate!(mul_left, h, g)
    else
        return MA.operate_to!(out, *, g, h)
    end
end

Then, for types that implement the AA mutable API, nothing changes and for types that implement the MA API, they are now usable with either the AA API and the MA API.
So if we implement the MA API for some types in AA, the users that are using the AA API won't see their code break.

[kalmarek]

@blegat I've just sent the version 0.5 for registration; since mutable API is specifically described as unstable I think it's fine to rebase this and release under 0.5.*