Oct. 28, 2010
3:50 a.m.
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