Just looking for one out of all the possible solutions.
@see "Convex Optimization, C.5 p. 681".
766866
private DoubleMatrix1D findOneRoot(DoubleMatrix2D A, DoubleMatrix1D b) throws Exception{
		return originalProblem.findEqFeasiblePoint(A, b);
	}
743471
static float[] randomPointInConvexPoly(float[] pts, int npts, float[] areas, float s, float t) {
        // Calc triangle araes
        float areasum = 0.0f;
        for (int i = 2; i < npts; i++) {
            areas[i] = triArea2D(pts, 0, (i - 1) * 3, i * 3);
            areasum += Math.max(0.001f, areas[i]);
        }
        // Find sub triangle weighted by area.
        float thr = s * areasum;
        float acc = 0.0f;
        float u = 1.0f;
        int tri = npts - 1;
        for (int i = 2; i < npts; i++) {
            float dacc = areas[i];
            if (thr >= acc && thr < (acc + dacc)) {
                u = (thr - acc) / dacc;
                tri = i;
                break;
            }
            acc += dacc;
        }

        float v = (float) Math.sqrt(t);

        float a = 1 - v;
        float b = (1 - u) * v;
        float c = u * v;
        int pa = 0;
        int pb = (tri - 1) * 3;
        int pc = tri * 3;

        return new float[] { a * pts[pa] + b * pts[pb] + c * pts[pc],
                a * pts[pa + 1] + b * pts[pb + 1] + c * pts[pc + 1],
                a * pts[pa + 2] + b * pts[pb + 2] + c * pts[pc + 2] };
    }
1135418
public static Optimizer<DifferentiableFunction> getRegularizedOptimizer(final Optimizer<DifferentiableFunction> opt,
            final double l1Lambda, final double l2Lambda) {
        if (l1Lambda == 0 && l2Lambda == 0) {
            return opt;
        }
        return new Optimizer<DifferentiableFunction>() {
            
            @Override
            public boolean minimize(DifferentiableFunction objective, IntDoubleVector point) {
                DifferentiableFunction fn = getRegularizedFn(objective, false, l1Lambda, l2Lambda);
                return opt.minimize(fn, point);
            }
            
        };
    }
232905
def opt(self, x_init, f_fp=None, f=None, fp=None):
        
        tnc_rcstrings = ['Local minimum', 'Converged', 'XConverged', 'Maximum number of f evaluations reached',
             'Line search failed', 'Function is constant']

        assert f_fp != None, "TNC requires f_fp"

        opt_dict = {}
        if self.xtol is not None:
            opt_dict['xtol'] = self.xtol
        if self.ftol is not None:
            opt_dict['ftol'] = self.ftol
        if self.gtol is not None:
            opt_dict['pgtol'] = self.gtol

        opt_result = optimize.fmin_tnc(f_fp, x_init, messages=self.messages,
                       maxfun=self.max_f_eval, **opt_dict)
        self.x_opt = opt_result[0]
        self.f_opt = f_fp(self.x_opt)[0]
        self.funct_eval = opt_result[1]
        self.status = tnc_rcstrings[opt_result[2]]
986146
def convex_conj(self):
        
        if self.scalar <= 0:
            raise ValueError('scaling with nonpositive values have no convex '
                             'conjugate. Current value: {}.'
                             ''.format(self.scalar))

        return self.scalar * self.functional.convex_conj * (1.0 / self.scalar)
262917
def solve_factorized_aug(z, Fval, LU, G, A):
    M, N=G.shape
    P, N=A.shape

    """"""
    m=M    

    """Primal variable"""
    x=z[0:N]

    """Multiplier for equality constraints"""
    nu=z[N:N+P]

    """Multiplier for inequality constraints"""
    l=z[N+P:N+P+M]

    """Slacks"""
    s=z[N+P+M:]

    """Dual infeasibility"""
    rd = Fval[0:N]
    
    """Primal infeasibility"""
    rp1 = Fval[N:N+P]
    rp2 = Fval[N+P:N+P+M]

    """Centrality"""
    rc = Fval[N+P+M:]

    """Sigma matrix"""
    SIG = diags(l/s, 0)

    """LU is actually the augmented system W"""
    W = LU

    b1 = -rd - mydot(G.T, mydot(SIG, rp2)) + mydot(G.T, rc/s)
    b2 = -rp1
    b = np.hstack((b1, b2))

    """Prepare iterative solve via MINRES"""
    sign = np.zeros(N+P)
    sign[0:N//2] = 1.0
    sign[N//2:] = -1.0
    T = diags(sign, 0)        
   
    """Change rhs"""
    b_new = mydot(T, b)

    dW = np.abs(W.diagonal())
    dPc = np.ones(W.shape[0])
    ind = (dW > 0.0)
    dPc[ind] = 1.0/dW[ind]
    Pc = diags(dPc, 0)    
    dxnu, info = minres(W, b_new, tol=1e-10, M=Pc)
    
    # dxnu = solve(J, b)
    dx = dxnu[0:N]
    dnu = dxnu[N:]

    """Obtain search directions for l and s"""
    ds = -rp2 - mydot(G, dx)
    # ds = s*ds
    # SIG = np.diag(l/s)
    dl = -mydot(SIG, ds) - rc/s

    dz = np.hstack((dx, dnu, dl, ds))
    return dz
650218
def solve_factorized_schur(z, Fval, LU, G, A):
    M, N=G.shape
    P, N=A.shape

    """"""
    m = M    

    """Primal variable"""
    x = z[0:N]

    """Multiplier for equality constraints"""
    nu = z[N:N+P]

    """Multiplier for inequality constraints"""
    l = z[N+P:N+P+M]

    """Slacks"""
    s = z[N+P+M:]

    """Dual infeasibility"""
    rd = Fval[0:N]
    
    """Primal infeasibility"""
    rp1 = Fval[N:N+P]
    rp2 = Fval[N+P:N+P+M]

    """Centrality"""
    rc = Fval[N+P+M:]

    """Sigma matrix"""
    SIG = diags(l/s, 0)

    """Assemble right hand side of augmented system"""
    r1 = rd + mydot(G.T, mydot(SIG, rp2)) - mydot(G.T, rc/s)
    r2 = rp1

    """Unpack LU-factors"""
    LU_S, LU_H = LU

    """Assemble right hand side for normal equation"""
    b = r2 - mydot(A, mysolve(LU_H, r1))  

    """Solve for dnu"""
    dnu = mysolve(LU_S, b)
       
    """Solve for dx"""
    dx = mysolve(LU_H, -(r1 + mydot(A.T, dnu)))    
    
    """Obtain search directions for l and s"""
    ds = -rp2 - mydot(G, dx)
    dl = -mydot(SIG, ds) - rc/s

    dz = np.hstack((dx, dnu, dl, ds))
    return dz
1080511
def solve(self, angles0, target):
        """"""
        return self.optimizer.optimize(np.array(angles0), target)
45530
def dumbguy_minimize(f, xmin, xmax, xstep):
    

    prev = f(xmin)
    this = f(xmin+xstep)
    for x in frange(xmin+xstep,xmax,xstep):
        next = f(x+xstep)

        # see if we're on top
        if this < prev and this < next: return x, this

        prev = this
        this = next

    return x, this
827941
public static LinearSolverDense<DMatrixRMaj> general(int numRows , int numCols ) {
        if( numRows == numCols )
            return linear(numRows);
        else
            return leastSquares(numRows,numCols);
    }
131525
def _create_optimizer(self, X, y, status):
        """"""
        if self.optimizer is None:
            self.optimizer = 'avltree'

        times, ranks = y

        if self.optimizer == 'simple':
            optimizer = SimpleOptimizer(X, status, self.alpha, self.rank_ratio, timeit=self.timeit)
        elif self.optimizer == 'PRSVM':
            optimizer = PRSVMOptimizer(X, status, self.alpha, self.rank_ratio, timeit=self.timeit)
        elif self.optimizer == 'direct-count':
            optimizer = LargeScaleOptimizer(self.alpha, self.rank_ratio, self.fit_intercept,
                                            SurvivalCounter(X, ranks, status, len(ranks), times), timeit=self.timeit)
        elif self.optimizer == 'rbtree':
            optimizer = LargeScaleOptimizer(self.alpha, self.rank_ratio, self.fit_intercept,
                                            OrderStatisticTreeSurvivalCounter(X, ranks, status, RBTree, times),
                                            timeit=self.timeit)
        elif self.optimizer == 'avltree':
            optimizer = LargeScaleOptimizer(self.alpha, self.rank_ratio, self.fit_intercept,
                                            OrderStatisticTreeSurvivalCounter(X, ranks, status, AVLTree, times),
                                            timeit=self.timeit)
        else:
            raise ValueError('unknown optimizer: {0}'.format(self.optimizer))

        return optimizer
631772
def minimize_l1(self, weights={}):
        

        if self._z is None:
            self._add_minimization_vars()

        objective = self._z.expr(
            (reaction_id, -weights.get(reaction_id, 1))
            for reaction_id in self._model.reactions)
        self._prob.set_objective(objective)

        self._solve()
128380
def gradient(self, P, Q, Y, i):
        """"""

        return 4 * sum([
            (P[i, j] - Q[i, j]) * (Y[i] - Y[j]) * (1 + np.linalg.norm(Y[i] - Y[j]) ** 2) ** -1 \
            for j in range(Y.shape[0])
        ])
1368470
def proximal(self):
        
        if self.pointwise_norm.exponent == np.inf:
            return proximal_convex_conj_l1(space=self.domain)
        elif self.pointwise_norm.exponent == 2:
            return proximal_convex_conj_l1_l2(space=self.domain)
        else:
            raise NotImplementedError('`proximal` only implemented for p = 1 '
                                      'or 2')
1244593
def _optimize(self, period, tmpid, use_gradient=True):
        """"""
        theta_0 = [self.y.min(), self.y.max() - self.y.min(), 0]
        result = minimize(self._chi2, theta_0, jac=bool(use_gradient),
                          bounds=[(None, None), (0, None), (None, None)],
                          args=(period, tmpid, use_gradient))
        return result.x
199118
def solve_kkt(Q_LU, d, G, A, S_LU, rx, rs, rz, ry):
    
    nineq, nz, neq, nBatch = get_sizes(G, A)

    invQ_rx = rx.btrisolve(*Q_LU)
    if neq > 0:
        h = torch.cat((invQ_rx.unsqueeze(1).bmm(A.transpose(1, 2)).squeeze(1) - ry,
                       invQ_rx.unsqueeze(1).bmm(G.transpose(1, 2)).squeeze(1) + rs / d - rz), 1)
    else:
        h = invQ_rx.unsqueeze(1).bmm(G.transpose(1, 2)).squeeze(1) + rs / d - rz

    w = -(h.btrisolve(*S_LU))

    g1 = -rx - w[:, neq:].unsqueeze(1).bmm(G).squeeze(1)
    if neq > 0:
        g1 -= w[:, :neq].unsqueeze(1).bmm(A).squeeze(1)
    g2 = -rs - w[:, neq:]

    dx = g1.btrisolve(*Q_LU)
    ds = g2 / d
    dz = w[:, neq:]
    dy = w[:, :neq] if neq > 0 else None

    return dx, ds, dz, dy
750377
function constrainedMDSLayout(areas, params) {
        params = params || {};
        var restarts = params.restarts || 10;

        // bidirectionally map sets to a rowid  (so we can create a matrix)
        var sets = [],
            setids = {},
            i;
        for (i = 0; i < areas.length; ++i) {
            var area = areas[i];
            if (area.sets.length == 1) {
                setids[area.sets[0]] = sets.length;
                sets.push(area);
            }
        }

        var matrices = getDistanceMatrices(areas, sets, setids),
            distances = matrices.distances,
            constraints = matrices.constraints;

        // keep distances bounded, things get messed up otherwise.
        // TODO: proper preconditioner?
        var norm = norm2(distances.map(norm2)) / distances.length;
        distances = distances.map(function (row) {
            return row.map(function (value) {
                return value / norm;
            });
        });

        var obj = function obj(x, fxprime) {
            return constrainedMDSGradient(x, fxprime, distances, constraints);
        };

        var best, current;
        for (i = 0; i < restarts; ++i) {
            var initial = zeros(distances.length * 2).map(Math.random);

            current = conjugateGradient(obj, initial, params);
            if (!best || current.fx < best.fx) {
                best = current;
            }
        }
        var positions = best.x;

        // translate rows back to (x,y,radius) coordinates
        var circles = {};
        for (i = 0; i < sets.length; ++i) {
            var set = sets[i];
            circles[set.sets[0]] = {
                x: positions[2 * i] * norm,
                y: positions[2 * i + 1] * norm,
                radius: Math.sqrt(set.size / Math.PI)
            };
        }

        if (params.history) {
            for (i = 0; i < params.history.length; ++i) {
                scale(params.history[i].x, norm);
            }
        }
        return circles;
    }