[Numpy-discussion] problems with numdifftools

[Nicolai Heitz]

m 27.10.2010 02:02, schrieb Sebastian Walter:

>  On Wed, Oct 27, 2010 at 12:59 AM, Pauli Virtanen<pav@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).
I was pretty sure that I put it there. Unfortunately it had a different
name there: Errors in calculation of the Hessian using numdifftools. But
I can't find the post by myself at the moment so maybe something went
wrong.
>>  [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.

Probably you are right. I converted it to natural units and it worked
out in the 2 ion case. Increasing the number of ions leads to problems
again and I have no idea where those problems come from.

>>>  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.

Is there by chance any possibility to make PYADOLC run on a (lame)
windows xp engine. If not what else would u recommend (besides switching
to Linux, what I am going to do soon).
>  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.
>
Your code example looks awesome and leads to the correct results. Thank
you very much. I try to make it work on my pc as well.
>  regards,
>  Sebastian
>
>>  --
>>  Pauli Virtanen
>>
>>  _______________________________________________
>>  NumPy-Discussion mailing list
>>  NumPy-Discussion@scipy.org
>>  http://mail.scipy.org/mailman/listinfo/numpy-discussion
>>
>  _______________________________________________
>  NumPy-Discussion mailing list
>  NumPy-Discussion@scipy.org
>  http://mail.scipy.org/mailman/listinfo/numpy-discussion