title: Understanding the Jacobian output of scipy.optimize.minimize
body: <p>I'm working with scipy.optimize.minimize to find the minimum of the RSS for a custom nonlinear function. I'll provide a simple linear example to illustrate what I am doing:</p>
  
  <pre><code>import numpy as np
  from scipy import optimize
  
  def response(X, b0, b1, b2):
      return b2 * X[1]**2 + b1 * X[0] + b0
  
  def obj_rss(model_params, y_true, X):
      return np.sum((y_true - response(X, *model_params))**2)  
  
  x = np.array([np.arange(0, 10), np.arange(10, 20)])
  r = 15. * x[1]**2 - 32. * x[0] + 10.
  
  init_guess = np.array([0., 50., 10.])
  res = optimize.minimize(obj_rss, init_guess, args=(r, x))
  
  print res
  </code></pre>
  
  <p>This yields the results:</p>
  
  <pre><code>      fun: 3.0218799331864133e-08
   hess_inv: array([[  7.50606278e+00,   2.38939463e+00,  -8.33333575e-02],
         [  2.38939463e+00,   8.02462363e-01,  -2.74621294e-02],
         [ -8.33333575e-02,  -2.74621294e-02,   9.46969972e-04]])
        jac: array([ -3.31359843e-07,  -5.42022462e-08,   2.34304025e-08])
    message: 'Optimization terminated successfully.'
       nfev: 45
        nit: 6
       njev: 9
     status: 0
    success: True
          x: array([ 10.00066577, -31.99978062,  14.99999243])
  </code></pre>
  
  <p>And we see that the fitted parameters 10, -32, and 15 are equivalent to those used to generate the actuals data. That's great. Now my question:</p>
  
  <p>I have the understanding that the Jacobian should be an m x n matrix where m is the number of records from the X input and n is the number of parameters. Clearly I don't have that in the results object. The results object yields an array that is referred to as the Jacobian in the documentation (<a href="http://docs.scipy.org/doc/scipy-0.15.1/reference/generated/scipy.optimize.minimize.html" rel="noreferrer">1</a> and <a href="http://docs.scipy.org/doc/scipy-0.15.1/reference/generated/scipy.optimize.OptimizeResult.html#scipy.optimize.OptimizeResult" rel="noreferrer">2</a>), but is only one-dimensional with a number of elements equal to the number of parameters. </p>
  
  <p>Further confusing the matter, when I use method='SLSQP', the Jacobian that is returned has one more element than that returned by other minimization algorithms. </p>
  
  <p>. . . </p>
  
  <p>My larger goal here is to be able to calculate either confidence intervals or standard errors, t-, and p-values for the fitted parameters, so if you think I'm way off track here, please let me know.</p>
  
  <p>EDIT:</p>
  
  <p>The following is intended to show how the SLSQP minimization algorithm yields different results in the Jacobian than the default minimization algorithm, which is one of BFGS, L-BFGS-B, or SLSQP, depending on if the problem has constraints (as mentioned in the documentation). The SLSQP solver is intended for use with constraints. </p>
  
  <pre><code>import numpy as np
  from scipy import optimize
  
  
  def response(X, b0, b1, b2):
      return b2 * X[1]**2 + b1 * X[0] + b0
  
  def obj_rss(model_params, y_true, X):
      return np.sum((y_true - response(X, *model_params))**2)
  
  
  x = np.array([np.arange(0, 10), np.arange(10, 20)])
  r = 15. * x[1]**2 - 32. * x[0] + 10.
  
  init_guess = np.array([0., 50., 10.])
  res = optimize.minimize(obj_rss, init_guess, method='SLSQP', args=(r, x))
  
  print res
  
  r_pred = response(x, *res.x)
  </code></pre>
  
  <p>Yields results:</p>
  
  <pre><code>     fun: 7.5269461938291697e-10
       jac: array([  2.94677643e-05,   5.52844499e-04,   2.59870917e-02,
           0.00000000e+00])
   message: 'Optimization terminated successfully.'
      nfev: 58
       nit: 10
      njev: 10
    status: 0
   success: True
         x: array([ 10.00004495, -31.9999794 ,  14.99999938])
  </code></pre>
  
  <p>One can see that there is an extra element in the Jacobian array that is returned from the SLSQP solver. I am confused where this comes from. </p>