# [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.
```