MAINT: Some improvements to optimize._numdiff

scipy/optimize/_numdiff.py

@@ -177,14 +177,14 @@ def group_columns(A, order=0):
 
 
 def approx_derivative(fun, x0, method='3-point', rel_step=None, f0=None,
-                      args=(), kwargs={}, bounds=(-np.inf, np.inf),
-                      sparsity=None):
+                      bounds=(-np.inf, np.inf), sparsity=None, args=(),
+                      kwargs={}):
     """Compute finite difference approximation of the derivatives of a
     vector-valued function.
 
     If a function maps from R^n to R^m, its derivatives form m-by-n matrix
-    called Jacobian, where an element (i, j) is a partial derivative of f[i]
-    with respect to x[j].
+    called the Jacobian, where an element (i, j) is a partial derivative of
+    f[i] with respect to x[j].
 
     Parameters
     ----------
@@ -215,36 +215,42 @@ def approx_derivative(fun, x0, method='3-point', rel_step=None, f0=None,
     f0 : None or array_like, optional
         If not None it is assumed to be equal to ``fun(x0)``, in  this case
         the ``fun(x0)`` is not called. Default is None.
-    args, kwargs : tuple and dict, optional
-        Additional arguments passed to `fun`. Both empty by default.
-        The calling signature is ``fun(x, *args, **kwargs)``.
     bounds : tuple of array_like, optional
         Lower and upper bounds on independent variables. Defaults to no bounds.
         Each bound must match the size of `x0` or be a scalar, in the latter
         case the bound will be the same for all variables. Use it to limit the
         range of function evaluation.
-    sparsity : None or (structure, groups) tuple, optional
-        Defines sparsity structure of Jacobian.
+    sparsity : {None, array_like, sparse matrix, 2-tuple}, optional
+        Defines a sparsity structure of the Jacobian matrix. If the Jacobian
+        matrix is known to have only few non-zero elements in each row, then
+        it's possible to estimate its several columns by a single function
+        evaluation [3]_. To perform such economic computations two ingredients
+        are required:
 
         * structure : array_like or sparse matrix of shape (m, n). A zero
-          element means that a corresponding element of Jacobian identically
-          equals to zero.
-        * groups : array_like of shape (n,). Precomputed columns grouping
-          for a given sparsity structure, function `group_columns` should be
-          used beforehand to compute it.
-
-        Note, that sparse differencing makes sense only for actually large
-        and sparse Jacobians. If None (default) standard dense differencing
-        will be used.
+          element means that a corresponding element of the Jacobian
+          identically equals to zero.
+        * groups : array_like of shape (n,). A column grouping for a given
+          sparsity structure, use `group_columns` to obtain it.
+
+        A single array or a sparse matrix is interpreted as a sparsity
+        structure, and groups are computed inside the function. A tuple is
+        interpreted as (structure, groups). If None (default), a standard
+        dense differencing will be used.
+
+        Note, that sparse differencing makes sense only for large Jacobian
+        matrices where each row contains few non-zero elements.
+    args, kwargs : tuple and dict, optional
+        Additional arguments passed to `fun`. Both empty by default.
+        The calling signature is ``fun(x, *args, **kwargs)``.
 
     Returns
     -------
     J : ndarray or csr_matrix
-        Finite difference approximation of derivatives in the form of Jacobian
-        matrix. If `sparse` is None then ndarray with shape (m, n) is
-        returned. Although if m=1 it is returned as a gradient with
-        shape (n,). If `sparse` is not None, csr_matrix with shape (m, n) is
-        returned.
+        Finite difference approximation of the Jacobian matrix. If `sparsity`
+        is None then ndarray with shape (m, n) is returned. Although if m=1 it
+        is returned as a gradient with shape (n,). If `sparsity` is not None,
+        csr_matrix with shape (m, n) is returned.
 
     See Also
     --------
@@ -276,7 +282,11 @@ def approx_derivative(fun, x0, method='3-point', rel_step=None, f0=None,
     .. [1] W. H. Press et. al. "Numerical Recipes. The Art of Scientific
            Computing. 3rd edition", sec. 5.7.
 
-    .. [2] B. Fornberg, "Generation of Finite Difference Formulas on
+    .. [2] A. Curtis, M. J. D. Powell, and J. Reid, "On the estimation of
+           sparse Jacobian matrices", Journal of the Institute of Mathematics
+           and its Applications, 13 (1974), pp. 117-120.
+
+    .. [3] B. Fornberg, "Generation of Finite Difference Formulas on
            Arbitrarily Spaced Grids", Mathematics of Computation 51, 1988.
 
     Examples
@@ -348,11 +358,17 @@ def fun_wrapped(x):
     if sparsity is None:
         return _dense_difference(fun_wrapped, x0, f0, h, use_one_sided, method)
     else:
-        structure, groups = sparsity
+        if not issparse(sparsity) and len(sparsity) == 2:
+            structure, groups = sparsity
+        else:
+            structure = sparsity
+            groups = group_columns(sparsity)
+
         if issparse(structure):
             structure = csc_matrix(structure)
         else:
             structure = np.atleast_2d(structure)
+
         groups = np.atleast_1d(groups)
         return _sparse_difference(fun_wrapped, x0, f0, h, use_one_sided,
                                   structure, groups, method)
@@ -478,11 +494,10 @@ def _sparse_difference(fun, x0, f0, h, use_one_sided,
     return csr_matrix(J)
 
 
-def check_derivative(fun, jac, x0, args=(), kwargs={},
-                     bounds=(-np.inf, np.inf), sparse_diff=False,
-                     sparsity=None):
+def check_derivative(fun, jac, x0, bounds=(-np.inf, np.inf), args=(),
+                     kwargs={}):
     """Check correctness of a function computing derivatives (Jacobian or
-    gradient) by comparison with finite difference approximation.
+    gradient) by comparison with a finite difference approximation.
 
     Parameters
     ----------
@@ -493,35 +508,19 @@ def check_derivative(fun, jac, x0, args=(), kwargs={},
     jac : callable
         Function which computes Jacobian matrix of `fun`. It must work with
         argument x the same way as `fun`. The return value must be array_like
-        with an appropriate shape.
+        or sparse matrix with an appropriate shape.
     x0 : array_like of shape (n,) or float
         Point at which to estimate the derivatives. Float will be converted
         to 1-d array.
-    args, kwargs : tuple and dict, optional
-        Additional arguments passed to `fun` and `jac`. Both empty by default.
-        The calling signature is ``fun(x, *args, **kwargs)`` and the same
-        for `jac`.
     bounds : 2-tuple of array_like, optional
         Lower and upper bounds on independent variables. Defaults to no bounds.
         Each bound must match the size of `x0` or be a scalar, in the latter
         case the bound will be the same for all variables. Use it to limit the
         range of function evaluation.
-    sparse_diff : boolean, optional
-        If `jac` returns a sparse matrix this flag indicates whether to
-        use sparse differencing for Jacobian estimation.
-    sparsity : None or (structure, groups) tuple, optional
-        If `jac` returns a sparse matrix and ``sparse_diff=True``, provide a
-        tuple with following elements:
-
-        * structure : array_like or sparse matrix of shape (m, n). A zero
-          element means that a corresponding element of Jacobian identically
-          equals to zero.
-        * groups : array_like of shape (n,). Precomputed columns grouping
-          for a given sparsity structure, function `group_columns` should be
-          used beforehand to compute it.
-
-        If None `structure` will be taken from evaluated `jac` and `groups`
-        will be computed.
+    args, kwargs : tuple and dict, optional
+        Additional arguments passed to `fun` and `jac`. Both empty by default.
+        The calling signature is ``fun(x, *args, **kwargs)`` and the same
+        for `jac`.
 
     Returns
     -------
@@ -557,22 +556,17 @@ def check_derivative(fun, jac, x0, args=(), kwargs={},
     2.4492935982947064e-16
     """
     J_to_test = jac(x0, *args, **kwargs)
-    if issparse(J_to_test) and sparse_diff:
-        if sparsity is None:
-            groups = group_columns(J_to_test)
-            sparsity = (J_to_test, groups)
-        J_diff = approx_derivative(fun, x0, args=args, kwargs=kwargs,
-                                   bounds=bounds, sparsity=sparsity)
+    if issparse(J_to_test):
+        J_diff = approx_derivative(fun, x0, bounds=bounds, sparsity=J_to_test,
+                                   args=args, kwargs=kwargs)
         J_to_test = csr_matrix(J_to_test)
         abs_err = J_to_test - J_diff
         i, j, abs_err_data = find(abs_err)
-        jac_diff_data = np.asarray(J_diff[i, j]).ravel()
+        J_diff_data = np.asarray(J_diff[i, j]).ravel()
         return np.max(np.abs(abs_err_data) /
-                      np.maximum(1, np.abs(jac_diff_data)))
+                      np.maximum(1, np.abs(J_diff_data)))
     else:
-        J_diff = approx_derivative(
-            fun, x0, args=args, kwargs=kwargs, bounds=bounds)
-        if issparse(J_to_test):
-            J_to_test = J_to_test.toarray()
+        J_diff = approx_derivative(fun, x0, bounds=bounds,
+                                   args=args, kwargs=kwargs)
         abs_err = np.abs(J_to_test - J_diff)
         return np.max(abs_err / np.maximum(1, np.abs(J_diff)))

scipy/optimize/tests/test__numdiff.py

@@ -443,6 +443,11 @@ def test_all(self):
                                   rel_step=rel_step, sparsity=(A, groups))
             assert_allclose(J.toarray(), self.J_true, rtol=1e-5)
 
+    def test_no_precomputed_groups(self):
+        A = self.structure(self.n)
+        J = approx_derivative(self.fun, self.x0, sparsity=A)
+        assert_allclose(J.toarray(), self.J_true, rtol=1e-6)
+
     def test_equivalence(self):
         structure = np.ones((self.n, self.n), dtype=int)
         groups = np.arange(self.n)
@@ -456,28 +461,14 @@ def test_check_derivative(self):
         def jac(x):
             return csr_matrix(self.jac(x))
 
-        accuracy = check_derivative(
-            self.fun, jac, self.x0,
-            bounds=(self.lb, self.ub), sparse_diff=True)
+        accuracy = check_derivative(self.fun, jac, self.x0,
+                                    bounds=(self.lb, self.ub))
         assert_(accuracy < 1e-9)
 
-        A = self.structure(self.n)
-        groups = group_columns(A)
-        accuracy = check_derivative(
-            self.fun, jac, self.x0, bounds=(self.lb, self.ub),
-            sparse_diff=True, sparsity=(A, groups)
-        )
+        accuracy = check_derivative(self.fun, jac, self.x0,
+                                    bounds=(self.lb, self.ub))
         assert_(accuracy < 1e-9)
 
-        accuracy = check_derivative(
-            self.fun, jac, self.x0,
-            bounds=(self.lb, self.ub), sparse_diff=False)
-        # Slightly worse accuracy because all elements are computed.
-        # Floating point issues make true 0 to some smal value, then
-        # it is divided by small step and as a result we have ~1e-9 element,
-        # which is actually should be zero.
-        assert_(accuracy < 1e-8)
-
 
 if __name__ == '__main__':
     run_module_suite()