[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....
Tue Oct 12 01:58:35 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: scipy.optimize | Version: 0.7.0
Keywords: |
----------------------------+-----------------------------------------------
Description changed by warren.weckesser:
Old description:
> #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
New description:
{{{
#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#comment:2>
SciPy <http://www.scipy.org>
SciPy is open-source software for mathematics, science, and engineering.
More information about the Scipy-tickets
mailing list