[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