Python cvxopt.solvers.qp() Examples

def radius(K):
    """evaluate the radius of the MEB (Minimum Enclosing Ball) of examples in
    feature space.

    Parameters
    ----------
    K : (n,n) ndarray,
        the kernel that represents the data.

    Returns
    -------
    r : np.float64,
        the radius of the minimum enclosing ball of examples in feature space.
    """
    K = validation.check_K(K).numpy()
    n = K.shape[0]
    P = 2 * matrix(K)
    p = -matrix(K.diagonal())
    G = -spdiag([1.0] * n)
    h = matrix([0.0] * n)
    A = matrix([1.0] * n).T
    b = matrix([1.0])
    solvers.options['show_progress']=False
    sol = solvers.qp(P,p,G,h,A,b)
    return abs(sol['primal objective'])**.5 

def fit(self, Xs, Xt):
        '''
        Fit source and target using KMM (compute the coefficients)
        :param Xs: ns * dim
        :param Xt: nt * dim
        :return: Coefficients (Pt / Ps) value vector (Beta in the paper)
        '''
        ns = Xs.shape[0]
        nt = Xt.shape[0]
        if self.eps == None:
            self.eps = self.B / np.sqrt(ns)
        K = kernel(self.kernel_type, Xs, None, self.gamma)
        kappa = np.sum(kernel(self.kernel_type, Xs, Xt, self.gamma) * float(ns) / float(nt), axis=1)

        K = matrix(K)
        kappa = matrix(kappa)
        G = matrix(np.r_[np.ones((1, ns)), -np.ones((1, ns)), np.eye(ns), -np.eye(ns)])
        h = matrix(np.r_[ns * (1 + self.eps), ns * (self.eps - 1), self.B * np.ones((ns,)), np.zeros((ns,))])

        sol = solvers.qp(K, -kappa, G, h)
        beta = np.array(sol['x'])
        return beta 

def solve_qp(P, q, G, h,
             A=None, b=None, sym_proj=False,
             solver='cvxopt'):
    if sym_proj:
        P = .5 * (P + P.T)
    cvxmat(P)
    cvxmat(q)
    cvxmat(G)
    cvxmat(h)
    args = [cvxmat(P), cvxmat(q), cvxmat(G), cvxmat(h)]
    if A is not None:
        args.extend([cvxmat(A), cvxmat(b)])
    sol = qp(*args, solver=solver)
    if not ('optimal' in sol['status']):
        raise ValueError('QP optimum not found: %s' % sol['status'])
    return np.array(sol['x']).reshape((P.shape[1],)) 

def optimize_portfolio(n, avg_ret, covs, r_min):
	P = covs
	# x = variable(n)
	q = matrix(numpy.zeros((n, 1)), tc='d')
	# inequality constraints Gx <= h
	# captures the constraints (avg_ret'x >= r_min) and (x >= 0)
	G = matrix(numpy.concatenate((
		-numpy.transpose(numpy.array(avg_ret)),
		-numpy.identity(n)), 0))
	h = matrix(numpy.concatenate((
		-numpy.ones((1,1))*r_min,
		numpy.zeros((n,1))), 0))
	# equality constraint Ax = b; captures the constraint sum(x) == 1
	A = matrix(1.0, (1,n))
	b = matrix(1.0)
	sol = solvers.qp(P, q, G, h, A, b)
	return sol

### setup the parameters 

def projection_in_norm(self, x, M):
        """
        Projection of x to simplex induced by matrix M. Uses quadratic programming.
        """
        m = M.shape[0]

        # Constrains matrices
        P = opt.matrix(2 * M)
        q = opt.matrix(-2 * M * x)
        G = opt.matrix(-np.eye(m))
        h = opt.matrix(np.zeros((m, 1)))
        A = opt.matrix(np.ones((1, m)))
        b = opt.matrix(1.)

        # Solve using quadratic programming
        sol = opt.solvers.qp(P, q, G, h, A, b)
        return np.squeeze(sol['x']) 

def update_w(self):
        """ alternating least squares step, update W under the convexity
        constraint """
        def update_single_w(i):
            """ compute single W[:,i] """
            # optimize beta     using qp solver from cvxopt
            FB = base.matrix(np.float64(np.dot(-self.data.T, W_hat[:,i])))
            be = solvers.qp(HB, FB, INQa, INQb, EQa, EQb)
            self.beta[i,:] = np.array(be['x']).reshape((1, self._num_samples))

        # float64 required for cvxopt
        HB = base.matrix(np.float64(np.dot(self.data[:,:].T, self.data[:,:])))
        EQb = base.matrix(1.0, (1, 1))
        W_hat = np.dot(self.data, pinv(self.H))
        INQa = base.matrix(-np.eye(self._num_samples))
        INQb = base.matrix(0.0, (self._num_samples, 1))
        EQa = base.matrix(1.0, (1, self._num_samples))

        for i in range(self._num_bases):
            update_single_w(i)

        self.W = np.dot(self.beta, self.data.T).T 

def update_h(self):
        """ alternating least squares step, update H under the convexity
        constraint """
        def update_single_h(i):
            """ compute single H[:,i] """
            # optimize alpha using qp solver from cvxopt
            FA = base.matrix(np.float64(np.dot(-self.W.T, self.data[:,i])))
            al = solvers.qp(HA, FA, INQa, INQb, EQa, EQb)
            self.H[:,i] = np.array(al['x']).reshape((1, self._num_bases))

        EQb = base.matrix(1.0, (1,1))
        # float64 required for cvxopt
        HA = base.matrix(np.float64(np.dot(self.W.T, self.W)))
        INQa = base.matrix(-np.eye(self._num_bases))
        INQb = base.matrix(0.0, (self._num_bases,1))
        EQa = base.matrix(1.0, (1, self._num_bases))

        for i in range(self._num_samples):
            update_single_h(i) 

def _QP(self, x, y):
        # In QP formulation (dual): m variables, 2m+1 constraints (1 equation, 2m inequations)
        m = len(y)
        print x.shape, y.shape
        P = self._kernel(x) * np.outer(y, y)
        P, q = matrix(P, tc='d'), matrix(-np.ones((m, 1)), tc='d')
        G = matrix(np.r_[-np.eye(m), np.eye(m)], tc='d')
        h = matrix(np.r_[np.zeros((m,1)), np.zeros((m,1)) + self.C], tc='d')
        A, b = matrix(y.reshape((1,-1)), tc='d'), matrix([0.0])
        # print "P, q:"
        # print P, q
        # print "G, h"
        # print G, h
        # print "A, b"
        # print A, b
        solution = solvers.qp(P, q, G, h, A, b)
        if solution['status'] == 'unknown':
            print 'Not PSD!'
            exit(2)
        else:
            self.alphas = np.array(solution['x']).squeeze() 

def projection_in_norm(self, x, M):
        """ Projection of x to simplex indiced by matrix M. Uses quadratic programming.
        """
        m = M.shape[0]

        P = matrix(2*M)
        q = matrix(-2 * M * x)
        G = matrix(-np.eye(m))
        h = matrix(np.zeros((m,1)))
        A = matrix(np.ones((1,m)))
        b = matrix(1.)

        sol = solvers.qp(P, q, G, h, A, b)
        return np.squeeze(sol['x']) 

def calculate(self, x_fit: np.ndarray) -> np.ndarray:
        if x_fit.ndim == 1:
            x_fit = x_fit.reshape(x_fit.shape[0], 1)
        solvers.options['show_progress'] = False

        M = self.X.collapse()

        N, p1 = M.shape
        nvars, p2 = x_fit.T.shape
        C = _numpy_to_cvxopt_matrix(x_fit)
        Q = C.T * C

        lb_A = -np.eye(nvars)
        lb = np.repeat(0, nvars)
        A = _numpy_None_vstack(None, lb_A)
        b = _numpy_None_concatenate(None, -lb)
        A = _numpy_to_cvxopt_matrix(A)
        b = _numpy_to_cvxopt_matrix(b)

        Aeq = _numpy_to_cvxopt_matrix(np.ones((1, nvars)))
        beq = _numpy_to_cvxopt_matrix(np.ones(1))

        M = np.array(M, dtype=np.float64)
        self.map = np.zeros((N, nvars), dtype=np.float32)
        for n1 in range(N):
            d = matrix(M[n1], (p1, 1), 'd')
            q = - d.T * C
            sol = solvers.qp(Q, q.T, A, b, Aeq, beq, None, None)['x']
            self.map[n1] = np.array(sol).squeeze()
        self.map = self.map.reshape(self.X.shape[:-1] + (x_fit.shape[-1],))

        return self.map 

def fit(self, X, V, k, s=None, method="qp", eps=None, tol=1e-4):  # pass index
        # the parameter k will be replaced by a connectivity matrix in the future.
        self.__reset__()
        # knn clustering
        nbrs = NearestNeighbors(n_neighbors=k, algorithm="ball_tree").fit(X)
        _, Idx = nbrs.kneighbors(X)
        # compute transition prob.
        n = X.shape[0]
        self.P = np.zeros((n, n))
        if method == "kernel":
            inv_s = np.linalg.inv(s)
            # compute density kernel
            if eps is not None:
                self.Kd = np.zeros((n, n))
                inv_eps = 1 / eps
                for i in range(n):
                    self.Kd[i, Idx[i]] = compute_density_kernel(
                        X[i], X[Idx[i]], inv_eps
                    )
                D = np.sum(self.Kd, 0)
        for i in range(n):
            y = X[i]
            v = V[i]
            if method == "qp":
                Y = X[Idx[i, 1:]]
                p = compute_markov_trans_prob(y, v, Y, s)
                p[p <= tol] = 0  # tolerance check
                self.P[Idx[i, 1:], i] = p
                self.P[i, i] = 1 - np.sum(p)
            else:
                Y = X[Idx[i]]
                # p = compute_kernel_trans_prob(y, v, Y, inv_s)
                k = compute_drift_kernel(y, v, Y, inv_s)
                if eps is not None:
                    k /= D[Idx[i]]
                p = k / np.sum(k)
                p[p <= tol] = 0  # tolerance check
                p = p / np.sum(p)
                self.P[Idx[i], i] = p 

def compute_markov_trans_prob(x, v, X, s=None, cont_time=False):
    from cvxopt import matrix, solvers

    n = X.shape[0]
    R = X - x
    # normalize R, v, and s
    Rn = np.array(R, copy=True)
    vn = np.array(v, copy=True)
    scale = np.abs(np.max(R, 0) - np.min(R, 0))
    Rn = Rn / scale
    vn = vn / scale
    if s is not None:
        sn = np.array(s, copy=True)
        sn = sn / scale
        A = np.hstack((Rn, 0.5 * Rn * Rn))
        b = np.hstack((vn, 0.5 * sn * sn))
    else:
        A = Rn
        b = vn

    H = A.dot(A.T)
    f = b.dot(A.T)
    if cont_time:
        G = -np.eye(n)
        h = np.zeros(n)
    else:
        G = np.vstack((-np.eye(n), np.ones(n)))
        h = np.zeros(n + 1)
        h[-1] = 1
    p = solvers.qp(matrix(H), matrix(-f), G=matrix(G), h=matrix(h))["x"]
    p = np.array(p).flatten()
    return p 

def markov_combination(x, v, X):
    from cvxopt import matrix, solvers

    n = X.shape[0]
    R = matrix(X - x).T
    H = R.T * R
    f = matrix(v).T * R
    G = np.vstack((-np.eye(n), np.ones(n)))
    h = np.zeros(n + 1)
    h[-1] = 1
    p = solvers.qp(H, -f.T, G=matrix(G), h=matrix(h))["x"]
    u = R * p
    return p, u 

def fit(self, x, y):
        from cvxopt import matrix, solvers
        solvers.options['show_progress'] = False

        check_classification_targets(y)
        x, y = check_X_y(x, y)
        x_s, x_u = x[y == +1, :], x[y == 0, :]
        n_s, n_u = len(x_s), len(x_u)

        p_p = self.prior
        p_n = 1 - self.prior
        p_s = p_p ** 2 + p_n ** 2
        k_s = self._basis(x_s)
        k_u = self._basis(x_u)
        d = k_u.shape[1]

        P = np.zeros((d + 2 * n_u, d + 2 * n_u))
        P[:d, :d] = self.lam * np.eye(d)
        q = np.vstack((
            -p_s / (n_s * (p_p - p_n)) * k_s.T.dot(np.ones((n_s, 1))),
            -p_n / (n_u * (p_p - p_n)) * np.ones((n_u, 1)),
            -p_p / (n_u * (p_p - p_n)) * np.ones((n_u, 1))
        ))
        G = np.vstack((
            np.hstack((np.zeros((n_u, d)), -np.eye(n_u), np.zeros((n_u, n_u)))),
            np.hstack((0.5 * k_u, -np.eye(n_u), np.zeros((n_u, n_u)))),
            np.hstack((k_u, -np.eye(n_u), np.zeros((n_u, n_u)))),
            np.hstack((np.zeros((n_u, d)), np.zeros((n_u, n_u)), -np.eye(n_u))),
            np.hstack((-0.5 * k_u, np.zeros((n_u, n_u)), -np.eye(n_u))),
            np.hstack((-k_u, np.zeros((n_u, n_u)), -np.eye(n_u)))
        ))
        h = np.vstack((
            np.zeros((n_u, 1)),
            -0.5 * np.ones((n_u, 1)),
            np.zeros((n_u, 1)),
            np.zeros((n_u, 1)),
            -0.5 * np.ones((n_u, 1)),
            np.zeros((n_u, 1))
        ))
        sol = solvers.qp(matrix(P), matrix(q), matrix(G), matrix(h))
        self.coef_ = np.array(sol['x'])[:d] 

def update_w(self):
        def updatesingleW(i):
        # optimize alpha using qp solver from cvxopt
            FA = base.matrix(np.float64(np.dot(-self.H, self.data[i,:].T)))
            al = solvers.qp(HA, FA, INQa, INQb)
            self.W[i,:] = np.array(al['x']).reshape((1,-1))

        # float64 required for cvxopt
        HA = base.matrix(np.float64(np.dot(self.H, self.H.T)))
        INQa = base.matrix(-np.eye(self._num_bases))
        INQb = base.matrix(0.0, (self._num_bases,1))

        map(updatesingleW, range(self._data_dimension)) 

def update_h(self):
        def updatesingleH(i):
            # optimize alpha using qp solver from cvxopt
            FA = base.matrix(np.float64(np.dot(-self.W.T, self.data[:,i])))
            al = solvers.qp(HA, FA, INQa, INQb)
            self.H[:,i] = np.array(al['x']).reshape((1,-1))

        # float64 required for cvxopt
        HA = base.matrix(np.float64(np.dot(self.W.T, self.W)))
        INQa = base.matrix(-np.eye(self._num_bases))
        INQb = base.matrix(0.0, (self._num_bases,1))

        map(updatesingleH, range(self._num_samples)) 

def fit_nnl2reg(self):
        try:
            from cvxopt import matrix, solvers
        except ImportError:
            raise ImportError("To infer titer models, you need a working installation of cvxopt")
        n_params = self.design_matrix.shape[1]
        P = matrix(np.dot(self.design_matrix.T, self.design_matrix) + self.lam_drop*np.eye(n_params))
        q = matrix( -np.dot( self.titer_dist, self.design_matrix))
        h = matrix(np.zeros(n_params)) # Gw <=h
        G = matrix(-np.eye(n_params))
        W = solvers.qp(P,q,G,h)
        return np.array([x for x in W['x']]) 

def _fit(self,X,Y):    
        self.X = X
        values = np.unique(Y)
        Y = [1 if l==values[1] else -1 for l in Y]
        self.Y = Y
        npos = len([1.0 for l in Y if l == 1])
        nneg = len([1.0 for l in Y if l == -1])
        gamma_unif = matrix([1.0/npos if l == 1 else 1.0/nneg for l in Y])
        YY = matrix(np.diag(list(matrix(Y))))

        Kf = self.__kernel_definition__()
        ker_matrix = matrix(Kf(X,X).astype(np.double))
        #KLL = (1.0 / (gamma_unif.T * YY * ker_matrix * YY * gamma_unif)[0])*(1.0-self.lam)*YY*ker_matrix*YY
        KLL = (1.0-self.lam)*YY*ker_matrix*YY
        LID = matrix(np.diag([self.lam * (npos * nneg / (npos+nneg))]*len(Y)))
        Q = 2*(KLL+LID)
        p = matrix([0.0]*len(Y))
        G = -matrix(np.diag([1.0]*len(Y)))
        h = matrix([0.0]*len(Y),(len(Y),1))
        A = matrix([[1.0 if lab==+1 else 0 for lab in Y],[1.0 if lab2==-1 else 0 for lab2 in Y]]).T
        b = matrix([[1.0],[1.0]],(2,1))
        
        solvers.options['show_progress'] = False#True
        solvers.options['maxiters'] = self.max_iter
        sol = solvers.qp(Q,p,G,h,A,b)
        self.gamma = sol['x']
        if self.verbose:
            print ('[KOMD]')
            print ('optimization finished, #iter = %d' % sol['iterations'])
            print ('status of the solution: %s' % sol['status'])
            print ('objval: %.5f' % sol['primal objective'])
            
        bias = 0.5 * self.gamma.T * ker_matrix * YY * self.gamma
        self.bias = bias
        self.is_fitted = True
        self.ker_matrix = ker_matrix
        return self 

def opt_margin(K,YY,init_sol=None):
	'''optimized margin evaluation'''
	n = K.shape[0]
	P = 2 * (YY * matrix(K) * YY)
	p = matrix([0.0]*n)
	G = -spdiag([1.0]*n)
	h = matrix([0.0]*n)
	A = matrix([[1.0 if YY[i,i]==+1 else 0 for i in range(n)],
				[1.0 if YY[j,j]==-1 else 0 for j in range(n)]]).T
	b = matrix([[1.0],[1.0]],(2,1))
	solvers.options['show_progress']=False
	sol = solvers.qp(P,p,G,h,A,b,initvals=init_sol)	
	margin2 = sol['primal objective']
	return margin2, sol['x'], sol 

def opt_margin(K, YY, init_sol=None):
    '''optimized margin evaluation'''
    n = K.shape[0]
    P = 2 * (YY * matrix(K.numpy()) * YY)
    p = matrix([0.0]*n)
    G = -spdiag([1.0]*n)
    h = matrix([0.0]*n)
    A = matrix([[1.0 if YY[i,i]==+1 else 0 for i in range(n)],
                [1.0 if YY[j,j]==-1 else 0 for j in range(n)]]).T
    b = matrix([[1.0],[1.0]],(2,1))
    solvers.options['show_progress']=False
    sol = solvers.qp(P,p,G,h,A,b,initvals=init_sol) 
    margin2 = sol['primal objective']
    return sol, margin2 

def opt_radius(K, init_sol=None): 
    n = K.shape[0]
    K = matrix(K.numpy())
    P = 2 * K
    p = -matrix([K[i,i] for i in range(n)])
    G = -spdiag([1.0] * n)
    h = matrix([0.0] * n)
    A = matrix([1.0] * n).T
    b = matrix([1.0])
    solvers.options['show_progress']=False
    sol = solvers.qp(P,p,G,h,A,b,initvals=init_sol)
    radius2 = (-p.T * sol['x'])[0] - (sol['x'].T * K * sol['x'])[0]
    return sol, radius2 

def margin(K,Y):
    """evaluate the margin in a classification problem of examples in feature space.
    If the classes are not linearly separable in feature space, then the
    margin obtained is 0.

    Note that it works only for binary tasks.

    Parameters
    ----------
    K : (n,n) ndarray,
        the kernel that represents the data.
    Y : (n) array_like,
        the labels vector.
    """
    K, Y = validation.check_K_Y(K, Y, binary=True)
    n = Y.size()[0]
    Y = [1 if y==Y[0] else -1 for y in Y]
    YY = spdiag(Y)
    P = 2*(YY*matrix(K.numpy())*YY)
    p = matrix([0.0]*n)
    G = -spdiag([1.0]*n)
    h = matrix([0.0]*n)
    A = matrix([[1.0 if Y[i]==+1 else 0 for i in range(n)],
                [1.0 if Y[j]==-1 else 0 for j in range(n)]]).T
    b = matrix([[1.0],[1.0]],(2,1))
    solvers.options['show_progress']=False
    sol = solvers.qp(P,p,G,h,A,b)
    return sol['primal objective']**.5 

def iwe_kernel_mean_matching(self, X, Z):
        """
        Estimate importance weights based on kernel mean matching.

        Parameters
        ----------
        X : array
            source data (N samples by D features)
        Z : array
            target data (M samples by D features)

        Returns
        -------
        iw : array
            importance weights (N samples by 1)

        """
        # Data shapes
        N, DX = X.shape
        M, DZ = Z.shape

        # Assert equivalent dimensionalities
        if not DX == DZ:
            raise ValueError('Dimensionalities of X and Z should be equal.')

        # Compute sample pairwise distances
        KXX = cdist(X, X, metric='euclidean')
        KXZ = cdist(X, Z, metric='euclidean')

        # Check non-negative distances
        if not np.all(KXX >= 0):
            raise ValueError('Non-positive distance in source kernel.')
        if not np.all(KXZ >= 0):
            raise ValueError('Non-positive distance in source-target kernel.')

        # Compute kernels
        if self.kernel_type == 'rbf':
            # Radial basis functions
            KXX = np.exp(-KXX / (2*self.bandwidth**2))
            KXZ = np.exp(-KXZ / (2*self.bandwidth**2))

        # Collapse second kernel and normalize
        KXZ = N/M * np.sum(KXZ, axis=1)

        # Prepare for CVXOPT
        Q = matrix(KXX, tc='d')
        p = matrix(KXZ, tc='d')
        G = matrix(np.concatenate((np.ones((1, N)), -1*np.ones((1, N)),
                                   -1.*np.eye(N)), axis=0), tc='d')
        h = matrix(np.concatenate((np.array([N/np.sqrt(N) + N], ndmin=2),
                                   np.array([N/np.sqrt(N) - N], ndmin=2),
                                   np.zeros((N, 1))), axis=0), tc='d')

        # Call quadratic program solver
        sol = solvers.qp(Q, p, G, h)

        # Return optimal coefficients as importance weights
        return np.array(sol['x'])[:, 0] 

def fit_nnl1reg(self):
        '''l1 regularization of titer drops with non-negativity constraints

        Returns
        -------
        TYPE
            Description
        '''
        try:
            from cvxopt import matrix, solvers
        except ImportError:
            raise ImportError("To infer titer models, you need a working installation of cvxopt")
        n_params = self.design_matrix.shape[1]
        n_genetic = self.genetic_params
        n_sera = len(self.sera)
        n_v = len(self.test_strains)

        # set up the quadratic matrix containing the deviation term (linear xterm below)
        # and the l2-regulatization of the avidities and potencies
        P1 = np.zeros((n_params,n_params))
        P1[:n_params, :n_params] = self.TgT
        for ii in range(n_genetic, n_genetic+n_sera):
            P1[ii,ii]+=self.lam_pot
        for ii in range(n_genetic+n_sera, n_params):
            P1[ii,ii]+=self.lam_avi
        P = matrix(P1)

        # set up cost for auxillary parameter and the linear cross-term
        q1 = np.zeros(n_params)
        q1[:n_params] = -np.dot(self.titer_dist, self.design_matrix)
        q1[:n_genetic] += self.lam_drop
        q = matrix(q1)

        # set up linear constraint matrix to enforce positivity of the
        # dTiters and bounding of dTiter by the auxillary parameter
        h = matrix(np.zeros(n_genetic))     # Gw <=h
        G1 = np.zeros((n_genetic,n_params))
        G1[:n_genetic, :n_genetic] = -np.eye(n_genetic)
        G = matrix(G1)
        W = solvers.qp(P,q,G,h)
        return np.array([x for x in W['x']])[:n_params]

##########################################################################################
# END GENERIC CLASS
##########################################################################################



##########################################################################################
# TREE MODEL
########################################################################################## 

def fit_l1reg(self):
        '''
        regularize genetic parameters with an l1 norm regardless of sign

        Returns
        -------
        TYPE
            Description
        '''
        try:
            from cvxopt import matrix, solvers
        except ImportError:
            raise ImportError("To infer titer models, you need a working installation of cvxopt")
        n_params = self.design_matrix.shape[1]
        n_genetic = self.genetic_params
        n_sera = len(self.sera)
        n_v = len(self.test_strains)

        # set up the quadratic matrix containing the deviation term (linear xterm below)
        # and the l2-regulatization of the avidities and potencies
        P1 = np.zeros((n_params+n_genetic,n_params+n_genetic))
        P1[:n_params, :n_params] = self.TgT
        for ii in range(n_genetic, n_genetic+n_sera):
            P1[ii,ii]+=self.lam_pot
        for ii in range(n_genetic+n_sera, n_params):
            P1[ii,ii]+=self.lam_avi
        P = matrix(P1)

        # set up cost for auxillary parameter and the linear cross-term
        q1 = np.zeros(n_params+n_genetic)
        q1[:n_params] = -np.dot( self.titer_dist, self.design_matrix)
        q1[n_params:] = self.lam_drop
        q = matrix(q1)

        # set up linear constraint matrix to regularize the HI parametesr
        h = matrix(np.zeros(2*n_genetic))   # Gw <=h
        G1 = np.zeros((2*n_genetic,n_params+n_genetic))
        G1[:n_genetic, :n_genetic] = -np.eye(n_genetic)
        G1[:n_genetic:, n_params:] = -np.eye(n_genetic)
        G1[n_genetic:, :n_genetic] = np.eye(n_genetic)
        G1[n_genetic:, n_params:] = -np.eye(n_genetic)
        G = matrix(G1)
        W = solvers.qp(P,q,G,h)
        return np.array([x for x in W['x']])[:n_params] 

def update(self, cov_mat, exp_rets):
        """
         Note: As the Sharpe ratio is not invariant with respect
         to leverage, it is not possible to construct non-trivial
         market neutral tangency portfolios. This is because for
         a positive initial Sharpe ratio the sharpe grows unbound
         with increasing leverage.

         Parameters
         ----------
         cov_mat: pandas.DataFrame
             Covariance matrix of asset returns.
         exp_rets: pandas.Series
             Expected asset returns (often historical returns).
         allow_short: bool, optional
             If 'False' construct a long-only portfolio.
             If 'True' allow shorting, i.e. negative weights.

         Returns
         -------
         weights: pandas.Series
             Optimal asset weights.
         """
        if not isinstance(cov_mat, pd.DataFrame):
            raise ValueError("Covariance matrix is not a DataFrame")

        if not isinstance(exp_rets, pd.Series):
            raise ValueError("Expected returns is not a Series")

        if not cov_mat.index.equals(exp_rets.index):
            raise ValueError("Indices do not match")

        n = len(cov_mat)

        P = opt.matrix(cov_mat.values)
        q = opt.matrix(0.0, (n, 1))

        # Constraints Gx <= h
        # exp_rets*x >= 1 and x >= 0
        G = opt.matrix(np.vstack((-exp_rets.values,
                                  -np.identity(n))))
        h = opt.matrix(np.vstack((-1.0,
                                  np.zeros((n, 1)))))

        # Solve
        optsolvers.options['show_progress'] = False
        sol = optsolvers.qp(P, q, G, h)

        if sol['status'] != 'optimal':
            warnings.warn("Convergence problem")

        weights = np.append(np.squeeze(sol['x']), [0.0])

        # Rescale weights, so that sum(weights) = 1
        weights /= weights.sum()
        return weights 

def cvxopt_solve_qp(P, q, G=None, h=None, A=None, b=None, solver=None,
                    initvals=None, verbose=False):
    """
    Solve a Quadratic Program defined as:

    .. math::

        \\begin{split}\\begin{array}{ll}
        \\mbox{minimize} &
            \\frac{1}{2} x^T P x + q^T x \\\\
        \\mbox{subject to}
            & G x \\leq h                \\\\
            & A x = h
        \\end{array}\\end{split}

    using `CVXOPT <http://cvxopt.org/>`_.

    Parameters
    ----------
    P : numpy.array, cvxopt.matrix or cvxopt.spmatrix
        Symmetric quadratic-cost matrix.
    q : numpy.array, cvxopt.matrix or cvxopt.spmatrix
        Quadratic-cost vector.
    G : numpy.array, cvxopt.matrix or cvxopt.spmatrix
        Linear inequality matrix.
    h : numpy.array, cvxopt.matrix or cvxopt.spmatrix
        Linear inequality vector.
    A : numpy.array, cvxopt.matrix or cvxopt.spmatrix
        Linear equality constraint matrix.
    b : numpy.array, cvxopt.matrix or cvxopt.spmatrix
        Linear equality constraint vector.
    solver : string, optional
        Set to 'mosek' to run MOSEK rather than CVXOPT.
    initvals : numpy.array, optional
        Warm-start guess vector.
    verbose : bool, optional
        Set to `True` to print out extra information.

    Returns
    -------
    x : array, shape=(n,)
        Solution to the QP, if found, otherwise ``None``.

    Note
    ----
    CVXOPT only considers the lower entries of `P`, therefore it will use a
    wrong cost function if a non-symmetric matrix is provided.
    """
    options['show_progress'] = verbose
    args = [cvxopt_matrix(P), cvxopt_matrix(q)]
    kwargs = {'G': None, 'h': None, 'A': None, 'b': None}
    if G is not None:
        kwargs['G'] = cvxopt_matrix(G)
        kwargs['h'] = cvxopt_matrix(h)
    if A is not None:
        kwargs['A'] = cvxopt_matrix(A)
        kwargs['b'] = cvxopt_matrix(b)
    sol = qp(*args, solver=solver, initvals=initvals, **kwargs)
    if 'optimal' not in sol['status']:
        return None
    return array(sol['x']).reshape((q.shape[0],)) 

def min_var_portfolio(cov_mat, allow_short=False):
    """
    Computes the minimum variance portfolio.

    Note: As the variance is not invariant with respect
    to leverage, it is not possible to construct non-trivial
    market neutral minimum variance portfolios. This is because
    the variance approaches zero with decreasing leverage,
    i.e. the market neutral portfolio with minimum variance
    is not invested at all.
    
    Parameters
    ----------
    cov_mat: pandas.DataFrame
        Covariance matrix of asset returns.
    allow_short: bool, optional
        If 'False' construct a long-only portfolio.
        If 'True' allow shorting, i.e. negative weights.

    Returns
    -------
    weights: pandas.Series
        Optimal asset weights.
    """
    if not isinstance(cov_mat, pd.DataFrame):
        raise ValueError("Covariance matrix is not a DataFrame")

    n = len(cov_mat)

    P = opt.matrix(cov_mat.values)
    q = opt.matrix(0.0, (n, 1))

    # Constraints Gx <= h
    if not allow_short:
        # x >= 0
        G = opt.matrix(-np.identity(n))
        h = opt.matrix(0.0, (n, 1))
    else:
        G = None
        h = None

    # Constraints Ax = b
    # sum(x) = 1
    A = opt.matrix(1.0, (1, n))
    b = opt.matrix(1.0)

    # Solve
    optsolvers.options['show_progress'] = False
    sol = optsolvers.qp(P, q, G, h, A, b)

    if sol['status'] != 'optimal':
        warnings.warn("Convergence problem")

    # Put weights into a labeled series
    weights = pd.Series(sol['x'], index=cov_mat.index)
    return weights