# [Scipy-svn] r6833 - in trunk/scipy/optimize: . tests

scipy-svn@scip... scipy-svn@scip...
Tue Oct 12 15:02:13 CDT 2010

```Author: warren.weckesser
Date: 2010-10-12 15:02:12 -0500 (Tue, 12 Oct 2010)
New Revision: 6833

Modified:
trunk/scipy/optimize/optimize.py
trunk/scipy/optimize/tests/test_optimize.py
Log:
BUG: optimize: Fixed formula in rosen_hess (ticket #1248)--thanks to John Gabriel.  Also added docstrings to the rosen* functions.

Modified: trunk/scipy/optimize/optimize.py
===================================================================
--- trunk/scipy/optimize/optimize.py	2010-10-10 19:42:10 UTC (rev 6832)
+++ trunk/scipy/optimize/optimize.py	2010-10-12 20:02:12 UTC (rev 6833)
@@ -23,7 +23,7 @@
__docformat__ = "restructuredtext en"

import numpy
-from numpy import atleast_1d, eye, mgrid, argmin, zeros, shape, empty, \
+from numpy import atleast_1d, eye, mgrid, argmin, zeros, shape, \
squeeze, vectorize, asarray, absolute, sqrt, Inf, asfarray, isinf
from linesearch import \
line_search_BFGS, line_search_wolfe1, line_search_wolfe2, \
@@ -64,11 +64,47 @@
else:
return numpy.sum(abs(x)**ord,axis=0)**(1.0/ord)

-def rosen(x):  # The Rosenbrock function
+def rosen(x):
+    """The Rosenbrock function.
+
+    The function computed is
+
+        sum(100.0*(x[1:] - x[:-1]**2.0)**2.0 + (1 - x[:-1])**2.0
+
+    Parameters
+    ----------
+    x : array_like, 1D
+        The point at which the Rosenbrock function is to be computed.
+
+    Returns
+    -------
+    f : float
+        The value of the Rosenbrock function
+
+    --------
+    rosen_der, rosen_hess, rosen_hess_prod
+    """
x = asarray(x)
return numpy.sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0,axis=0)

def rosen_der(x):
+    """The derivative (i.e. gradient) of the Rosenbrock function.
+
+    Parameters
+    ----------
+    x : array_like, 1D
+        The point at which the derivative is to be computed.
+
+    Returns
+    -------
+    der : 1D numpy array
+        The gradient of the Rosenbrock function at `x`.
+
+    --------
+    rosen, rosen_hess, rosen_hess_prod
+    """
x = asarray(x)
xm = x[1:-1]
xm_m1 = x[:-2]
@@ -80,16 +116,51 @@
return der

def rosen_hess(x):
+    """The Hessian matrix of the Rosenbrock function.
+
+    Parameters
+    ----------
+    x : array_like, 1D
+        The point at which the Hessian matrix is to be computed.
+
+    Returns
+    -------
+    hess : 2D numpy array
+        The Hessian matrix of the Rosenbrock function at `x`.
+
+    --------
+    rosen, rosen_der, rosen_hess_prod
+    """
x = atleast_1d(x)
H = numpy.diag(-400*x[:-1],1) - numpy.diag(400*x[:-1],-1)
diagonal = numpy.zeros(len(x), dtype=x.dtype)
-    diagonal[0] = 1200*x[0]-400*x[1]+2
+    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 + numpy.diag(diagonal)
return H

def rosen_hess_prod(x,p):
+    """Product of the Hessian matrix of the Rosenbrock function with a vector.
+
+    Parameters
+    ----------
+    x : array_like, 1D
+        The point at which the Hessian matrix is to be computed.
+    p : array_like, 1D, same size as `x`.
+        The vector to be multiplied by the Hessian matrix.
+
+    Returns
+    -------
+    v : 1D numpy array
+        The Hessian matrix of the Rosenbrock function at `x` multiplied
+        by the vector `p`.
+
+    --------
+    rosen, rosen_der, rosen_hess
+    """
x = atleast_1d(x)
Hp = numpy.zeros(len(x), dtype=x.dtype)
Hp[0] = (1200*x[0]**2 - 400*x[1] + 2)*p[0] - 400*x[0]*p[1]

Modified: trunk/scipy/optimize/tests/test_optimize.py
===================================================================
--- trunk/scipy/optimize/tests/test_optimize.py	2010-10-10 19:42:10 UTC (rev 6832)
+++ trunk/scipy/optimize/tests/test_optimize.py	2010-10-12 20:02:12 UTC (rev 6833)
@@ -11,7 +11,7 @@
"""

from numpy.testing import assert_raises, assert_almost_equal, \
-        assert_, TestCase, run_module_suite
+        assert_equal, assert_, TestCase, run_module_suite

from scipy import optimize
from numpy import array, zeros, float64, dot, log, exp, inf, sin, cos
@@ -345,5 +345,17 @@
if ef > 1e-8:
raise err

+
+class TestRosen(TestCase):
+
+    def test_hess(self):
+        """Compare rosen_hess(x) times p with rosen_hess_prod(x,p) (ticket #1248)"""
+        x = array([3, 4, 5])
+        p = array([2, 2, 2])
+        hp = optimize.rosen_hess_prod(x, p)
+        dothp = np.dot(optimize.rosen_hess(x), p)
+        assert_equal(hp, dothp)
+
+
if __name__ == "__main__":
run_module_suite()

```