[Numpy-discussion] problems with numdifftools
Sebastian Walter
sebastian.walter@gmail....
Thu Oct 28 14:34:58 CDT 2010
hmm, I have just realized that I forgot to upload the new version to pypi:
it is now available on
http://pypi.python.org/pypi/algopy
On Thu, Oct 28, 2010 at 10:47 AM, Sebastian Walter
<sebastian.walter@gmail.com> wrote:
> 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