[sympy] partial derivatives, please explain

[janwillem]

I need some explanation on the workings of SymPy. As an example the
following script:
import sympy
X, F, B = sympy.symbols('XFB')
Y = X / F - B #eqn 1
DY = sympy.Matrix(sympy.diff(Y, (X, F, B))).T
print DY

I had expected (eqn 2): [1/F, -X/F**2, -1]
But get: [D(-B + X/F, X), D(-B + X/F, F), D(-B + X/F, B)]

Not after doing a trick of which I cannot remember why I tried it, I
get the desired result
DY = DY.subs({X:X, F:F, B:B})
From the doc I had thought that DY.doin() would work but that gives
"raise AttributeError()".
So there is obviously something I do not understand, please some help
Janwillem

--
You received this message because you are subscribed to the Google Groups "sympy" group.
To post to this group, send email to sy...@googlegroups.com.
To unsubscribe from this group, send email to sympy+un...@googlegroups.com.
For more options, visit this group at http://groups.google.com/group/sympy?hl=en.

[Bastian Weber]

Hi Jan,

janwillem wrote:
> I need some explanation on the workings of SymPy. As an example the
> following script:
> import sympy
> X, F, B = sympy.symbols('XFB')
> Y = X / F - B #eqn 1
> DY = sympy.Matrix(sympy.diff(Y, (X, F, B))).T
> print DY
>
> I had expected (eqn 2): [1/F, -X/F**2, -1]
> But get: [D(-B + X/F, X), D(-B + X/F, F), D(-B + X/F, B)]
>
> Not after doing a trick of which I cannot remember why I tried it, I
> get the desired result
> DY = DY.subs({X:X, F:F, B:B})
>>From the doc I had thought that DY.doin() would work but that gives
> "raise AttributeError()".

> So there is obviously something I do not understand, please some help
> Janwillem
>

.doit is a method of derivative-objects ( like D(-B + X/F, X) ) and not
of Matrix-objects (like DY)

What you could do is:

DY.applyfunc(lambda arg: arg.doit())

But maybe thats not very intuitve (at least as long as you are not
familiar with pythons lambda functions).

An alternative way to calculate your desired result would be

DY = Matrix([Y]).jacobian([X,F,B])


The reason why your subs-trick works is out of my knowledge.


Regards,
Bastian.

[Sebastian]

On 05/04/2010 10:32 AM, janwillem wrote:
> I need some explanation on the workings of SymPy. As an example the
> following script:
> import sympy
> X, F, B = sympy.symbols('XFB')
> Y = X / F - B #eqn 1
> DY = sympy.Matrix(sympy.diff(Y, (X, F, B))).T
>

Do sympy.diff(Y, (X, F, B), evaluate=True) instead for now. But you're
right - this shouldn't be necessary. Somehow the vectorization of the
diff method seems to fail and set evaluate to False.

> print DY
>
> I had expected (eqn 2): [1/F, -X/F**2, -1]
> But get: [D(-B + X/F, X), D(-B + X/F, F), D(-B + X/F, B)]
>
> Not after doing a trick of which I cannot remember why I tried it, I
> get the desired result
> DY = DY.subs({X:X, F:F, B:B})
>

This works because subs evaluates the expression in the end (which
wasn't done before).

> From the doc I had thought that DY.doin() would work but that gives
> "raise AttributeError()".
>

First, the name of the method is doit. Second, DY is a python tuple
which doesn't have this method.

> So there is obviously something I do not understand, please some help
> Janwillem
>
>

Sebastian

[janwillem]

Thanks for the explanation. I moved the doit()s to later in the script
where everything boils down to scalars and that works.

I had tried jacobean but sometimes that gives me AttributeError. Here
is the real thing I was working on:

def lpu_nonlin(Y, X, C):
"""
LPU for non-linear measurement models according to
Giovanni Mana and Francesca Pennecchi
Metrologia 44 (2007) 246-251
Inputs:
Y measurement model e.g. Y = (X[0] / X[1]) - X[2]
X vector of input quantities
C covariance matrix
Returns symbolic equations for:
expectation using LPU with 3rd degree Taylor series
variance using LPU with 3rd degree Taylor series
expectation using LPU with 1st degree Taylor series
variance using LPU with 1st degree Taylor series
"""
#Jacobean vector
J = sympy.Matrix([Y]).jacobian(X).T
#Hessian matrix
H = sympy.hessian(Y, X)

#Components of Mana & Pennacchi equations 2 and 3
JTC = J.T * C
JTCJ = (J.T * C * J)[0, 0]
HC = H * C
TrHC = HC.trace()
HC2 = HC**2
TrHC2 = HC2.trace()
dTrHC = sympy.Matrix([sympy.diff(TrHC, X)]).T
#dTrHC = sympy.Matrix([TrHC]).jacobean(X).T

#Expectation, Mana & Pennacchi eqn 2
exp3_y = (Y + TrHC / 2).doit()
#variance, Mana & Pennacchi eqn 3
var3_y = (JTCJ + TrHC2/2 + (JTC * dTrHC)[0, 0]).doit()

#according to GUM framework LPU
exp1_y = Y
var1_y = JTCJ.doit()
return exp3_y, var3_y, exp1_y, var1_y

For getting J I already changed from diff to jacobean as you suggested
but doing the same a few lines further on with TrHC fails (the
commented version, TrHC is a scalar function of X and C).

Why can that be? And why is the syntax of hessian and jacobean so
different?
Thanks again, Janwillem


On May 4, 11:55 am, Bastian Weber <bastian.we...@gmx-topmail.de>
wrote:

[Bastian Weber]

After a short look I guess its a simple typo:
you wrote 'jacobean' instead of 'jacobian'

> And why is the syntax of hessian and jacobean so
> different?

My guess: the hessian is typically calculated for a scalar function
while the jacobian is usually calculated for a vector field (i.e a
n-by-1 matrix or its transpose). Therefore 'jacobian' is a method of the
Matrix-class while hessian is a function that lives directly in the
sympy namespace.

As I said: Thats just a guess.
Regards. Bastian.

[Aaron S. Meurer]

On May 4, 2010, at 4:25 AM, Sebastian wrote:

> On 05/04/2010 10:32 AM, janwillem wrote:
>> I need some explanation on the workings of SymPy. As an example the
>> following script:
>> import sympy
>> X, F, B = sympy.symbols('XFB')
>> Y = X / F - B #eqn 1
>> DY = sympy.Matrix(sympy.diff(Y, (X, F, B))).T
>>
> Do sympy.diff(Y, (X, F, B), evaluate=True) instead for now. But you're
> right - this shouldn't be necessary. Somehow the vectorization of the
> diff method seems to fail and set evaluate to False.

I created issue 1929 for this.

Aaron Meurer

[b45ch1]

May I ask in what context you require the matrix derivatives.
Depending on what you are after an Algorithmic Differentiation (AD)
tool could also be a viable alternative.

With an AD tool you could differentiate functions (i.e. compute
gradient, Jacobian, Hessian, higher order tensors) very easily and
efficiently.

An example of what functions are possible:

def f(A,x):
for n in range(30):
y = dot(x.T,dot(A,x))
A = inv(A) - dot(x,x.T) * y

return trace(A)


best regards,
Sebastian

[janwillem]

Yes a typo. Dyslexia doesn't always come in handy.
I think that using jacobian is the elegant solution.
The not evaluating the diff is essentially only an issue when I wanted
to texify the intermediate results of a derivation and compare with
the original publication.
Now everything works and after lambdify my script works like a charm.
Thanks, Janwillem

[Bastian Weber]

Hi,

is there any (easy) way to get an idea of how complex a given expression is?

I saw that "implementing complexity measures" was among the ideas for
the GSoC2009. I dont know how sophisticated this was planned to be.

My naive approach is the following:

def depth(expr):
if expr.is_Atom:
return 0
else:
return max(map(depth, expr.args))+1


(I have a tree with many levels of subbranches in mind.)

Any comments?

The background is: I am trying to investigate why expand() hangs or
needs very long time if applied to determinants of polynomial matrices
(just 4-by-4).


Regards,
Bastian

[Aaron S. Meurer]

You might want to look at the compare_ode_sol() (and the key equivalent rewrite I have at http://github.com/asmeurer/sympy/commit/2c5167990473ca4079f334214ffe0c898b18cbda). It is specifically for solutions given from dsolve(), but the heuristic involves seeing if the expression contains unevaluated Integral's, is or isn't solved for f(x), and finally, the length of str(expr).

len(str(expr)) is probably the simplest way to check this. I never considered looking at the depth of the expression, though.

By the way, you can probably make that operation faster first, by using the Poly class, and second, by somehow making det() do intermediate simplifications (things should actually already be better in this respect with the new polys).

Aaron Meurer

[Ondrej Certik]

On Tue, May 4, 2010 at 2:28 PM, Bastian Weber
<bastia...@gmx-topmail.de> wrote:
> Hi,
>
> is there any (easy) way to get an idea of how complex a given expression is?
>
> I saw that "implementing complexity measures" was among the ideas for
> the GSoC2009. I dont know how sophisticated this was planned to be.
>
> My naive approach is the following:
>
> def depth(expr):
>    if expr.is_Atom:
>        return 0
>    else:
>        return max(map(depth, expr.args))+1
>
>
> (I have a tree with many levels of subbranches in mind.)
>
> Any comments?

One undocumented function is count_ops:

In [1]: e = x**2+y+z

In [2]: e.count_ops(False)
Out[2]: 3

In [3]: e = x**2+y+z+sin(z)

In [4]: e.count_ops(False)
Out[4]: 5

In [5]: e = x**2+y+z+sin(z)**2

In [6]: e.count_ops(False)
Out[6]: 6

See also the tests in sympy/core/tests/test_count_ops.py

I think it's similar to your function too. It'd be cool to maybe write
some documentation and example doctests for this, and/or integrate it
with your depth() function.

Ondrej

[asmeurer]

I have a fix for the unevaluated derivatives bug at
http://code.google.com/p/sympy/issues/detail?id=1929.

In my branch:

In [1]: import sympy

In [2]: X, F, B = sympy.symbols('XFB')

In [3]: Y = X / F - B #eqn 1

In [4]: DY = sympy.Matrix(sympy.diff(Y, (X, F, B))).T

In [5]: print DY
[1/F, -X/F**2, -1]

If either of you want to review this, I can push this in.

Aaron Meurer

> > For more options, visit this group athttp://groups.google.com/group/sympy?hl=en.

>
> --
> You received this message because you are subscribed to the Google Groups "sympy" group.
> To post to this group, send email to sy...@googlegroups.com.
> To unsubscribe from this group, send email to sympy+un...@googlegroups.com.

> For more options, visit this group athttp://groups.google.com/group/sympy?hl=en.

[Bastian Weber]

Aaron S. Meurer wrote:
> You might want to look at the compare_ode_sol() (and the key
> equivalent rewrite I have at
> http://github.com/asmeurer/sympy/commit/2c5167990473ca4079f334214ffe0c898b18cbda).
> It is specifically for solutions given from dsolve(), but the
> heuristic involves seeing if the expression contains unevaluated
> Integral's, is or isn't solved for f(x), and finally, the length of
> str(expr).
>
> len(str(expr)) is probably the simplest way to check this. I never
> considered looking at the depth of the expression, though.
>

Currently I have just basic arithmetics, so I guess the heuristic would
reduce to len(str(expr)). That came to my mind earlier but I tried to
avoid that because computing str(expr) seems to be rather costly. My
experience is: printing the str representation of large expressions
(such as determinants) takes much more time than calculating them.

Just a (non representative) example: Computing that an expression has
"depth" 23 take less than a second. But computing str(expr) took too
long for my patience. (i.e. more than about 30 seconds)


For curiosity I did the following:

I used the function:


def count_branches(ex):
if ex.is_Atom:
return 1

else:
return sum(map(count_branches, ex.args))


to compute how many lowest level branches the tree has.
Result for my expression is 157130. (!)

The expression has many Rational-instances as atoms, so the average
stringlength of an atom in my case is 20. This means more than 2 MB of
str representation. When I relate this to the computation time of str(.)
of moderate expressions I get the feeling that this will take very long.

> By the way, you can probably make that operation faster first, by
> using the Poly class, and second, by somehow making det() do
> intermediate simplifications (things should actually already be
> better in this respect with the new polys).

Using the Poly class is a very good hint! In my first test speedup is
very notable.

But I really think intermediate simplifications are the key.
Unfortunately I have no idea how to do that.

But it is not that important at the moment.

I will concentrate on how I can avoid the complexity from the beginning.

[Aaron S. Meurer]

So the problem is large rational coefficients. You might want to try setting the environment variable SYMPY_GROUND_TYPES=gmpy (you will need to install gmpy). Assuming you are working in master with the new polys, this should make things faster, especially if you are using Poly()'s. ZZ.dtype will tell you the ground type. Integer is SymPy (slowest), int is Python (next slowest), and mpz is gmpy (the fastest).

>
>
>
>
>> By the way, you can probably make that operation faster first, by
>> using the Poly class, and second, by somehow making det() do
>> intermediate simplifications (things should actually already be
>> better in this respect with the new polys).
>
> Using the Poly class is a very good hint! In my first test speedup is
> very notable.
>
> But I really think intermediate simplifications are the key.
> Unfortunately I have no idea how to do that.

You also might try the non-default determinant method, i.e. det(method="berkowitz"). Actually, for me that method returns less simple results (more factored) than the default one, but you might do it anyway just to see what you get.

Aaron Meurer

[Bastian Weber]

>
> One undocumented function is count_ops:
>
> In [1]: e = x**2+y+z
>
> In [2]: e.count_ops(False)
> Out[2]: 3
>
> In [3]: e = x**2+y+z+sin(z)
>
> In [4]: e.count_ops(False)
> Out[4]: 5
>
> In [5]: e = x**2+y+z+sin(z)**2
>
> In [6]: e.count_ops(False)
> Out[6]: 6

Thats very cool.
Its what I was looking for.

>
> See also the tests in sympy/core/tests/test_count_ops.py
>
> I think it's similar to your function too. It'd be cool to maybe write
> some documentation and example doctests for this, and/or integrate it
> with your depth() function.

cout_ops(False) does almost exactly the same as my count_branches from
my last mail. (I should check for new mail before sending a mail to the
list...)


The main difference is that count_ops seems to be faster.

I dont see any useful possibility to integrate the depth function with
count_ops and I doubt that anybody might want to use the former.

Concerning the docs / doctests for count_ops I do not understand the
current situation:


in ipython

Symbol.count_ops.__doc__ == None

evaluates to True

but on the other hand, if I go to the source file, I can read the following:

""" Return the number of operations in expressions.

Examples:
>>> from sympy.abc import a, b, x
>>> from sympy import sin
>>> (1+a+b**2).count_ops()
POW + 2*ADD
>>> (sin(x)*x+sin(x)**2).count_ops()
2 + ADD + MUL + POW

"""


__version__ = '0.7.0-git'

I am confused about this.

Bastian.

[Ondrej Certik]

Ah, I see. It's because it's only documented in Basic, but you use it
for Symbol.

The count_ops in all classes except Basic should be renamed to _count_ops.

The docstring in Basic should also include the symbolic=False option
(that should probably be the default).

Ondrej

[Aaron S. Meurer]

Symbol should be inheriting the method from Basic though, and it's docstring. For example,

In [8]: Symbol.expand.__doc__
Out[8]:

Expand an expression using hints.

See the docstring in function.expand for more information.

In [11]: Expr.expand.__doc__ is Symbol.expand.__doc__
Out[11]: True

Maybe it is because it comes from Basic and not Expr.

Aaron Meurer

[Ondrej Certik]

Yes, but Symbol doesn't reimplement expand(), but the count_ops are
reimplemented in couple places.

[smichr]

commit 1923 (waiting for review) makes symbolic (renamed visual) False
by default. count_ops is much faster than computing len(str(expr)).
With count_ops you are just adding up integers.

The reason that different objects implement count_ops is for the
visual feature:

>>> count_ops(3+4*x,visual=1)
ADD + MUL
>>>

Bastian: regarding the hang in expansion...please try it on branch
1766 at my smichr account at github. Ronan timed it recently and it
looks like it's doing well.

[Vinzent Steinberg]

On 4 Mai, 23:58, "Aaron S. Meurer" <asmeu...@gmail.com> wrote:
> So the problem is large rational coefficients.  You might want to try setting the environment variable SYMPY_GROUND_TYPES=gmpy (you will need to install gmpy).   Assuming you are working in master with the new polys, this should make things faster, especially if you are using Poly()'s.  ZZ.dtype will tell you the ground type.  Integer is SymPy (slowest), int is Python (next slowest), and mpz is gmpy (the fastest).  

I think gmpy is only faster for large integers.

Vinzent

[Bastian Weber]

Hi at all,

Aaron S. Meurer wrote:
>
> So the problem is large rational coefficients. You might want to try
> setting the environment variable SYMPY_GROUND_TYPES=gmpy (you will
> need to install gmpy). Assuming you are working in master with the
> new polys, this should make things faster, especially if you are
> using Poly()'s. ZZ.dtype will tell you the ground type. Integer is
> SymPy (slowest), int is Python (next slowest), and mpz is gmpy (the
> fastest).

Thank you for the suggestion. I'll keep that in mind. But currently I
dont have the time for this. Because for the current task this would not
help. The reason is: at the end I want to cope with symbolic
coefficients. I just used random integers for testing. These result in
rational coefficients because I apply the sympy.gcdex algorithm
repeatedly. I think I have to concentrate on this point. If I could
simplify the results here (and I think I have some degrees of freedom),
I can avoid the huge expressions in the coefficients.

The overall goal is to compute an unimodular polynomial matrix U
(unimodular means U.det() is a degree zero polynomial) which transforms
a given Matrix in hermitian normal form.
To be more precise: The computation of U is easy but I also need its
inverse V. The computation of V still is quite fast. But the Problem is:
the elements of V seem to be very complex. As stated earlier count_ops
results in ca. 1.5e5 (for one element). Printing an element of V takes
too long. (But it would not give much insight I guess).


While playing around with those matrices: I came across the following
question/issue:

It happens that U.det() is 1. So I can compute the inverse just by
U.adjugate(). Now I would think that U.inverse_ADJ() leads to the same
result, but it does not. It returns a more complex one.
Reason: it uses .berkowitz_det() which in my case is a quite complex
expression. Only when it is expanded it simplifies to 1.

So I would suggest in /matrices/matrices.py to replace

d = self.berkowitz_det()

by

d = expand(self.berkowitz_det())

and/or to introduce an optional method argument like .det, .adjugate and
related methods already have.

smichr wrote:
> Bastian: regarding the hang in expansion...please try it on branch
> 1766 at my smichr account at github. Ronan timed it recently and it
> looks like it's doing well.
>

Thanks for the suggestion. I'll try that at the weekend


Regarding the more or less the original topic:
Would it be useful to test if .count_ops (or some other (fast) measure
for complexity) exceeds a pre-determined value and in this case print a
warning that the operation (e.g. expand) could take very long? I
sometimes feel like missing this. How about other users?


Regards,
Bastian.

[Vinzent Steinberg]

On 5 mai, 23:43, Bastian Weber <bastian.we...@gmx-topmail.de> wrote:
> Regarding the more or less the original topic:
> Would it be useful to test if .count_ops (or some other (fast) measure
> for complexity) exceeds a pre-determined value and in this case print a
> warning that the operation (e.g. expand) could take very long? I
> sometimes feel like missing this. How about other users?

Such a warning would probably be ok for the interactive usage, but not
for the library. The user is responsible if he does something that
takes a lot of time. Maybe we could add a warning to expand()'s
docstring.

Vinzent