Automatic conversion of ^ to repeated * without applying promotion rule

I’m writing a symbolic differentiation library for Julia (see https://www.microsoft.com/en-us/research/wp-content/uploads/2016/02/main-65.pdf if you are interested) and ran into a problem with promote rules and ^. In an expression like a^7, where a is one of my symbolic numbers, I expected the promotion rule to be applied to the exponent to turn it into a symbolic constant returning a symbolic expression a^7. Instead it seems that before the promote rule is called a^7 is transformed into ((a * (a * a)) * ((a * a) * (a * a))). This is a harder expression to differentiate efficiently. My code doesn’t control when the promote rule to convert from numbers to symbolic numbers gets called.

Is there any way around this behavior?

In an expression like a^7, where a is one of my symbolic numbers, I expected the promotion rule to be applied to the exponent to turn it into a symbolic constant returning a symbolic expression a^7. Instead it seems that before the promote rule is called a^7 is transformed into ((a * (a * a)) * ((a * a) * (a * a))).

I’m not sure what this has to do with promotion rules? Not every arithmetic operation promotes its operands to a common type.

If a is a Number, then a^7 by default calls power_by_squaring(a, 7) — promotion rules are never invoked. To override it, define a method Base.:^(a::MySymbolicNumber, p::Integer) that does what you want.

If you want, you can also have it do something different (e.g. more optimized) for literal integer powers than for variable integer exponents.

3 Likes

I think you are hit by the quirk that a^7 and begin x=7; a^7 end are very different things in julia (look at Meta.@lower).

Once you successfully handled the begin x=7; a^x case, you can just Base.literal_pow(typeof(^), var::YourType, ::Val{N}) where N = var^N.

Thank you for point this out. Did you mean

a^7 and begin x=7; a^x

``

instead of

``

a^7 and begin x=7; a^7

``

Is this covered in the Julia documentation somewhere? I would have thought they were the same.

``

-brian

I don’t see how that would affect him here. a^7 lowers to Base.literal_pow(^, a, Val{7}()), but that calls ^(a, 7) by default. So overriding Base.:^ as I suggested should be sufficient.

It sounds like the problem was that he was expecting a^7 to call ^(promote(a,7)...), i.e. to promote the operands to a common type, when in fact for integer exponents it calls Base.power_by_squaring without promotion.

It’s in the documentation for ^, in addition to the source that I linked, but there is an open issue for more docs. However, as I said above, it doesn’t seem like it should affect you here.

1 Like

I have solved the problem (so far anyway) by overriding various versions of Base.:^ as suggested by stevengj. Thanks for explaining this.

It seems logically consistent though that a promotion rule would be applied to a^7, even if it does nothing for normal Number types. It would make it simpler to introduce your own types that interact well with the existing Julia Number
types and operators. It was surprising when it didn’t because promotion worked seamlessly for the other binary operators, +,/,*,- . I naively expected it would work for ^ as well. Perhaps adding this promotion rule would cause other weird bugs and that’s
why it wasn’t done.

Thanks everybody for your help.

Not at all.

In an exponentiation operations like a^7, not only is it not necessary to promote a and 7 to a common type, often it is not even sensible to do so. For example, consider the case where a is a matrix, or a dimensionful quantity, or even a string (since * is string concatenation in Julia).

Second, even in cases where you could promote operands to a common type, realize that the promotion should really be seen as a fallbackif there is no more-specialized method for a OP b, it tries promoting a and b to a common type and then calling OP. But in lots of cases there will be more specialized methods. For example, if you were computing (3+4im) + 7, you could promote to a common type as (3+4im) + (7+0im), but there is actually a specialized method for complex + real that only adds the real parts, saving an addition.

4 Likes

This is great advice, would you consider making a small PR to the manual?