[for discussion] an abstract symbolic type and chain rule code

[tshort]

John,

This is based on the julia-dev discussion here:

https://groups.google.com/d/msg/julia-dev/3rqlCOf2TCw/cpRk1e7fKrMJ

This pull request includes:

• A new abstract type for symbolic variables, so packages that manipulate expressions can use placeholder types as well as just symbols to represent variables. Here is an example:
• An adaptation of Miles Lubin's chain rule code from here:

https://github.com/IainNZ/NLTester/blob/master/julia/nlp.jl#L81

Here is an example:

julia> x = BasicVariable(:x)
x

julia> y = BasicVariable(:y)
y

julia> x^3 + sin(x^2 + 1)
SymbolicExpression :(+(^($(x),3),sin(+(^($(x),2),1))))

julia> differentiate(x^3 + sin(x^2 + 1), x)
:(+(*(3,^($(x),2)),*(cos(+(^($(x),2),1)),*(2,$(x)))))

julia> chainRule(x^3 + sin(x^2 + 1), x)
:(+(*(3,1,^($(x),2)),*(*(2,1,^($(x),1)),cos(+(^($(x),2),1)))))

Here are some items for discussion:

• In the big picture, is this the right place for this, and is this the direction we want to go?
• In my adaptation of Miles' code, I split up chainRule to operate by the symbol representing the function, so there's one function for :+ another for :cos and so on. It seems like it works well. It relies on types being able to be parametrized by a symbol.
• chainRule and differentiate seem redundant. We should probably go with one or the other. I like the name differentiate better.
• chainRule doesn't have a simplifying step.
• In adapting my symbolic code from Sims, I came across some issues with John's code. The biggest is that comparison operations like x == 3 return an expression if x is a symbolic variable. I ended up using isequal when I wanted to compare the actual contents rather than return an expression.
• I did not add the capability to do :x^2 + cos(:z). That would be possible, but that'd be classified as monkey patching unless we got buy-in from Julia-core. Because we can't operate on Expr's, I added a SymbolicExpression type.

[tshort]

Pinging @mlubin, @IainNZ, @lindahua, @stevengj, and @dmbates.

[mlubin]

There's something strange about using operator overloading to build up expression trees when the Julia parser can already give you them through macros. While operator overloading is the only possible approach in many languages, it's especially slow and creates lots of temporary objects. We can completely avoid it by using some Julia magic:

function processExpr(x::Expr)
    if x.head == :call
        quoted = Expr(:quote,x.args[1])
        code = :(Expr(:call,$quoted))
        for y in x.args[2:end]
            push!(code.args,processExpr(y))
        end
        return code
    else
        return x
    end
end

processExpr(x::Any) = x

macro sexpr(x)
    esc(processExpr(x))
end

Then

julia> x = BasicVariable(:x)
x

julia> y = BasicVariable(:y)
y

julia> @sexpr x^3 + sin(x^2 + 1)
:(+(^($(x),3),sin(+(^($(x),2),1))))

This is a simplified version of what's implemented in https://github.com/IainNZ/NLTester/blob/master/julia/nlp.jl. (Note I didn't actually test the above code, and the escaping behavior isn't precisely correct.) The macro-fu is a bit heavy but what the macro does is generate code that regenerates the input expression with the values of all symbols spliced in. This approach will give improvements of orders of magnitude over using operator overloading.

The current approach also defines the complete set of operator overloads for anything that inherits from Symbolic. It could be reasonable in some cases to use operator overloading when building these expressions isn't a performance bottleneck, but it's not essential for using the chain rule, and there may be applications (like mine) that don't want the operator overloads defined at all.

[tshort]

Good points, Miles. I'll play around with implementing that approach. In my
application, generating the expressions isn't time critical (compared to
the solution phase), but your approach does have some advantages,
especially the simplicity.

On Fri, Mar 29, 2013 at 6:07 PM, Miles Lubin notifications@github.comwrote:

There's something strange about using operator overloading to build up
expression trees when the Julia parser can already give you them through
macros. While this is the only possible approach in many languages, it's
especially slow and creates lots of temporary objects. We can completely
avoid this by using some Julia magic:

function processExpr(x::Expr)
if x.head == :call
quoted = Expr(:quote,x.args[1])
code = :(Expr(:call,$quoted))
for y in x.args[2:end]
push!(code.args,processExpr(y))
end
return code
else
return x
endend
processExpr(x::Any) = x
macro sexpr(x)
processExpr(x)end

Then

julia> x = BasicVariable(:x)x
julia> y = BasicVariable(:y)y
julia> @sexpr x^3 + sin(x^2 + 1):(+(^($(x),3),sin(+(^($(x),2),1))))

This is a simplified version of what's implemented in
https://github.com/IainNZ/NLTester/blob/master/julia/nlp.jl. (Note I
didn't actually test the above code.) The macro-fu is a bit heavy but what
the macro does is generate code that regenerates the input expression with
the values of all symbols spliced in. This approach will give improvements
of orders of magnitude over using operator overloading.

The current approach also defines the complete set of operator overloads
for anything that inherits from Symbolic. It could be reasonable in some
cases to use operator overloading when building these expressions isn't a
performance bottleneck, but it's not essential for using the chain rule,
and there may be applications (like mine) that don't want the operator
overloads defined at all.

—
Reply to this email directly or view it on GitHubhttps://github.com//pull/7#issuecomment-15663031
.

[tshort]

Miles, your @sexpr approach does look like the way to go. Everything is simpler.

[lindahua]

In Devectorize.jl, I implemented the texpr function that composes an instance of TExpr from Julia's expression. While TExpr was used as the basis for generating devectorized codes, and SymbolicExpression here is to be used to get derivatives, I think the part of translating Julia Expression to symbolic expression tree can be implemented in a similar way.

I would be nice to extract this common part to a package (e.g. SymbolicExpressions.jl) and merge the efforts there, so as to provide support to other packages that rely on symbolic expressions. (Now, I have known five packages rely on abstract expressions in one or the other way: Calculus, NLP, Devectorize, Sims, and DeMat)

Actually, @StefanKarpinski once mentioned this idea in julia-dev (https://groups.google.com/forum/?fromgroups=#!searchin/julia-dev/common$20infrastructure/julia-dev/xIrPPMY-6qM/h1vhxjF7hwAJ)

cc @ViralBShah

[lindahua]

Is vector calculus on the plan? If vector calculus is included, this may get real use in many practical problems.

For example, in a regression problem, I would write
f(w' x) + (w'w) / 2
Instead of writing

f(w1 * x1 + w2 * x2 + .... w1000 * x1000) + ...

[StefanKarpinski]

I would be happy with such a package being called simply Symbolic. Seems like an excellent piece of foundational infrastructure for other packages to build on.

[johnmyleswhite]

I'll be away all weekend. I'm happy to look over this next week and merge a finished product or introduce a dependency on a new Symbolic package.

[mlubin]

@StefanKarpinski: What's the right way to escape expressions in processExpr? Is it enough to replace the Expr call with Base.Expr so that it's reasonable to escape the whole thing?

@tshort Note you can also make multiple frontend macros like:

macro differentiate(x,wrt)
    :(differentiate($(esc(processExpr(x))),$(esc(wrt)))
end

So that you can do

@differentiate(x^3 + sin(x^2 + 1),x)

[tshort]

John's differentiation code was more complete than that in chainRule. I'm trying to merge them, but it's a work in progress.

[tshort]

@lindahua, what do you think should be in here for a SymbolicExpression? I originally had the following in here:

type SymbolicExpression <: Symbolic
    ex::Expr
end
sexpr(hd::Symbol, args::ANY...) = SymbolicExpression(Expr(hd, args...))
typealias SExpr Union(Expr, SymbolicExpression)

Is the idea that you could use SExpr anywhere you would use Expr for symbolic operations? Or, are you after something different?

[tshort]

John, what's in here is now more or less ready for review. I basically ended up reworking differentiate based on your code and what Miles did. I reworked simplify with the same approach.

Of course, more discussion is welcome on other symbolic features to include.

[johnmyleswhite]

Ok. This all looks good to me. I'll make a more thorough review tomorrow morning and then merge this.

[mlubin]

Should this be if xp != 0? There's also an instance of this for the bessel functions.

[tshort]

Good catch, Miles. Thanks for the review. I'll fix this tomorrow.

[mlubin]

I made a couple inline comments, but otherwise looks great.

[johnmyleswhite]

I'm happy enough with this that I'm going to merge it. Any more changes you'd like to make?

[tshort]

I don't have more immediate changes, John.

Note that I did remove your derivative(ex:Expr, target) code. It seemed
broken (or at least very limiting) to me. It might be good to have
machinery that will convert expressions to functions, including derivatives
and other expressions. It might take some creativity to do this flexibly.

On Mon, Apr 1, 2013 at 12:32 PM, John Myles White
notifications@github.comwrote:

I'm happy enough with this that I'm going to merge it. Any more changes
you'd like to make?

—
Reply to this email directly or view it on GitHubhttps://github.com//pull/7#issuecomment-15723133
.

[johnmyleswhite]

It was very limited. This new code is a good impetus to rethinking what I had done.