[Numpy-discussion] problems with numdifftools

Sebastian Walter sebastian.walter@gmail....
Thu Oct 28 03:47:28 CDT 2010


On Wed, Oct 27, 2010 at 10:50 PM, Nicolai Heitz <nicolaiheitz@gmx.de> wrote:
> 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).

1) PYADOLC
A windows version has been requested several times now. But until
recently ADOL-C wasn't available as windows version.
So yes, in principle it should be possible to get it to work on windows:

You will need
1) boost:python
http://www.boost.org/doc/libs/1_44_0/libs/python/doc/index.html
2) ADOL-C sources http://www.coin-or.org/projects/ADOL-C.xml
3) scons http://www.scons.org/
on windows.

 If you want to give it a try I could help to get it to work.

2) Alternatives:
You can also try the ALGOPY which is pure Python and is known to work
on Linux and Windows. The installation is also very easy (setup.py
build or setup.py install)
I have added your problem to the ALGOPY documentation:
http://packages.python.org/algopy/examples/hessian_of_potential_function.html
The catch is that ALGOPY is not as mature as PYADOLC. However, if you
are careful to write clean code it should work reliably.


Sebastian

>>  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
>
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>


More information about the NumPy-Discussion mailing list