Julia: 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]

https://github.com/JuliaLang/julia/issues/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.

[simonbyrne]

@lindahua My statement wasn't intended to be an argument against the introduction of in-place operators (your packages have clearly demonstrated how useful they are), but merely pointing out that introducing new in-place operators that are distinct from the "new-variable-and-change binding" operators that currently exist may have some unintended effects in how julia is used.

[cdsousa]

+1 for the .= + and .= * ternary operator with varargs versions.

[cdsousa]

Ok, let me add fuel to the fire :)

Perhaps the output parameters can be passed wrapped in some new type, let's say Output:

immutable Output{T} p::T end

Then, for a common function, e.g., the sum of two arrays,

function +{T}(a::Array{T, 1}, b::Array{T, 1})
...
end

one could define an in-place assignment sum function as

function +{T}(c::Output{Array{T, 1}}, a::Array{T, 1}, b::Array{T, 1})
    c = c.p  # BIG NOTE: it would be much nicer if c could be used directly...
    l = min(length(c), length(a), length(b))
    for i=1:l c[i] = a[i] + b[i] end
    return c
end

that would be called in the in-place assignment expression (with := or whatever the chosen operator )

c := a + b

which the interpreter would translate to

+(Output(c), a, b)

That can be simulated with the following @ipa (in-place assign) macro:

macro ipa(assignment)
    @assert assignment.head === :(:=)
    call = assignment.args[2]
    @assert call.head === :call
    e = assignment.args[1]
    output = :(Output($(assignment.args[1])))
    call.args = [call.args[1], output, call.args[2:]]
    eval(call)
end

for which the example

a = [1, 2, 3]
b = [4, 5, 6]
t = Array(Int, 3)
dump(t)
c = t
c = a + b # c is now referencing a new allocated array
dump(t)
c = t
@ipa c := a + b # array referenced by c is modified
dump(t)

prints

Array(Int64,(3,)) [14339488,14338960,16551936]     <-- garbage
Array(Int64,(3,)) [14339488,14338960,16551936]     <-- still garbage
Array(Int64,(3,)) [5,7,9]

The big advantage of this approach is that it can be used for any function (even accepting varargs), needing only the write of an in-place function version.

[nalimilan]

I'm a big fan of .=, .+= and similar functions too. In many cases element-wise in-place operations will avoid many memory issues when dealing with large datasets.

Though I wonder whether += couldn't be made to behave the same way. The problem raised above is that when x is mutable, x += y would not behave like x = x + y. But wouldn't the compiler be able to make them behave the same in all cases? I mean, when you do something like this function test(x::Array) = (x = x + 1), what happens is that x in the local scope points to a copy of the caller's x. Why couldn't the compiler copy x when it encounters x += y and that x isn't in the local scope?

If you did that, then .+= or (more logically) +=! could be defined to be _mutating_ operators, if that's deemed useful.

[JeffBezanson]

You have to make a distinction between mutation as an optimization, and semantic mutation (code that explicitly says something must be mutated, no matter what). Of course it is good to do x = x+y in place automatically if there is only 1 reference to x. I only object to += unconditionally mutating its left-hand side (for some argument types!)

There could be a version of + that has an output keyword argument: +(a::Array, b::Array; output=similar(a)). If the compiler sees x = x+y and it can prove x safe to overwrite, it could rewrite this to +(x, y, output=x). So far this is what I can come up with for generalizing "3-address code" (a compiler's bread and butter) to arbitrary types.

[lindahua]

As one can do x[:] += y, this doesn't seem to me a pressing issue.

[StefanKarpinski]

Another argument against having += do in-place array mutation:

v = zeros(Int,5)
v += 0.5 # InexactError

[andreasnoack]

As just mentioned in #8699, we might pay a big price for not allowing different behavior of immutable and non-immutables here.

However, we should think a little about commutativity. For matrices, in place linear algebra works for triangular multiplication and solves, but also when applying rotation, reflector and Fourier transform matrices as well as factorized matrices. General matrices should throw an error.

I think the rules would have to be something like

• A*=B is A=B*A for triangular B.
• A\=B is A=B\A for triangular B.
• A/=B is A=A/B for triangular B.

and these rules would also work for BigFloats.

Next problem would be to think about if we can and want to handle the gemm operation C=α*A*B+βC in place. It would be cool if that also generalized to BigFloats with α and β being Float64s, but I don't think it is possible to do that without temporary allocation.

[johnmyleswhite]

Personally, I'm not a fan of the way that *= conflates syntactic sugar with in-place operations. I think we need to focus on a generic approach to in-place operations without using *= operations as a crutch. In particular, I worry that the *= approach will solve enough problems that we'll do a bad job solving the other cases when in-place operations need special syntax.

[toivoh]

*= switching the order of the factors seems very dangerous and confusing to me. I also agree with @johnmyleswhite.

[StefanKarpinski]

@johnmyleswhite has hit the nail on the head here – letting *= do mutation as a performance hack is just sweeping a more general problem under the rug.

[simonbyrne]

I agree with @johnmyleswhite, though I suspect the general solution is far from trivial.

Think about a basic Newton's method algorithm

function opt(x)
    while true
        H = hess(x)
        g = grad(x)
        delta = H\g
        x -= delta
        norm(delta) < 1e-8 && return x
    end
end

Ideally, we would want hess, grad and \, and - to reuse H, g, delta and x, and any internal work arrays that they use, across iterations. It might also be possible to use the same array for g and delta (depending on the algorithm used for \), but I don't think this will result in a big change here.

[JeffBezanson]

A factor of 4 in performance, to me, is not a "big price" compared to the _disaster_ that this kind of premature optimization would cause. The long term goal should be to remove ! functions. Clever _new_ syntax for mutation might be fair game, but we need real optimizations. Admittedly it is very hard.

[andreasnoack]

focus on a generic approach to in-place operations

I order to get anywhere, I think we'll have to define some kind of class or characteristics of what could potentially be done in place with a kind of assignment syntax. Talking of the _general_ problem will lead nowhere. Multiplication and division of triangular matrices is one class where it might be possible to define reasonable update behavior with =.

conflates syntactic sugar with in-place operations

Suntactic sugar is one interpretation of a*=b. It doesn't have to be the right one. Couldn't the "rule" for the syntax be: update a with the value of a*b.

@JeffBezanson _Premature optimization_? This issue is more than three years old and I still have to write Fortran style code with nested loops and BLAS calls when doing _technical computing_ in Julia in order to get a speed anywhere near Fortran.

[simonbyrne]

To follow up to my earlier comment, one thing I wanted to do with InplaceOps.jl was to be able to indicate reuse of arrays in loops (see this post).

Although that may help with some issues, it would still be very fiddly. One approach I was thinking of is that in the _method definitions_ you would indicate which arrays are reusable, e.g.

function \{T}(A::Matrix{T}, b::Vector{T})
    @reusable y = Array(T,size(A,1))
    @reusable F = factorize(A)
    gemv!(A, b, y) #I know this isn't correct, but you get the idea
end 

Then when \ is called in a loop, if julia can determine that F and y are not used outside the scope of the loop, the arrays could be reused (obviously this would require checking to see if the array sizes get changed between calls).

Now, I'm sure there's a lot more required to get this to work, and I'm not a computer scientist, so I have no idea if something like this is at all feasible.

[jiahao]

In my experiments with IterativeSolvers.jl, the cost of not writing explicit in-place operations is closer to 50x-500x, not 4x.

[JeffBezanson]

This discussion got mixed up with the one about BigInts. I'm not against mutating and reusing arrays. The only changes I oppose are (1) making numbers mutable, (2) making += mutating for some types and not others. In fact one possible design is to make += _always_ a mutating operator, in which case it can't be implemented for numbers. On the whole I wouldn't prefer that, but it would be "safe" by my standards.

[ViralBShah]

Matlab does a pretty good job of reusing temporary arrays in loop iterations. I wonder if they have done this through escape analysis or if it is a well tuned GC or something else.

[ViralBShah]

I think that with various packages addressing this, and the various ! methods in Base, we can close this issue.

[kmsquire]

A little disappointed. I kinda like Jeff's := proposal (although I vaguely recall problems with it discussed in other issues?)

[JeffBezanson]

There are probably issues more worthy of closing than this one. ! methods are not very satisfactory, and I think it's likely we should try to do something lower level about in-place-ness.

[StefanKarpinski]

Agree on keeping this open. Still prefer .= to :=.

[ViralBShah]

I guess this issue started with the ! methods but then moved on to operators. Would be nice to have a separate issue - but perhaps it is simpler to just leave this open.

[ViralBShah]

With . operators meaning element-wise / broadcast operations, it seems like overloading the syntax too much to also include in-place operations.

[StefanKarpinski]

That's the point – it's an elementwise assignment.

[kmsquire]

I'm okay with .= as well. :-)

[Glen-O]

I may be the only one, but I feel like "A+=B" should inherently be in-place. That's why you'd use it over A=A+B - at least, until I'd recently learned about a speed issue with .*=, I didn't realise that it didn't default to in-place where possible.

I'm not sure I like .= in concept, because it feels like you should then be able to type A+.=B to perform matrix addition in-place, with elementwise addition being A.+.=B... but then it becomes confusing, with +=, +.=, .+=, and .+.=. I can see the argument for the dot meaning "elementwise", but it seems like it's a little weird as a notation.

As an alternative, how about a special notation, something like =<, to say "assign in place", with the operation failing and throwing an error if it the allocated space doesn't match or the variable is of an immutable type? Then all of the usual operators feel "normal" attached to it - A+=<B means "take A, and add B to it in place", for instance.

Not only does this make mutation more explicit, but provides a neat way to tell if there's a mistake in the code. And then .= could be used for actual elementwise assignment - that is, A.=1 would construct a matrix similar to A filled with 1 (basically, A=fill(1,size(A))), and then A.=<1 would be fill!(A,1). And if A were 2x2, and you wrote A.=[1 2], then you would get A==[1 2;1 2], because elementwise operations broadcast.

I'd also like to support what @jiahao said about the performance hit with non-in-place operations. I found the GeneticAlgorithms package insufficient for my need, so I wrote my own GA code... and when I went through and pre-allocated all my arrays and then mutated them instead of just assigning new arrays to the same variable, it had a massive impact. Easily more than 10x speedup. But it became a bit messy because I had to explicitly go through and do a heap of switching to forms like A[:]=B[:] (not actually from the code). The ability to write this as A=<B would make this so much easier.

Incidentally, why _are_ numbers immutable? Is there some performance or design issue with re-using the memory that a number variable is stored in? I'm not suggesting they should be mutable, just wanting to understand the reasoning.

[StefanKarpinski]

Think about what it would mean if you could change the value of a number.

[lucasb-eyer]

I also expected += to be in-place and got bitten by it not being so. It's not just about performance, it introduced an actual bug in my code.

But if += did different things for scalars than for arrays, it'd become useless in generic code.

I think at the very least, this should be made much more explicit in the documentation, especially given that += _is_ in-place in numpy. The note on this page was the only reference I found in the docs, and it's very implicit about the implications.

[KristofferC]

Doesn't this part make it clear?

An updating operator rebinds the variable on the left-hand side.

[Glen-O]

@StefanKarpinski - If you're answering the "why are they immutable" question, I should make more clear that I mean variables that are number types, not actual numbers. Is there a benefit to always rebinding numerical variables instead of updating them in-place?

Also, what is your opinion of an explicit "in place assignment" operator (even if it always throws an error for attempts to in-place assign numbers)?

[lucasb-eyer]

@KristofferC Technically, yes. Practically, it wasn't enough to let "A = ones(4); B = A; B += 3 doesn't change values in A" sink in - for me, anyways. This is a big difference to Numpy, but also typical C++ matrix classes. Maybe adding this as a point to the Noteworthy differences to Python might help mitigate this a bit.

[jiahao]

@lucasb-eyer the difference from Python is worth documenting, if you'd like to submit a PR

[johnmyleswhite]

Is there a benefit to always rebinding numerical variables instead of updating them in-place?

@Glen-O: The fact that number r-values are immutable _is the primary reason_ why you can optimize code that involves scalar computations. If that's unclear to you, you might benefit from learning more about how Julia's compiler generates machine code. Make sure you understand (1) what l-values and r-values are, (2) what stack and heap allocation are, and (3) how immutability and stack allocation are related.

[StefanKarpinski]

Bindings to numbers aren't immutable or rather const unless you declare them to be, so I'm kind of unclear on what you mean, @Glen-O.

[lucasb-eyer]

@jiahao Glad you agree, I'll prepare a PR in a few days, got a talk to prepare first :fearful:

[Glen-O]

@StefanKarpinski - I'm not referring to the bindings, but what you're binding to. Please note that I'm a mathematician, not a computer scientist. What I'm trying to understand is, in what way is it better to keep re-binding a numeric variable to a new location in order to give it a new value, rather than changing the value in-place? If you have a=1, and you want to increment a to equal 2, why is it better to create a new value of 2 in a new section of memory and then point a at it, rather than simply changing the value in the bit of memory that currently holds the value 1 so that it now holds the value 2? Note that I'm not specifically referring to function calls, but generally. If I run a=1;@time for i=1:1000000 a=a+1;end (on 0.3), it tells me that 15991840 bytes allocated. And if I run it a few times over, it has to run gc.

I do realise that it will do much better inside a function, but I just can't see why it should ever need to allocate so much for what amounts to incrementing a value. I was hoping for a quick, lay explanation since there is clearly some reason why it was chosen to do it like that.

Mind you, it's just a sidebar, anyway - I'm far more interested in whether or not my suggestion of a dedicated "in-place assignment" operator, separate from elementwise considerations and only usable on mutable objects with constant footprint (assigning a 5x4 integer array to a variable that currently holds a 5x4 integer array, for instance - possibly also working if it's a 20x1 integer array?), is of interest.

[nalimilan]

@Glen-O This is completely unrelated to the issue at stake. The allocation happens because in the global scope a might change type during the loop (e.g. it might become an array). Always run benchmarks inside functions, as the FAQ explains.

[Glen-O]

@nalimilan - let me try asking a different way, since I think you got caught on the specific example, which I used to demonstrate that it's rebinding even though it clearly wasn't going to change type. So here's how I'll ask it - is there a reason why the "in place assignment operator" as I suggested above shouldn't be able to assign to an Int variable in-place if you're trying to give it a new Int value (that is, why shouldn't it be mutable when it remains the same type)?

Alternatively, if it's because it can't be certain that the type will remain the same, then why doesn't this fix it? a=0;@time for i=1:1000000 a::Int=(a+1)::Int;end. Now it has specifically been told that a and a+1 are both Int type, so it should be able to safely mutate, right?

Incidentally, if this discussion should end here (I can understand if you wish to focus on the specifics of the in-place assignment operator considerations, rather than this particular detail of in-place vs immutable), I'm happy to move it to julia-users - just say the word.

[hayd]

Now it has specifically been told that a and a+1 are both Int type, so it should be able to safely mutate, right?

Which is exactly what the compiler does in the function scope (but Ints are never going to be mutable in julia itself.) This would make a good Stack Overflow/julia-users question, but as mentioned above it's off topic for this issue.

[ihnorton]

if this discussion should end here

Yes. It is _really great_ to have people trying to understand the why of these design decisions, but GitHub can't effectively facilitate those discussions or archive them usefully and searchably. Questions like this are perfect for the StackOverflow model because the answers are timeless (and relevant beyond Julia). Thanks!

[NHDaly]

The biggest thing to me is that Arrays are supposed to be mutable containers of possibly immutable values. So it feels like an element-wise operation should _mutate_ the container, since you are only touching the elements.

A clarifying example:

It seems wrong that these two do different things:

julia> x = [ 1 1 1 ]
julia> y = x        # setup x === y

julia> for i in 1:length(x)
         x[i] += 1  # x still === y
       end
julia> println("$x : $y")
[2 2 2] : [2 2 2]

julia> x .+= 1      # x now !== y
julia> println("$x : $y")
[3 3 3] : [2 2 2]

Because .+= promises to be element-wise, it really feels like it should behave the same as x[i] += for each element.

As I was coding, I made the mistake of using += 1, expecting it to update in-place like numpy. When it didn't I realized, oh I guess that makes sense -- it's modifying the whole Array. So I switched to .+= 1 expecting that to fix the problem and was really surprised when it didn't.

The same arguments above apply to .=, I think.

[NHDaly]

Also this is pretty jarring. It makes sense now that I've read this whole thread, but still pretty freaky:

julia> tup = Any[[1 1 1], [2 2 2]]
julia> tup[1] .*= 0     # 1. Modifies tup.
julia> z = tup[2]       #    z === tup[2]
julia> z[1] *= 0        # 2. Modifies tup.
julia> z .*= 0          # 3. Doesn't modify tup.
julia> println(tup)
Any[
[0 0 0],     # Modified in statement 1.
[0 2 2]]     # Modified in statement 2.

1. and 3. seem so similar: you are updating a matrix element-wise. It's confusing that only 1. modifies the matrix in-place, especially given that z is just another name for tup[2].

2. and 3. seem similar too, (like my previous post) 3. seems just to be a loop across all the elements, and it feels like it should do the same thing as 2. applied to each element.

[vtjnash]

especially given that z is just another name for tup[2]

this may be pedantic, but this is not really a correct statement. this sentence says that z => tup[2] => #obj#. but in actuality, the correct statement would be that "z and tup[2] were both names for the same object", or:

z      =>  #ob#
tup[2] ====^ 

this also explains why changes to tup[2] or z have no impact on the other, while changes to the contents referenced by these names (denotated #obj#) are visible, regardless of the name used to reference #obj#

[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
https://github.com/JuliaLang/julia/issues/249#issuecomment-173639124.

[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.