[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