Python cvxopt.spdiag() 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 initialize_optimization(self):
        YY          = spdiag([1 if y==self.Y[0] else -1 for y in self.Y])
        weights     = uniform_vector(self.n_kernels)
        ker_matrix  = self.func_form(self.KL, weights)
        alpha,r2    = opt_radius(ker_matrix)
        gamma,m2    = opt_margin(ker_matrix, YY)
        obj         = (r2 / m2) / len(self.Y)

        #caching
        self.cache.YY = YY
        self.cache.alpha = alpha
        self.cache.gamma = gamma

        return Solution(
            weights=weights, 
            objective=obj,
            ker_matrix=ker_matrix,
            ) 

def diag(self, dim):
        if self.size != (1, 1):
            raise Exception('not implemented')
        selfcopy = self.copy()
        idx = cvx.spdiag([1.] * dim)[:].I

        for k in self.factors:
            tc = 'z' if self.factors[k].typecode=='z' else 'd'
            selfcopy.factors[k] = spmatrix(
                [], [], [], (dim**2, self.factors[k].size[1]),tc=tc)
            for i in idx:
                selfcopy.factors[k][i, :] = self.factors[k]
        if not self.constant is None:
            tc = 'z' if self.constant.typecode=='z' else 'd'
            selfcopy.constant = cvx.matrix(0., (dim**2, 1),tc=tc)
            for i in idx:
                selfcopy.constant[i] = self.constant[0]
        else:
            selfcopy.constant = None
        selfcopy._size = (dim, dim)
        selfcopy.string = 'diag(' + selfcopy.string + ')'
        return selfcopy 

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 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 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_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 _combine_kernels(self):
        assert len(self.Y.unique()) == 2
        Y = [1 if y==self.classes_[1] else -1 for y in self.Y]
        n_sample = len(self.Y)
        ker_matrix = matrix(self.func_form(self.KL).numpy())
        YY = spdiag(Y)
        #KLL = (1.0-self.lam)*YY*ker_matrix*YY
        #LID = spdiag([self.lam]*n_sample)
        #Q = 2*(KLL+LID)
        Q = 2 * ((1.0-self.lam)*YY*ker_matrix*YY + spdiag([self.lam]*n_sample))
        p = matrix([0.0]*n_sample)
        G = -spdiag([1.0]*n_sample)
        h = matrix([0.0]*n_sample,(n_sample,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
        solvers.options['maxiters'] = 200
        sol = solvers.qp(Q,p,G,h,A,b)
        gamma = sol['x']

        yg = gamma.T * YY
        weights = [(yg*matrix(K.numpy())*yg.T)[0] for K in self.KL]
         
        norm2 = sum([w for w in weights])
        
        weights = torch.tensor([w / norm2 for w in weights])
        ker_matrix = self.func_form(self.KL, weights)
        return Solution(
            weights=weights,
            objective=None,
            ker_matrix=ker_matrix,
            ) 

def diag_vect(exp):
    """
    Returns the vector with the diagonal elements of the matrix expression ``exp``

    **Example**

    >>> import picos as pic
    >>> prob=pic.Problem()
    >>> X=prob.add_variable('X',(3,3))
    >>> pic.tools.diag_vect(X)
    # (3 x 1)-affine expression: diag(X) #

    """
    from .expression import AffinExp
    (n, m) = exp.size
    n = min(n, m)
    idx = cvx.spdiag([1.] * n)[:].I
    expcopy = AffinExp(exp.factors.copy(), exp.constant, exp.size,
                       exp.string)
    proj = spmatrix([1.] * n, range(n), idx,
                        (n, exp.size[0] * exp.size[1]))
    for k in exp.factors.keys():
        expcopy.factors[k] = proj * expcopy.factors[k]
    if not exp.constant is None:
        expcopy.constant = proj * expcopy.constant
    expcopy._size = (n, 1)
    expcopy.string = 'diag(' + exp.string + ')'
    return expcopy 

def eval(self, ind=None):
        val = self.exp.eval(ind)

        if not isinstance(val, cvx.base.matrix):
            val = cvx.matrix(val)
        p = float(self.numerator) / float(self.denominator)
        if self.M is None:
            ev = np.linalg.eigvalsh(np.matrix(val))
            return sum([vi**p for vi in ev])
        else:
            Mval = self.M.eval(ind)
            U, S, V = np.linalg.svd(val)
            Xp = cvx.matrix(U) * cvx.spdiag([s**p for s in S]) * cvx.matrix(V)
            return np.trace(Mval * Xp) 

def diag(exp, dim=1):
    r"""
    if ``exp`` is an affine expression of size (n,m),
    ``diag(exp,dim)`` returns a diagonal matrix of size ``dim*n*m`` :math:`\times` ``dim*n*m``,
    with ``dim`` copies of the vectorized expression ``exp[:]`` on the diagonal.

    In particular:

      * when ``exp`` is scalar, ``diag(exp,n)`` returns a diagonal
        matrix of size :math:`n \times n`, with all diagonal elements equal to ``exp``.

      * when ``exp`` is a vector of size :math:`n`, ``diag(exp)`` returns the diagonal
        matrix of size :math:`n \times n` with the vector ``exp`` on the diagonal


    **Example**

    >>> import picos as pic
    >>> prob=pic.Problem()
    >>> x=prob.add_variable('x',1)
    >>> y=prob.add_variable('y',1)
    >>> pic.tools.diag(x-y,4)
    # (4 x 4)-affine expression: Diag(x -y) #
    >>> pic.tools.diag(x//y)
    # (2 x 2)-affine expression: Diag([x;y]) #

    """
    from .expression import AffinExp
    if not isinstance(exp, AffinExp):
        mat, name = _retrieve_matrix(exp)
        exp = AffinExp({}, constant=mat[:], size=mat.size, string=name)
    (n, m) = exp.size
    expcopy = AffinExp(exp.factors.copy(), exp.constant, exp.size,
                       exp.string)
    idx = cvx.spdiag([1.] * dim * n * m)[:].I
    for k in exp.factors.keys():
        # ensure it's sparse
        mat = cvx.sparse(expcopy.factors[k])
        I, J, V = list(mat.I), list(mat.J), list(mat.V)
        newI = []
        for d in range(dim):
            for i in I:
                newI.append(idx[i + n * m * d])
        expcopy.factors[k] = spmatrix(
            V * dim, newI, J * dim, ((dim * n * m)**2, exp.factors[k].size[1]))
    expcopy.constant = cvx.matrix(0., ((dim * n * m)**2, 1))
    if not exp.constant is None:
        for k, i in enumerate(idx):
            expcopy.constant[i] = exp.constant[k % (n * m)]
    expcopy._size = (dim * n * m, dim * n * m)
    expcopy.string = 'Diag(' + exp.string + ')'
    return expcopy 

def __xor__(self, fact):
        """hadamard (elementwise) product"""
        selfcopy = self.copy()
        if isinstance(fact, AffinExp):
            if fact.isconstant():
                fac, facString = cvx.sparse(fact.eval()), fact.string
            else:
                if self.isconstant():
                    return fact ^ self
                else:
                    raise Exception('not implemented')
        else:
            fac, facString = _retrieve_matrix(fact, self.size[0])
        if fac.size == (1, 1) and selfcopy.size[0] != 1:
            fac = fac[0] * cvx.spdiag([1.] * selfcopy.size[0])
        if self.size == (1, 1) and fac.size[1] != 1:
            oldstring = selfcopy.string
            selfcopy = selfcopy.diag(fac.size[1])
            selfcopy.string = oldstring
        if selfcopy.size[0] != fac.size[0] or selfcopy.size[1] != fac.size[1]:
            raise Exception('incompatible dimensions')
        mm, nn = selfcopy.size
        bfac = spmatrix([], [], [], (mm * nn, mm * nn))
        for i, j, v in zip(fac.I, fac.J, fac.V):
            bfac[j * mm + i, j * mm + i] = v
        for k in selfcopy.factors:
            newfac = bfac * selfcopy.factors[k]
            selfcopy.factors[k] = newfac
        if selfcopy.constant is None:
            newfac = None
        else:
            newfac = bfac * selfcopy.constant
        selfcopy.constant = newfac
        """
                #the following removes 'I' from the string when a matrix is multiplied
                #by the identity. We leave the 'I' when the factor of identity is a scalar
                if len(facString)>0:
                        if facString[-1]=='I' and (len(facString)==1
                                 or facString[-2].isdigit() or facString[-2]=='.') and (
                                 self.size != (1,1)):
                                facString=facString[:-1]
		"""
        sstring = selfcopy.affstring()
        if len(facString) > 0:
            if ('+' in sstring) or ('-' in sstring):
                sstring = '( ' + sstring + ' )'
            if ('+' in facString) or ('-' in facString):
                facString = '( ' + facString + ' )'

            selfcopy.string = facString + '∘' + sstring

        return selfcopy 

def __gt__(self, exp):

        if isinstance(exp, AffinExp):
            if exp.size != (1, 1):
                raise Exception(
                    'lower bound of a sum_k_smallest must be scalar')
            from .problem import Problem
            Ptmp = Problem()
            if self.eigenvalues:
                n = self.exp.size[0]
                I = new_param('I', cvx.spdiag([1.] * n))
                if self.k == n:
                    return (I | self.exp) < exp
                elif self.k == 1:
                    cons = self.exp >> exp * I
                    cons.myconstring = self.string + '>=' + exp.string
                    return cons
                else:
                    s = Ptmp.add_variable('s', 1)
                    Z = Ptmp.add_variable('Z', (n, n), 'symmetric')
                    Ptmp.add_constraint(Z >> 0)
                    Ptmp.add_constraint(-self.exp << Z + s * I)
                    Ptmp.add_constraint(-exp > (I | Z) + (self.k * s))
            else:
                n = self.exp.size[0] * self.exp.size[1]
                if self.k == 1:
                    cons = self.exp > exp
                    cons.myconstring = self.string + '>=' + exp.string
                    return cons
                elif self.k == n:
                    return (1 | self.exp) > exp
                else:
                    lbda = Ptmp.add_variable('lambda', 1)
                    mu = Ptmp.add_variable('mu', self.exp.size, lower=0)
                    Ptmp.add_constraint(-self.exp < lbda + mu)
                    Ptmp.add_constraint(self.k * lbda + (1 | mu) < -exp)

            return Sumklargest_Constraint(
                exp,
                self.exp,
                self.k,
                self.eigenvalues,
                False,
                Ptmp,
                self.string +
                '>' +
                exp.string)

        else:  # constant
            term, termString = _retrieve_matrix(exp, (1, 1))
            exp1 = AffinExp(
                factors={}, constant=term, size=(
                    1, 1), string=termString)
            return self > exp1