[Scipy-tickets] [SciPy] #1248: Error in scipy.optimize.rosen_hess(x) , in H[0, 0] position, when x[0] != 1

SciPy Trac scipy-tickets@scipy....
Mon Jul 26 08:19:59 CDT 2010


#1248: Error in scipy.optimize.rosen_hess(x) ,   in H[0,0] position, when x[0] !=
1
-------------------------+--------------------------------------------------
 Reporter:  johngabriel  |       Owner:  somebody
     Type:  defect       |      Status:  new     
 Priority:  normal       |   Milestone:  0.9.0   
Component:  Other        |     Version:  0.7.0   
 Keywords:               |  
-------------------------+--------------------------------------------------
 #Error in scipy.optimize.rosen_hess,  in H[0,0] position, when x[0] != 1

 #Error is illustrated below, where:-
 #   "rosen_hess(x) dot product with p,
 #    does not give same result as rosen_hess_prod(x, p)"

 """
 PythonWin 2.6.2 (r262:71605, Apr 14 2009, 22:40:02)
 [MSC v.1500 32 bit (Intel)] on win32.
 >>> import scipy as sp
 >>> sp.version
 <module 'scipy.version' from
 'C:\Python26\lib\site-packages\scipy\version.pyc'>
 >>> sp.version.version
 '0.7.2'
 >>> from scipy.optimize import rosen_hess, rosen_hess_prod
 >>> import numpy as np
 >>> x = np.array([3, 4, 5]) # test vector
 >>> p = np.array([2, 2, 2]) # weights

 # rosen_hess_prod is correct
 >>> rosen_hess_prod(x, p)
 array([16004, 29204, -2800])

 # rosen_hess dot product with p, should be the same as above,
 # but is NOT correct in [0]th position.
 >>> np.dot(rosen_hess(x), p)
 array([ 1604, 29204, -2800])

 # Size of error in [0]th position, due to rosen_hess.
 >>> 16004 - 1604
 14400

 # Explains size of error is in H[0,0] term of rosen_hess,
 # where there is a missing square of x[0] term in line 4 of code.
 >>> 1200 * (x[0]**2 - x[0]) * p[1]
 14400
 >>>

 """
 # Note that the error shown above would not come to light
 # if rosen_hess is tested with the unit vector,
 # or indeed any vector with 1.0 in the [0]th position.

 #A Corrected version of rosen_hess follows:-

 def rosen_hess(x):
     # ****** THIS IS NOW CORRECTED. H[0,0] is now correct.
     # See corrected 4th line of code*****
     # Returns a symetric ndarray of all of the second partial derivatives,
     # i.e. the Hessian matrix, H, of the Rosenbrock function,
     # w.r.t the unknowns x[0], x[1], x[2] ...
     # and again w.r.t x[0], x[1], x[2] ...
     x = np.atleast_1d(x)
     H = np.diag(-400*x[:-1],1) - np.diag(400*x[:-1],-1)
     diagonal = np.zeros(len(x), dtype=x.dtype)

     # Original code , with 1200*x[0]   ,i.e. missing "square" of x[0]
     #  diagonal[0] = 1200*x[0]-400*x[1]+2

     # Corrected code, with 1200*x[0]**2, J.R.Gabriel, 25th July, 2010
     diagonal[0] = 1200*x[0]**2-400*x[1]+2

     diagonal[-1] = 200
     diagonal[1:-1] = 202 + 1200*x[1:-1]**2 - 400*x[2:]
     H = H + np.diag(diagonal)
     return H

 # Test Section is only run if this is the main program
 if __name__ == '__main__':
     from scipy.optimize import rosen_hess_prod
     import numpy as np
     x = np.array([3, 4, 5]) # test vector
     p = np.array([2, 2, 2]) # weights

     # rosen_hess_prod is correct
     hp1 = rosen_hess_prod(x, p)     # array([16004, 29204, -2800])

     # corrected rosen_hess dot product with p, now gives the same as
 above,
     # in [0]th position.
     hp2 = np.dot(rosen_hess(x), p)  # array([16004, 29204, -2800])

     assert np.allclose(hp1, hp2)

     print "hp1 =", hp1
     print "hp2 =", hp2

-- 
Ticket URL: <http://projects.scipy.org/scipy/ticket/1248>
SciPy <http://www.scipy.org>
SciPy is open-source software for mathematics, science, and engineering.


More information about the Scipy-tickets mailing list