in-place assignment operator?

[StefanKarpinski]

Examples:

• sort!
• conj!
• transpose! — interesting question: is the operator version .'!?
• ctranspose! — likewise: is the operator version '!?

Please add others as you think of them. Obviously, they need to be operations that cannot change the type of the array. @ViralBShah: are there in-place FFTW operations that we should expose?

[ViralBShah]

FFT

various linalg functions from LAPACK

Also array arithmetic operations such as += etc.

[JeffBezanson]

• map! — Knowing the result type is super easy :)

[JeffBezanson]

We might want to change the way update operators like += are handled. Currently, x+=1 is replaced very early with x=x+1. Instead, we could lower it to x = (x+=1) where the definition of += defaults to +. Then we can make it a mutating operator for arrays.

[StefanKarpinski]

That definitely seems like a good thing to me. We want to be able to do things like X += 1 where X is a matrix in-place.

[JeffBezanson]

Actually this is a bit tricky. If we have a fallback definition of += so that it works for anything that defines +, then all mutable types are forced to implement all update operators, since otherwise they will work but not be mutating as expected. I think it would be better for the fallbacks to apply only to Number.
Do we have any current uses of these operators on whole arrays? I believe we don't.

[JeffBezanson]

For example, consider A\=b. In general you can't mutate A in such a way that it contains the result of A\b.

[JeffBezanson]

I'm worried that this makes code less generic. If f(x) contains x+=1, it suddenly becomes a mutating function if an array is passed in.
At the same time, it's not as flexible as it could be since you really want a three-argument version where the result of A+B is stored in C. How about some other syntax like

C := A+B

calls

+=(C, A, B)

?

[JeffBezanson]

Note the three argument version is needed for calling GEMM.

[StefanKarpinski]

I seems like a bit of a train wreck. If += is syntax then it works for immutables like numbers because x += y just means x = x + y, which reassigns x, giving it a new value. In order for += to work for arrays, it can't mean that because you want X += Y to mean something like +=(X,Y) where += is function that can mutate X. But then it will do absolutely nothing for immutables since they're immutable. The other issue is the f(x) containing x+=1 issue you mentioned. I'm not clear on how C := A+B solves that... is it because C is expected to be pre-allocated rather than created? I guess that makes sense, but only for arrays. I would be somewhat more inclined to call that .= since it does element-wise =.

[JeffBezanson]

Basically I don't want to have to think about whether to write x=x+1 or x+=1. We all habitually write x+=1 when possible, and I don't want to stop and think "oh, could x ever be an array?"

This is solved by using any syntax other than += for mutating arrays. C .= A+B is a ternary operator like ?:, that calls some function. C has to be pre-allocated because it is passed to this function.

[StefanKarpinski]

So this would mean that we can write things like A .= A+B and get efficient in-place element-wise addition of the elements of B to A. Although, this of course makes me want to be able to write something like A .+= B. Not entirely sure if I'm serious or not.

[StefanKarpinski]

Note the three argument version is needed for calling GEMM.

This is a really excellent reason to have this.

[StefanKarpinski]

How about having ternary forms like A .= B + C which means .+=(A,B,C), but also binary forms like A .+= B which means .+=(A,B) and can be defined to call .+=(A,A,B)? The .= form naturally generalizes to a varargs version: A .= B + C + D + E. Not entirely sure if this generalization makes sense for the binary form.

[StefanKarpinski]

Changing the title from the original "in-place versions of functions" since there are a lot of those now and people keep adding more. Now we just need to decide if we want some kind of in-place assignment operator or not.

[ViralBShah]

I should point out that in places where I have needed in place assignment, I have been able to use stuff like gemm! and passing the result as part of the input, or using loops and assigning into arrays. The code is not as elegant, but it's not too bad either.

[c42f]

Hi guys, this is an interesting discussion, and a thorny problem. Here's a naive suggestion for the in-place operators like +=. A lot of the problems seem to stem from needing different behaviour for mutable vs immutable. Given that += is special syntax already, can you just bite the bullet and put in a branch which checks for immutability? I'm imagining something like

if isimmutable(x)
    x = x + y
else
    addassign!(x, y)
end

with some sensible default implementation of addassign! which might depend on the existance of an in-place .=, thus:

function addassign!(x,y)
    x .= x + y
end

Of course, this depends on the compiler being able to optimize the branch away for good performance. I don't have enough experience with julia's dynamic type system to know whether this is plausible in the cases you care about performance.

[diegozea]

#3424 is related to this

[diegozea]

! in Julia's function names means mutation and . in functions names means element-wise. But given != is the different comparison operator, ! can not be used here. Is there a reason for use . instead of other symbols?

[StefanKarpinski]

This issue is not that it's impossible to achieve different behavior for mutable vs. immutable – the issue is that it's not desirable because that sabotages generically written code. For example, if you write x += y in some generic code, thinking that it means x = x + y, which is what it does mean for integers or floating point numbers, you expect this to have no effect on the caller's value of x. Then if you call the code on something mutable like a vector, you're inadvertently mutating the caller's value of x.

[c42f]

Hmm. On further thought what I proposed is rather similar to Jeff's x = (x+=1) idea, but worse. Sorry.

I remain hopeful that there's a solution which could make += work as a mutating update for arrays. It's very ingrained, and in practise I've never found it to be a problem in numpy which has semantics such that += is not mutating for numbers, but is for arrays.

[lindahua]

+1 for the element-wise assignment operator .= and the syntax a .+= b. This is consistent and clean.

[johnmyleswhite]

+1 for element-wise assignment via .=. If we can combine that with some of the devectorize macro's tricks, I think we'd be close to a general solution.

[StefanKarpinski]

@c42f, regarding NumPy, the thing is that one doesn't tend to write a lot of highly generic code in NumPy, and certainly not code that might be used on both mutable and immutable types of values. Then again, it's not 100% clear to me that this is actually a sufficiently realistic situation for Julia, but it still makes me uncomfortable enough that I don't want to do this without considering all the alternatives.

[simonbyrne]

While I appreciate the technical difference between += and a hypothetical .+= (or would +!= be more appropriate?), there is the potential that it might lead to code like:

if isimmutable(x)
    x += 1
else
    x .+= 1
end

in which we have the same result, but with much uglier code.

[lindahua]

I suspect how often people would ever write codes like

if immutable(x)
    x += 1
else
    x .+= 1
end

If you really want to write something as above and want the code beautiful, you can still write x += 1 without the if statement. If you really want inplace operation when x is an array. Then, you have to write

if isa(x, AbstractArray)
     for i in 1 : length(x)
         x[i] += 1
     end
else
     x += 1
end

I don't think this is any prettier than the code above.

The introduction of .= (and .+=) is to provide people additional options when they need. It does not break any existing code and suddenly make any existing way of writing codes infeasible. If one have a fantastic way of writing things without such operators, it is completely fine for him to continue to write codes in his way.

I have plenty of examples in writing numeric computation, machine learning, and image processing algorithms that need inplace updating of arrays --- to the extent that I feel I have to create a package and introduce functions like add! and multiply!. The introduction of .+= and .*= would have a huge benefit here.

I have a strong opinion about this, as I do such computation in a daily basis.

[StefanKarpinski]

@NHDaly, you're interpreting x .+= y as x .(+=) y whereas it means x (.+)= y i.e. x = x .+ y.

[NHDaly]

I see... To be honest that is terribly unclear. I think the language would
be better off if we simply removed these operators rather than have them
operate as they do now. It is not clear at all from context how they
behave. I think that is too subtle a distinction, and doesn't save much
typing.

I'm excited to see how this turns out! :)
On Jan 21, 2016 9:08 AM, "Stefan Karpinski" [email protected]
wrote:

@NHDaly https://github.com/NHDaly, you're interpreting x .+= y as x
.(+=) y whereas it means x (.+)= y i.e. x = x .+ y.

—
Reply to this email directly or view it on GitHub
#249 (comment).

[315234]

This immutable/mutable stuff is just way too subtle a distinction for most people. It is not intuitive that += and friends, especially operating on arrays, will create an entirely new object and destroy the old one. It is counterintuitive that doing an operation element-wise on an array (a[:] .*= a) has a different effect of using array expressions (a *= a). In a numerical language where performance is supposed to be a focus, it is almost never anyone's intention to implicitly allocate a new block of memory and throw out the old one. The end result is that Julia has array expressions but you basically can't use them without sacrificing performance.

[carnaval]

y = a
[...]
x = a
[...]
x += 1
[...]
print(y)

y = a
[...]
x = a
[...]
x = x+1
[...]
print(y)

I don't think that having those two pieces of code have different answers is at all intuitive

[StefanKarpinski]

This doesn't warrant any further discussion at this point and the mutable vs immutable distinction is completely orthogonal – the mutability of a type has no effect on the meaning of operations on it. We have a large number of in-place operators at this point, which is what this issue was originally opened to address.

[JaredCrean2]

Is this a definitive decision against including an in-place assignment operator in the language? It seems like the original question got lost in the mutable vs immutable discussion.

[johnmyleswhite]

I don't think anyone has a credible proposal for what such an operator would do. Everyone agrees the situation could be improved, but no seems to know how to do it.

[JaredCrean2]

Then shouldn't the issue be left open until someone figures it out? If the issue is closed, there is a chance no one will revisit it in the future.

[johnmyleswhite]

People can reasonably disagree about that point, but I personally think that's not an accurate assessment of how work gets done by most of the Julia developers. I prefer the open issues to be clearly actionable. But you don't need an issue to be free to work on something.

[JaredCrean2]

I guess that's reasonable. I definitely don't want this issue to be forgotten because it would make a big difference in much of the code I write.

[StefanKarpinski]

This particular issue has gotten too long and muddle at this point. It was also originally about adding more in-place operations, which we now have lots of. The issue of in-place operation will not be forgotten, I can assure you.

[stevengj]

Fixed by #17510, although we won't get the full power of fusing in-place assignment until .+ etcetera turn into fusing calls ala #16285.

[NHDaly]

Fascinating. Thanks for the update!

[bramtayl]

It might be useful for

f!(args...)

to be parsed to

inplace(x, f, args...)

and we could take advantage of new closure types to do:

inplace(x, ::typeof(f), args...) = f_in_place(x, args...)

where f_in_place is an internally defined in-place operator.

[KristofferC]

What's the gain? Just can't know for a general f how its f! should look so you need to define both anyway.

[bramtayl]

Yes, I do think you need to define both anyway. But having ! exposed to the parser would allow it to hook into any forthcoming more general loop fusion syntax (see #16285)

[GiggleLiu]

Now is 2018, having reached some conclusions yet? Somehow the gotcha 6 in this post still exists.

a = randn(4)
@allocated a .+= 2

will still produce 48!
I spent a lot of time handling this unexpect behavior. It is definitely a killer of julia beginers intended to write high performance packages.

Maybe 7 years is too short to make such a "huge" improvement?

[KristofferC]

julia> const a = randn(4);

julia> @allocated a .+= 2
0

Maybe reading the very first section of the performance tips is too much to ask?

[GiggleLiu]

@KristofferC
So sorry for posing such stupid benchmark. In fact, the real problem does not have such const concern

(I am trying to improve the performance of SparseMatrixCSC * Diagonal)

using SparseArrays
using LinearAlgebra
import LinearAlgebra:rmul!, mul!, lmul!
import Base: *

function rmul!(X::SparseMatrixCSC, A::Diagonal)
    nA = size(A, 1)
    mX, nX = size(X)
    nX == nA || throw(DimensionMismatch())
    @inbounds for j = 1:nA
        X.nzval[X.colptr[j]:X.colptr[j+1]-1] .*= A.diag[j]
    end
    X
end

using BenchmarkTools
const N = 1000
const dg = Diagonal(randn(ComplexF64, 1000))
const sp = SparseMatrixCSC(dg)
@benchmark rmul!(sp, dg)

It has 1000 allocations, although using code_warntype, I will see the type is stable. On the other side, writing the loop explicitly, I will get 10x speed up and 0 allocation!

Could you please explain it? (I promise I have already read the performance tips many times, lol)

[KristofferC]

It seems the view created from the broadcasted assignment is not elided and if the length of nzrange(X, j) is small, and this can be a significant slowdown.

A MWE is

julia> f(a) = (a[1:100] .*= 2; a)
f (generic function with 1 method)

julia> a = ones(1000);

julia> @btime f($a);
  112.160 ns (1 allocation: 896 bytes)

[StefanKarpinski]

That seems worth opening a separate specific issue about, no?

[mbauman]

This is an example of a place where call-site inlining is essential. It felt a little wrong to always inline copyto! — that's where the function break is here IIRC.