[Numpy-discussion] problems with numdifftools

[Sebastian Walter]

On Wed, Oct 27, 2010 at 12:59 AM, Pauli Virtanen <pav at iki.fi> wrote:
> Tue, 26 Oct 2010 14:24:39 -0700, Nicolai Heitz wrote:
>> >  http://mail.scipy.org/mailman/listinfo/scipy-user
>>
>> I contacted them already but they didn't responded so far and I was
>> forwarded to that list which was supposed to be more appropriated.
>
> I think you are thinking here about some other list -- scipy-user
> is the correct place for this discussion (and I don't remember seeing
> your mail there).
>
> [clip]
>> 1) Can I make it run/fix it, so that it is also going to work for the SI
>> scaling?
>
> Based on a brief look, it seems that uniform scaling will not help you,
> as you have two very different length scales in the problem,
>
>        1/sqrt(m w^2) >> C
>
> If you go to CM+relative coordinates you might be able to scale them
> separately, but that's fiddly and might not work for larger N.
>
> In your problem, things go wrong when the ratio between the
> length scales approaches 1e-15 which happens to be the machine epsilon.
> This implies that the algorithm runs into some problems caused by the
> finite precision of floating-point numbers.
>
> What exactly goes wrong and how to fix it, no idea --- I didn't look into
> how Numdifftools is implemented.
>
>> 2) How can I be sure that increasing the number of ions or adding a
>> somehow more complicated term to the potential energy is not causing the
>> same problems even in natural units?
>>
>> 3) In which range is numdifftools working properly.
>
> That depends on the algorithm and the problem. Personally, I wouldn't
> trust numerical differentiation if the problem has significantly
> different length scales, it is important to capture all of them
> accurately, and it is not clear how to scale them to the same size.
> Writing ND software that works as expected all the time is probably
> not possible even in theory.
>
> Numerical differentiation is not the only game in the town. I'd look
> into automatic differentiation (AD) -- there are libraries available
> for Python also for that, and it is numerically stable.
>
> E.g.
>
> http://en.wikipedia.org/wiki/Automatic_differentiation#Software
>
> has a list of Python libraries. I don't know which of them would be
> the best ones, though.
>

they all have their pro's and con's.
Being (co-)author of some of these tools, my personal and very biased advice is:
if you are on Linux, I would go for PYADOLC. it provides bindings to a
feature-rich and well-tested C++ library.
However, the installation is a little tricker than a "setup.py build"
since you will need to compile ADOL-C and get Boost::Python to work.
PYADOLC can also differentiate much more complicated code than your
example in a relatively efficient manner.

For your example the code looks like:

--------------------- code ---------------------------

....

c=classicalHamiltonian()
xopt = optimize.fmin(c.potential, c.initialposition(), xtol = 1e-10)

import adolc;  import numpy

# trace the computation
adolc.trace_on(0)
x = adolc.adouble(c.initialposition())
adolc.independent(x)
y = c.potential(x)
adolc.dependent(y)
adolc.trace_off()

hessian = adolc.hessian(0, xopt)
eigenvalues = numpy.linalg.eigh(hessian)[0]
normal_modes = c.normal_modes(eigenvalues)
print 'hessian=\n',hessian
print 'eigenvalues=\n',eigenvalues
print 'normal_modes=\n',normal_modes
--------------------- code ---------------------------

and you get as output
Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 81
         Function evaluations: 153
hessian=
[[  5.23748399e-12  -2.61873843e-12]
 [ -2.61873843e-12   5.23748399e-12]]
eigenvalues=
[  2.61874556e-12   7.85622242e-12]
normal_modes=
[  6283185.30717959  10882786.30440101]


Also, you should use an eigenvalue solver for symmetric matrices, e.g.
numpy.linalg.eigh.


regards,
Sebastian

> --
> Pauli Virtanen
>
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion at scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>