in-place vectorized functions & operators

[Stefan Karpinski]

I was thinking through our discussion yesterday about wanting to have things like `X .+= Y` have the effect of doing:

for i=1:numel(X)

  X[i] += Y[i]

end

I've also been thinking about (and implementing) various in-place versions of functions for when they're applied to pre-allocated arrays, typically things with names that end in "!", like `rand!(X)`. Attempting to synthesize these two and come up with a single sane approach has led me to something I'm not entirely sure is sane but seems to do what we want. It's summarized in terse form in this gist: https://gist.github.com/c7a4e43561f2cbc59ef1.

[Stefan Karpinski]

A little explanation: the first column is what you write — a method call or an operator syntax; the second column is what it means — in terms of explicitly functional behavior or in terms of de-sugared syntax; the third column is commentary.

[Jeff Bezanson]

This is pretty good. We probably need something pretty close to this.

1. If indexing is involved, would we do this transformation:

y[I] .= -x

to

-!(sub(y,I), x)

?

2. I assume .= only works with certain operators on the right, so "y
.= x" can't be used to copy?

3. It's too bad these are redundant with a[:,:] = b+c except with
respect to performance. That almost makes me think it's not a good
thing to be doing. Maybe this can be a feature of indexed assignment
somehow?

[Viral Shah]

This operator sounds like APL.

-viral

[Alan Edelman]

i kind of like this

[Stefan Karpinski]

This is pretty good. We probably need something pretty close to this.

Yeah, I think this is sort of the logical conclusion if we decide this is a good idea.

 

1. If indexing is involved, would we do this transformation:

y[I] .= -x

to

-!(sub(y,I), x)

?

Seems sensible. I can't see why not. Is there an issue with doing this?

2. I assume .= only works with certain operators on the right, so "y
.= x" can't be used to copy?

Just for cuteness, I was thinking that we might as well make `Y .= X` mean `copy_to(Y,X)`. I'm also slightly tempted to rename map_to and copy_to to map! and copy! but I'm not sure that improves clarity in any way. Maybe having map! and copy! be aliases for map_to and copy_to would be a good thing.

 

3. It's too bad these are redundant with a[:,:] = b+c except with
respect to performance. That almost makes me think it's not a good
thing to be doing. Maybe this can be a feature of indexed assignment
somehow?

That could be another way to go. Although indexed assign is ugly to write — and less good for polymorphism because you have to use the correct number of : for the rank of the array. Or use linear indexing on all of them: `a[:] = b[:] + c[:]`. But that seems far less nice to write than `a .= b + c`. If we do implement this, we could automatically stick a . in front of any un-dotted assignment operator following an indexing expression. So `a[:,:] = b + c` actually means `a[:,:] .= b + c` which in turn actually means `+!(sub(a,1:size(a,1),1:size(a,2)),b,c)`. That way indexed assignment automatically gets better performance with the same effect.

[Stefan Karpinski]

Yeah, I had the same thought while working this out. However, I think the upshot, which is that you can write things like

X .+= 1

Y .+= X

Y .= X1 * X2

and have it work efficiently without temporary arrays, is fairly sane and desirable. The internal syntax is really just so that we can define methods that implement this stuff — and will mostly actually be done with macros that you give the core operator name to. The user would never have to write anything like `*!(Y,X1,X2)`. That's just what the AST looks like. And the parser and everything else doesn't give a shit if function names look like line noise.

On Fri, Dec 9, 2011 at 11:23 PM, Viral Shah <vi...@mayin.org> wrote:

[Stefan Karpinski]

The f! naming consistency that leads to the APL-like operator names is a good thing because it lets the same macro do more work. Here's what, e.g. the vectorize_unary macro might look like:

macro vectorize_unary(f)

  f! = symbol("$(f)!")

  quote

    function ($f!)(Y::Array, X::Array)

      for i = 1:numel(X)

        X[i] = ($f)(Y[i])

      end

      return Y

    end

    ($f!)(X::Array) = ($f!)(X,X)

    ($f)(X::Array) = ($f!)(similar(X,map_type(f,X)),X)

  end

end

The consistent naming of functions and operators means you can call this equally well as `@vectorize_unary log` or as `@vectorize_unary (-)` and have it work. Sure you get weirdly named operations like `-!` but since you never have to use that name with the proposed syntax, who cares? For non-operators, you get functions like `log!` which can be used to apply the log function in place across an array, replacing each value with its logarithm, which is handy.

[Viral Shah]

What if the syntax were 

X! += 1

-viral

[Stefan Karpinski]

We can't do that because `X!` is a legitimate name.

[Stefan Karpinski]

The `.=` name makes a lot of sense because does vectorized assignment.

[Gustavo Goretkin]

Can it be (or is it the case now) that the compiler converts something like

a = exp(a)

into in-place assignment as opposed to making an intermediate array?

[Viral Shah]

It will make a new copy as of now. The suggestion (as in the trailing mail) is to use the following syntax:

a .= exp(a)

The thing I don't like about the . syntax is that we now have two ways of doing in-place operations - namely by suffixing ! after functions or by prefixing . in front of operators. I would have liked it if there were one way.

Also, for things like a = exp(a), it would be nice if the compiler can just do it automatically, since I am already telling the compiler that by using the same name in the input and output.

Perhaps even fancier usage could be where the compiler can do analysis on all the temporary creation and decide in which cases it needs to allocate new stuff and which stuff it can overwrite the input. Thus, in some cases, the compiler would be able to call foo! automatically instead of foo. It may also mean that the compiler needs to differentiate between mutable and immutable types.

-viral

[Stefan Karpinski]

This could be done but the compiler needs to know a whole lot to safely make this transformation:

1. that the incoming `a` references an array
2. that `a` is the only reference to that array
3. that the outgoing `a` will also be an array with the same size and element type
4. that `exp!` exists and does the same thing as `exp` but in-place

Being able to know all of this is sufficiently hard that it's probably almost never going to pay off. It seems better to just give the programmer explicit power to do in-place operations — hence functions like `exp!`.

[Stefan Karpinski]

The thing I don't like about the . syntax is that we now have two ways of doing in-place operations - namely by suffixing !  after functions or by prefixing . in front of operators. I would have liked it if there were one way.

With my proposal `Y .= X1 + X2` and such is just syntax for `+!(Y,X1,X2)` so there really is just one way. We could write something like `Y =! X1 + X2` but that doesn't seem better to me — it just looks like a weird not-equals symbol. Is the issue with the choice of `.=` as an operator or with the whole idea of having syntax for something like in-place addition?

Also, for things like a = exp(a), it would be nice if the compiler can just do it automatically, since I am already telling the compiler that by using the same name in the input and output.

Perhaps even fancier usage could be where the compiler can do analysis on all the temporary creation and decide in which cases it needs to allocate new stuff and which stuff it can overwrite the input. Thus, in some cases, the compiler would be able to call foo! automatically instead of foo. It may also mean that the compiler needs to differentiate between mutable and immutable types.

As my other email pointed out, this is very hard to do. It also seems to me to go against the grain of Julia as a language: we tend toward designing things so that you can write code that can turned into an efficient implementation in a fairly obvious way, rather than expecting the compiler to do insanely clever things that might transform something that can't obviously be done efficiently into an efficient implementation. This would be introducing an unprecedented expectation of high-level compiler intelligence. I feel like languages that do that either a) never get the promised level of intelligence and performance, or b) do get it, but at the expense that it is no longer very clear what the system is doing and why any particular code might be slow or fast. By having explicitly in-place operations, we knew exactly what is going on and can decide for ourselves whether or not an operation should be done in-place or not. I much prefer that to some promise of a compiler that might someday be able to generate fast, in-place code for me.

[Viral Shah]

> With my proposal `Y .= X1 + X2` and such is just syntax for `+!(Y,X1,X2)` so there really is just one way. We could write something like `Y =! X1 + X2` but that doesn't seem better to me — it just looks like a weird not-equals symbol. Is the issue with the choice of `.=` as an operator or with the whole idea of having syntax for something like in-place addition?

I love the idea of in-place operations. My objection is just stylistic, and about the choice of the ".=" operator. I agree that =! looks like a weird not-equals symbol. I don't yet have a better solution, but I feel that the whole syntax could be a little more unified - the current approach with ! and .= feels like a bit of a hack.

Do we need to think through something at a more fundamental level?

Also, what does "a .*= b" mean? Is it "a = a .* b" or in-place "*!(a, a, b)"? Should we then have "a ..*= b" for in-place ".*".

-viral

[Stefan Karpinski]

I love the idea of in-place operations. My objection is just stylistic, and  about the choice of the ".=" operator. I agree that =! looks like a weird not-equals symbol. I don't yet have a better solution, but I feel that the whole syntax could be a little more unified - the current approach with ! and .= feels like a bit of a hack.

Do we need to think through something at a more fundamental level?

Also, what does "a .*= b" mean? Is it "a = a .* b" or in-place "*!(a, a, b)"? Should we then have "a ..*= b" for in-place ".*".

I actually like the .= operator since I think it's pretty intuitive what it means: vectorized assignment into an array. However, I could live without the mutating versions like `Y .+= X` since that can be written as `Y .= Y + X` which has a clear meaning and doesn't introduce all the APL-like weird operators. And you're absolutely right that `a .*= b` is ambiguous.

More fundamentally, one thing this doesn't address is if you write something like

Y .= X1 + X2 - X3

Since the RHS has mixed operations, you can't have this handled be any single in-place operation implementation. In the current proposal, it would mean

+!(Y, X2, X2-X3)

which requires a temporary array for the result of `X2-X3`, defeating the purpose of vectorized assignment. The obvious desirable meaning is

for i = 1:numel(Y)

  Y[i] = X1[i] + X2[i] - X3[i]

end

So maybe that's what it should mean: use the RHS as an expression template and do vectorized assignment to each slot of the LHS, substituting for each variable, the appropriate element of that variable. You could even do things like

w .= f(v)

where f is a vector of functions and it would work. Of course, if f was just a function, then it would get applied every time.

[Viral Shah]

Often, the expressions you describe show up with indexing:

Y = 0.5*X[I] + 2*X + 0.5*X[J]

You are essentially using the .= syntax here as a compiler hint for devectorization. I think the right thing is for the compiler to detect these cases and devectorize. I believe Fortran has had this for the longest time. The cases where I and J are ranges can be optimized quite aggressively.

One usually has to decide which cases should be devectorized, since sometimes the vectorized versions can be faster, especially if you can use the vector registers. We don't yet have vectorization, and we I'm not sure we can use all the SSE and AVX stuff yet anyways.

-viral

[Stefan Karpinski]

This is not a vectorization hint. It's saying something completely different. Writing `Y = X` means bind the value bound to X  to the variable Y, implicitly declaring the variable Y if it hasn't already been declared. Using .= has nothing to do with binding anything to the variable Y or declaring it — it's a way of putting specific data into a pre-existing array that Y is already bound to. If Y doesn't already exist or isn't an array, or is the wrong size, then it's an error. In particular, `Y .= X` nevers declares a new variable because `Y` already needs to exist in order for there to be a place to put the results.

Doing automatic vectorization sounds nice in theory, but it's extremely hard for reasons I've already pointed out that you haven't at all addressed. In a language like Fortran, it's possible because the language is completely statically typed (zero polymorphism), with a fairly simple type system (no user-defined types even), and operations defined only on those built-in types. In Julia none of these things are true and it makes something like automatic vectorization very, very hard, if not impossible.

[Jeff Bezanson]

Maybe we could somehow consider .= to be syntax for a loop, basically asking to have a for loop around the statement. This could be helpful in cases where vectorization fails to generalize, like a[x]=b where a is a hash table and x is a vector. a[x].=b would insert many keys, while the first one would use the whole vector as one key.

[Stefan Karpinski]

Yeah, I like this general direction. Using .= or equivalent for dgemm is actually a bit weird because it's not really a vectorized operation — it's really a special function that writes to one of its arguments. So not handling dgemm with this syntax seems like a good thing to me. Another hint that it's not great for dgemm is that it doesn't generalize to n arguments.

[Stefan Karpinski]

Maybe we could somehow consider .= to be syntax for a loop, basically asking to have a for loop around the statement. This could be helpful in cases where vectorization fails to generalize, like a[x]=b where a is a hash table and x is a vector. a[x].=b would insert many keys, while the first one would use the whole vector as one key.

So this would mean something like

for i = numel(x)

  a[x[i]] = b[i]

end

? Seems reasonable. Hard to figure out the right rules for this though, but could be a great way to tersely write loops.