Improvement needed for inplace arithmetics

[lindahua]

Following is a simple benchmark on inplace addition:

# Benchmark of inplace add
function add1!(x::Array, y::Array)
    for i in 1 : length(x)
        x[i] += y[i]
    end
    x
end

function add2!(x::Array, y::Array)
    x[:,:] += y
end

function add3!(x::Array, y::Array)
    x += y  # this is not what is desired, 
            # it creates a new array instead of updating
            # I put it here simply for performance comparison
end

x = rand(1000, 1000)
y = rand(1000, 1000)

# warming
add1!(x, y)
add2!(x, y)
add3!(x, y)

# benchmark
@time for i in 1 : 100 add1!(x, y) end  # 0.123 sec
@time for i in 1 : 100 add2!(x, y) end  # 1.036 sec
@time for i in 1 : 100 add3!(x, y) end  # 0.474 sec

The Julia expressions in add2! and add3! are way too slow than they should be. I write several functions like add!, sub!, and mul! in MLBase.jl as a stopgap. However, I believe Julia itself should provide an efficient & concise way to do this, instead of forcing people to write loops whenever they want to add a matrix to another.

What about introducing a bunch of inplace arithmetic functions such as add!, subtract!, multiply!, etc into Base?

[JeffBezanson]

Well, now that this issue is getting filed about once a month (#3217) we may finally be forced to do something about it. To me add! is not really acceptable, but this is a hard problem.

Maybe x[ ] += y could be special syntax the same way x[ ] = y is, basically implying a loop.

[lindahua]

Sure, x[] += y is much better syntax. When this is ready, corresponding functions in MLBase.jl can be deprecated.

Obviously, x[] += y does not have desired behavior at the moment.

x = [1, 2]
x[] = 10  # ==> [10, 2]
x[] += 10 # ==> [20, 2]

[timholy]

function add4!(x::Array, y::Array)
    broadcast!(.+, x, x, y)
end

is almost as fast as add1!.

[lindahua]

broadcast is great. This issue basically requests some nice syntax for this (e.g. x[] += y).

[diegozea]

What about using ! like in mutating functions, something like:

x +=! y

[kmsquire]

x[:] += y seems more natural to me than x[] += y

[JeffBezanson]

My intent with x[ ] += y was that something would go between the brackets.

[kmsquire]

Thanks for clarifying!

[lindahua]

In current behavior, x[:] += y would yield ERROR: argument dimensions must match when y has more than one dimension.

One can write x[:,:] += y. But it has two problems:

1. As shown above, this is one order of magnitude slower than the for loop. Something needs improvement here.
2. It is difficult to write generic code where y can have arbitrary number of dimensions.

[toivoh]

How about making e.g.

A[inds...] += X

invoke

updateindex!(:+, A, X, inds...)

instead of

A[inds...] = (A[inds...] + X)

I get a feeling that with Jeff's work on overloading getfield, it could become possible to dispatch on the operator symbol :+ as the first argument to updateindex! in the above example.

[toivoh]

Thinking a little more about it, I guess that the first argument would have to be the actual function (e.g. +) so that it would be possible to provide a sensible fallback method such as

updateindex!(op::Function, A, X, inds) = (A[inds...] = op(A[inds...])

That would also mean that invocations of updateindex! would see if someone is using a different binding for e.g. +. Now it seems that this would have to wait for specialization on function typed arguments, and perhaps even dispatch on specific functions as arguments.

[vtjnash]

this looks like a better candidate for "Arraymaddon" (perhaps cross-ref to #5117)

[stevengj]

This should be closed by #17623, at which point x .+= y will be fusing and in-place (along with x .+= sin.(y.^2) .+ 7 etc.).

[KristofferC]

Can you do e.g. x[1:5] += 3 with .+? Maybe that is out of scope for this issue?

Just noticed @toivoh's comment in #3424 (comment) is scarily similar to what I did in https://github.com/KristofferC/UpdateIndex.jl.

[stevengj]

@KristofferC, yes you can do x[1:5] .= 3 already and soon x[1:5] .+= 3 will work fully in-place too.

[stevengj]

Or, sorry, no, x[1:5] .+= 3 will expand to broadcast!(x -> x + 3, view(x, 1:5), x[1:5]), which is not yet fully in-place because x[1:5] still creates a copy. Once slicing creates a view, it will be fully in place.

[StefanKarpinski]

Will be addressed by #17623, which is on 0.6 already.

[KristofferC]

Once slicing creates a view, it will be fully in place.

This is likely to be never though, right? Are we giving up on a nice syntax to update just a part of an array?

[StefanKarpinski]

I still think this is addressed by the broadcast changes since the LHS of .= and co is now defined to be a view. This was a breaking change but already happened in 0.5.

[stevengj]

@StefanKarpinski, x[1:5] .+= z is defined to be x[1:5] .= x[1:5] .+ z. The left-hand-side is updated in-place (via broadcast! on a view), but x[1:5] on the right-hand side still creates a copy.

[stevengj]

But I think @view x[1:5] .+= z should work totally in-place? No, currently @view can only be used for reference expressions, so you have to do @view(x[1:5]) .+= z. I wish I could do @views x[1:5] .+= y[3:7] .* z[7:11] and have every reference expression converted to a view.

[stevengj]

I wish I could do @views x[1:5] .+= y[3:7] .* z[7:11] and have every reference expression converted to a view.

Note to future readers: nowadays, we indeed have a @views macro that works exactly like this.