The Hessian function that will be passed to the optimizer,
        if needed.
179950
def hessfunc(self, p):
        
        self._set_stochastics(p)
        for i in xrange(self.len):

            di = self.diff(i)
            self.hess[i, i] = self.diff(i, 2)

            if i < self.len - 1:

                for j in xrange(i + 1, self.len):
                    dij = self.diff2(i, j)

                    self.hess[i, j] = dij
                    self.hess[j, i] = dij

        return -1. * self.hess
959730
def hessian(self, x, y, kwargs, diff=diff):
        
        alpha_ra, alpha_dec = self.alpha(x, y, kwargs)

        alpha_ra_dx, alpha_dec_dx = self.alpha(x + diff, y, kwargs)
        alpha_ra_dy, alpha_dec_dy = self.alpha(x, y + diff, kwargs)

        dalpha_rara = (alpha_ra_dx - alpha_ra)/diff
        dalpha_radec = (alpha_ra_dy - alpha_ra)/diff
        dalpha_decra = (alpha_dec_dx - alpha_dec)/diff
        dalpha_decdec = (alpha_dec_dy - alpha_dec)/diff

        f_xx = dalpha_rara
        f_yy = dalpha_decdec
        f_xy = dalpha_radec
        f_yx = dalpha_decra
        return f_xx, f_xy, f_yx, f_yy
791235
def hessian(self, x, y, amp, sigma, e1, e2, center_x=0, center_y=0):
        
        alpha_ra, alpha_dec = self.derivatives(x, y, amp, sigma, e1, e2, center_x, center_y)
        diff = self._diff
        alpha_ra_dx, alpha_dec_dx = self.derivatives(x + diff, y, amp, sigma, e1, e2, center_x, center_y)
        alpha_ra_dy, alpha_dec_dy = self.derivatives(x, y + diff, amp, sigma, e1, e2, center_x, center_y)

        f_xx = (alpha_ra_dx - alpha_ra) / diff
        f_xy = (alpha_ra_dy - alpha_ra) / diff
        # f_yx = (alpha_dec_dx - alpha_dec)/diff
        f_yy = (alpha_dec_dy - alpha_dec) / diff
        return f_xx, f_yy, f_xy
424955
def optimize(self, angles0, target):
        """"""
        def new_objective(angles):
            return self.f(angles, target)

        return scipy.optimize.minimize(
            new_objective,
            angles0,
            **self.optimizer_opt).x
571153
def _minimize_and_clip(optimizer, objective, var_list, clip_val=10):
    
    gradients = optimizer.compute_gradients(objective, var_list=var_list)
    for i, (grad, var) in enumerate(gradients):
        if grad is not None:
            gradients[i] = (tf.clip_by_norm(grad, clip_val), var)
    return gradients
585197
def _get_optimizer(self):
        
        optim = tf.train.AdagradOptimizer(self.learning_rate)
        gradients = optim.compute_gradients(self.cost)
        if self.log_dir:
            for name, (g, v) in zip(['W', 'C', 'bw', 'bc'], gradients):
                tf.summary.histogram("{}_grad".format(name), g)
                tf.summary.histogram("{}_vals".format(name), v)
        return optim.apply_gradients(gradients)
681440
def precompute_optimzation_Y(laplacian_matrix, n_samples, relaxation_kwds):
    """"""
    relaxation_kwds.setdefault('presave',False)
    relaxation_kwds.setdefault('presave_name','pre_comp_current.npy')
    relaxation_kwds.setdefault('verbose',False)
    if relaxation_kwds['verbose']:
        print ('Making Lk and nbhds')
    Lk_tensor, nbk, si_map = \
        compute_Lk(laplacian_matrix, n_samples, relaxation_kwds['subset'])
    if relaxation_kwds['presave']:
        raise NotImplementedError('Not yet implemented presave')
    return { 'Lk': Lk_tensor, 'nbk': nbk, 'si_map': si_map }
1080511
def solve(self, angles0, target):
        """"""
        return self.optimizer.optimize(np.array(angles0), target)
1046776
def hessian(self, x, y, theta_E, center_x=0, center_y=0):
        
        x_shift = x - center_x
        y_shift = y - center_y
        R = (x_shift*x_shift + y_shift*y_shift)**(3./2)
        if isinstance(R, int) or isinstance(R, float):
            prefac = theta_E / max(0.000001, R)
        else:
            prefac = np.empty_like(R)
            r = R[R>0]  #in the SIS regime
            prefac[R==0] = 0.
            prefac[R>0] = theta_E / r

        f_xx = y_shift*y_shift * prefac
        f_yy = x_shift*x_shift * prefac
        f_xy = -x_shift*y_shift * prefac
        return f_xx, f_yy, f_xy
717756
def hessian(self, x, y, n_sersic, R_sersic, k_eff, center_x=0, center_y=0):
        
        x_ = x - center_x
        y_ = y - center_y
        r = np.sqrt(x_**2 + y_**2)
        if isinstance(r, int) or isinstance(r, float):
            r = max(self._s, r)
        else:
            r[r < self._s] = self._s
        d_alpha_dr = self.d_alpha_dr(x, y, n_sersic, R_sersic, k_eff, center_x, center_y)
        alpha = -self.alpha_abs(x, y, n_sersic, R_sersic, k_eff, center_x, center_y)

        #f_xx_ = d_alpha_dr * calc_util.d_r_dx(x_, y_) * x_/r + alpha * calc_util.d_x_diffr_dx(x_, y_)
        #f_yy_ = d_alpha_dr * calc_util.d_r_dy(x_, y_) * y_/r + alpha * calc_util.d_y_diffr_dy(x_, y_)
        #f_xy_ = d_alpha_dr * calc_util.d_r_dy(x_, y_) * x_/r + alpha * calc_util.d_x_diffr_dy(x_, y_)

        f_xx = -(d_alpha_dr/r + alpha/r**2) * x_**2/r + alpha/r
        f_yy = -(d_alpha_dr/r + alpha/r**2) * y_**2/r + alpha/r
        f_xy = -(d_alpha_dr/r + alpha/r**2) * x_*y_/r

        return f_xx, f_yy, f_xy
38561
def gradient(self):
        """"""
        functional = self
        c = self.scale

        class RosenbrockGradient(Operator):

            """The gradient operator of the Rosenbrock functional."""

            def __init__(self):
                """Initialize a new instance."""
                super(RosenbrockGradient, self).__init__(
                    functional.domain, functional.domain, linear=False)

            def _call(self, x, out):
                """Apply the gradient operator to the given point."""
                for i in range(1, self.domain.size - 1):
                    out[i] = (2 * c * (x[i] - x[i - 1]**2) -
                              4 * c * (x[i + 1] - x[i]**2) * x[i] -
                              2 * (1 - x[i]))
                out[0] = (-4 * c * (x[1] - x[0] ** 2) * x[0] +
                          2 * (x[0] - 1))
                out[-1] = 2 * c * (x[-1] - x[-2] ** 2)

            def derivative(self, x):
                """The derivative of the gradient.

                This is also known as the Hessian.
                """
                # TODO: Implement optimized version of this that does not need
                # a matrix.
                shape = (functional.domain.size, functional.domain.size)
                matrix = np.zeros(shape)

                # Straightforward computation
                for i in range(0, self.domain.size - 1):
                    matrix[i, i] = (2 * c + 2 + 12 * c * x[i] ** 2 -
                                    4 * c * x[i + 1])
                    matrix[i + 1, i] = -4 * c * x[i]
                    matrix[i, i + 1] = -4 * c * x[i]
                matrix[-1, -1] = 2 * c
                matrix[0, 0] = 2 + 12 * c * x[0] ** 2 - 4 * c * x[1]
                return MatrixOperator(matrix, self.domain, self.range)

        return RosenbrockGradient()
997844
def hessian(self):
        
        return [[[sympy.diff(partial_dv, param) for param in self.params]
                 for partial_dv in comp] for comp in self.jacobian]
421358
def opt(self, x_init, f_fp=None, f=None, fp=None):
        
        rcstrings = ['','Maximum number of iterations exceeded', 'Gradient and/or function calls not changing']

        opt_dict = {}
        if self.xtol is not None:
            print("WARNING: bfgs doesn't have an xtol arg, so I'm going to ignore it")
        if self.ftol is not None:
            print("WARNING: bfgs doesn't have an ftol arg, so I'm going to ignore it")
        if self.gtol is not None:
            opt_dict['gtol'] = self.gtol

        opt_result = optimize.fmin_bfgs(f, x_init, fp, disp=self.messages,
                                            maxiter=self.max_iters, full_output=True, **opt_dict)
        self.x_opt = opt_result[0]
        self.f_opt = f_fp(self.x_opt)[0]
        self.funct_eval = opt_result[4]
        self.status = rcstrings[opt_result[6]]
836155
def eval_hessian(self, *args, **kwargs):
        
        # Evaluate the hessian model and use the resulting Ans namedtuple as a
        # dict. From this, take the relevant components.
        eval_hess_dict = self.hessian_model(*args, **kwargs)._asdict()
        hess = [[[np.broadcast_to(eval_hess_dict.get(D(var, p1, p2), 0),
                                  eval_hess_dict[var].shape)
                    for p2 in self.params]
                for p1 in self.params]
            for var in self
        ]
        # Use numpy to broadcast these arrays together and then stack them along
        # the parameter dimension. We do not include the component direction in
        # this, because the components can have independent shapes.
        for idx, comp in enumerate(hess):
            hess[idx] = np.stack(np.broadcast_arrays(*comp))

        Ans = variabletuple('Ans', self.keys())
        return Ans(*hess)
434713
def optimize(self, x0, target):
        """"""
        x = x0
        for i in range(self.maxiter):
            g = self.g(x, target)
            h = self.h(x, target)
            if i == 0:
                alpha = 0
                m = g
            else:
                alpha = - np.dot(m, np.dot(h, g)) / np.dot(m, np.dot(h, m))
                m = g + np.dot(alpha, m)
            t = - np.dot(m, g) / np.dot(m, np.dot(h, m))
            delta = np.dot(t, m)
            x = x + delta
            if np.linalg.norm(delta) < self.tol:
                break
        return x
899703
public boolean optimize(ResidualFunction function, DMatrixRMaj parameters )
    {
        configure(function,parameters.getNumElements());

        // save the cost of the initial parameters so that it knows if it improves or not
        double previousCost = initialCost = cost(parameters);

        // iterate until the difference between the costs is insignificant
        double lambda = initialLambda;

        // if it should recompute the Jacobian in this iteration or not
        boolean computeHessian = true;

        for( int iter = 0; iter < maxIterations; iter++ ) {
            if( computeHessian ) {
                // compute some variables based on the gradient
                computeGradientAndHessian(parameters);
                computeHessian = false;

                // check for convergence using gradient test
                boolean converged = true;
                for (int i = 0; i < g.getNumElements(); i++) {
                    if( Math.abs(g.data[i]) > gtol ) {
                        converged = false;
                        break;
                    }
                }
                if( converged )
                    return true;
            }

            // H = H + lambda*I
            for (int i = 0; i < H.numRows; i++) {
                H.set(i,i, Hdiag.get(i) + lambda);
            }

            // In robust implementations failure to solve is handled much better
            if( !CommonOps_DDRM.solve(H, g, negativeStep) ) {
                return false;
            }

            // compute the candidate parameters
            CommonOps_DDRM.subtract(parameters, negativeStep, candidateParameters);

            double cost = cost(candidateParameters);
            if( cost <= previousCost ) {
                // the candidate parameters produced better results so use it
                computeHessian = true;
                parameters.set(candidateParameters);

                // check for convergence
                // ftol <= (cost(k) - cost(k+1))/cost(k)
                boolean converged = ftol*previousCost >= previousCost-cost;

                previousCost = cost;
                lambda /= 10.0;

                if( converged ) {
                    return true;
                }
            } else {
                lambda *= 10.0;
            }

        }
        finalCost = previousCost;
        return true;
    }
707887
def optimize
      # Ensure pending variables and constraints are added
      update

      # Run the solver and save the status for later
      iocp = GLPK::IOCP.new
      GLPK.glp_init_iocp iocp
      iocp[:presolve] = GLPK::GLP_ON
      status = GLPK.glp_intopt(@ptr, iocp)

      save_solution status
    end