""" Convex programming solver. A primal-dual interior-point solver written in Python and interfaces for quadratic and geometric programming. Also includes an interface to the quadratic programming solver from MOSEK. """ # Copyright 2012-2023 M. Andersen and L. Vandenberghe. # Copyright 2010-2011 L. Vandenberghe. # Copyright 2004-2009 J. Dahl and L. Vandenberghe. # # This file is part of CVXOPT. # # CVXOPT is free software; you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation; either version 3 of the License, or # (at your option) any later version. # # CVXOPT is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program. If not, see . import sys if sys.version > '3': long = int __all__ = [] options = {} def cpl(c, F, G = None, h = None, dims = None, A = None, b = None, kktsolver = None, xnewcopy = None, xdot = None, xaxpy = None, xscal = None, ynewcopy = None, ydot = None, yaxpy = None, yscal = None, **kwargs): """ Solves a convex optimization problem with a linear objective minimize c'*x subject to f(x) <= 0 G*x <= h A*x = b. f is vector valued, convex and twice differentiable. The linear inequalities are with respect to a cone C defined as the Cartesian product of N + M + 1 cones: C = C_0 x C_1 x .... x C_N x C_{N+1} x ... x C_{N+M}. The first cone C_0 is the nonnegative orthant of dimension ml. The next N cones are second order cones of dimension mq[0], ..., mq[N-1]. The second order cone of dimension m is defined as { (u0, u1) in R x R^{m-1} | u0 >= ||u1||_2 }. The next M cones are positive semidefinite cones of order ms[0], ..., ms[M-1] >= 0. Input arguments (basic usage). c is a dense 'd' matrix of size (n,1). F is a function that handles the following arguments. F() returns a tuple (mnl, x0). mnl is the number of nonlinear inequality constraints. x0 is a point in the domain of f. F(x) returns a tuple (f, Df). f is a dense 'd' matrix of size (mnl, 1) containing f(x). Df is a dense or sparse 'd' matrix of size (mnl, n), containing the derivatives of f at x: Df[k,:] is the transpose of the gradient of fk at x. If x is not in dom f, F(x) returns None or (None, None). F(x, z) with z a positive 'd' matrix of size (mnl,1), returns a tuple (f, Df, H). f and Df are defined as above. H is a dense or sparse 'd' matrix of size (n,n). The lower triangular part of H contains the lower triangular part of sum_k z[k] * Hk where Hk is the Hessian of fk at x. If F is called with two arguments, it can be assumed that x is dom f. If Df and H are returned as sparse matrices, their sparsity patterns must be the same for each call to F(x) or F(x,z). dims is a dictionary with the dimensions of the components of C. It has three fields. - dims['l'] = ml, the dimension of the nonnegative orthant C_0. (ml >= 0.) - dims['q'] = mq = [ mq[0], mq[1], ..., mq[N-1] ], a list of N integers with the dimensions of the second order cones C_1, ..., C_N. (N >= 0 and mq[k] >= 1.) - dims['s'] = ms = [ ms[0], ms[1], ..., ms[M-1] ], a list of M integers with the orders of the semidefinite cones C_{N+1}, ..., C_{N+M}. (M >= 0 and ms[k] >= 0.) The default value of dims is {'l': G.size[0], 'q': [], 's': []}. G is a dense or sparse 'd' matrix of size (K,n), where K = ml + mq[0] + ... + mq[N-1] + ms[0]**2 + ... + ms[M-1]**2. Each column of G describes a vector v = ( v_0, v_1, ..., v_N, vec(v_{N+1}), ..., vec(v_{N+M}) ) in V = R^ml x R^mq[0] x ... x R^mq[N-1] x S^ms[0] x ... x S^ms[M-1] stored as a column vector [ v_0; v_1; ...; v_N; vec(v_{N+1}); ...; vec(v_{N+M}) ]. Here, if u is a symmetric matrix of order m, then vec(u) is the matrix u stored in column major order as a vector of length m**2. We use BLAS unpacked 'L' storage, i.e., the entries in vec(u) corresponding to the strictly upper triangular entries of u are not referenced. h is a dense 'd' matrix of size (K,1), representing a vector in V, in the same format as the columns of G. A is a dense or sparse 'd' matrix of size (p,n). The default value is a sparse 'd' matrix of size (0,n). b is a dense 'd' matrix of size (p,1). The default value is a dense 'd' matrix of size (0,1). It is assumed that rank(A) = p and rank([H; A; Df; G]) = n at all x in dom f. The other arguments are normally not needed. They make it possible to exploit certain types of structure, as described further below. Output arguments. Returns a dictionary with keys 'status', 'x', 'snl', 'sl', 'znl', 'zl', 'y', 'primal objective', 'dual objective', 'gap', 'relative gap', 'primal infeasibility', 'dual infeasibility', 'primal slack', 'dual slack'. The 'status' field has values 'optimal' or 'unknown'. If status is 'optimal', x, snl, sl, y, znl, zl are an approximate solution of the primal and dual optimality conditions f(x) + snl = 0, G*x + sl = h, A*x = b Df(x)'*znl + G'*zl + A'*y + c = 0 snl >= 0, znl >= 0, sl >= 0, zl >= 0 snl'*znl + sl'* zl = 0. If status is 'unknown', x, snl, sl, y, znl, zl are the last iterates before termination. They satisfy snl > 0, znl > 0, sl > 0, zl > 0, but are not necessarily feasible. The values of the other fields are defined as follows. - 'primal objective': the primal objective c'*x. - 'dual objective': the dual objective L(x,y,znl,zl) = c'*x + znl'*f(x) + zl'*(G*x-h) + y'*(A*x-b). - 'gap': the duality gap snl'*znl + sl'*zl. - 'relative gap': the relative gap, defined as gap / -primal objective if the primal objective is negative, gap / dual objective if the dual objective is positive, and None otherwise. - 'primal infeasibility': the residual in the primal constraints, defined as || (f(x) + snl, G*x + sl - h, A*x-b) ||_2 divided by max(1, || (f(x0) + 1, G*x0 + 1 - h, A*x0 - b) ||_2 ) where x0 is the point returned by F(). - 'dual infeasibility': the residual in the dual constraints, defined as || c + Df(x)'*znl + G'*zl + A'*y ||_2 divided by max(1, || c + Df(x0)'*1 + G'*1 ||_2 ). - 'primal slack': the smallest primal slack, min( min_k sl_k, sup {t | sl >= te} ) where e = ( e_0, e_1, ..., e_N, e_{N+1}, ..., e_{M+N} ) is the identity vector in C. e_0 is an ml-vector of ones, e_k, k = 1,..., N, is the unit vector (1,0,...,0) of length mq[k], and e_k = vec(I) where I is the identity matrix of order ms[k]. - 'dual slack': the smallest dual slack, min( min_k zl_k, sup {t | zl >= te} ). If the exit status is 'optimal', then the primal and dual infeasibilities are guaranteed to be less than solvers.options['feastol'] (default 1e-7). The gap is less than solvers.options['abstol'] (default 1e-7) or the relative gap is less than solvers.options['reltol'] (defaults 1e-6). Termination with status 'unknown' indicates that the algorithm failed to find a solution that satisfies the specified tolerances. In some cases, the returned solution may be fairly accurate. If the primal and dual infeasibilities, the gap, and the relative gap are small, then x, y, snl, sl, znl, zl are close to optimal. Advanced usage. Three mechanisms are provided to express problem structure. 1. The user can provide a customized routine for solving linear equations ('KKT systems') [ sum_k zk*Hk(x) A' GG' ] [ ux ] [ bx ] [ A 0 0 ] [ uy ] = [ by ] [ GG 0 -W'*W ] [ uz ] [ bz ] where GG = [ Df(x); G ], uz = (uznl, uzl), bz = (bznl, bzl). z is a positive vector of length mnl and x is a point in the domain of f. W is a scaling matrix, i.e., a block diagonal mapping W*u = ( Wnl*unl, W0*u_0, ..., W_{N+M}*u_{N+M} ) defined as follows. - For the nonlinear block (Wnl): Wnl = diag(dnl), with dnl a positive vector of length mnl. - For the 'l' block (W_0): W_0 = diag(d), with d a positive vector of length ml. - For the 'q' blocks (W_{k+1}, k = 0, ..., N-1): W_{k+1} = beta_k * ( 2 * v_k * v_k' - J ) where beta_k is a positive scalar, v_k is a vector in R^mq[k] with v_k[0] > 0 and v_k'*J*v_k = 1, and J = [1, 0; 0, -I]. - For the 's' blocks (W_{k+N}, k = 0, ..., M-1): W_k * u = vec(r_k' * mat(u) * r_k) where r_k is a nonsingular matrix of order ms[k], and mat(x) is the inverse of the vec operation. The optional argument kktsolver is a Python function that will be called as g = kktsolver(x, z, W). W is a dictionary that contains the parameters of the scaling: - W['dnl'] is a positive 'd' matrix of size (mnl, 1). - W['dnli'] is a positive 'd' matrix with the elementwise inverse of W['dnl']. - W['d'] is a positive 'd' matrix of size (ml, 1). - W['di'] is a positive 'd' matrix with the elementwise inverse of W['d']. - W['beta'] is a list [ beta_0, ..., beta_{N-1} ] - W['v'] is a list [ v_0, ..., v_{N-1} ] - W['r'] is a list [ r_0, ..., r_{M-1} ] - W['rti'] is a list [ rti_0, ..., rti_{M-1} ], with rti_k the inverse of the transpose of r_k. The call g = kktsolver(x, z, W) should return a function g that solves the KKT system by g(x, y, z). On entry, x, y, z contain the righthand side bx, by, bz. On exit, they contain the solution, with uz scaled: W*uz is returned instead of uz. In other words, on exit x, y, z are the solution of [ sum_k zk*Hk(x) A' GG'*W^{-1} ] [ ux ] [ bx ] [ A 0 0 ] [ uy ] = [ by ]. [ GG 0 -W' ] [ uz ] [ bz ] 2. The linear operators Df*u, H*u, G*u and A*u can be specified by providing Python functions instead of matrices. This can only be done in combination with 1. above, i.e., it also requires the kktsolver argument. If G is a function, the call G(u, v, alpha, beta, trans) should evaluate the matrix-vector products v := alpha * G * u + beta * v if trans is 'N' v := alpha * G' * u + beta * v if trans is 'T'. The arguments u and v are required. The other arguments have default values alpha = 1.0, beta = 0.0, trans = 'N'. If A is a function, the call A(u, v, alpha, beta, trans) should evaluate the matrix-vectors products v := alpha * A * u + beta * v if trans is 'N' v := alpha * A' * u + beta * v if trans is 'T'. The arguments u and v are required. The other arguments have default values alpha = 1.0, beta = 0.0, trans = 'N'. If Df is a function, the call Df(u, v, alpha, beta, trans) should evaluate the matrix-vectors products v := alpha * Df(x) * u + beta * v if trans is 'N' v := alpha * Df(x)' * u + beta * v if trans is 'T'. If H is a function, the call H(u, v, alpha, beta) should evaluate the matrix-vectors product v := alpha * H * u + beta * v. 3. Instead of using the default representation of the primal variable x and the dual variable y as one-column 'd' matrices, we can represent these variables and the corresponding parameters c and b by arbitrary Python objects (matrices, lists, dictionaries, etc). This can only be done in combination with 1. and 2. above, i.e., it requires a user-provided KKT solver and a function description of the linear mappings. It also requires the arguments xnewcopy, xdot, xscal, xaxpy, ynewcopy, ydot, yscal, yaxpy. These arguments are functions defined as follows. If X is the vector space of primal variables x, then: - xnewcopy(u) creates a new copy of the vector u in X. - xdot(u, v) returns the inner product of two vectors u and v in X. - xscal(alpha, u) computes u := alpha*u, where alpha is a scalar and u is a vector in X. - xaxpy(u, v, alpha = 1.0, beta = 0.0) computes v := alpha*u + v for a scalar alpha and two vectors u and v in X. If Y is the vector space of primal variables y: - ynewcopy(u) creates a new copy of the vector u in Y. - ydot(u, v) returns the inner product of two vectors u and v in Y. - yscal(alpha, u) computes u := alpha*u, where alpha is a scalar and u is a vector in Y. - yaxpy(u, v, alpha = 1.0, beta = 0.0) computes v := alpha*u + v for a scalar alpha and two vectors u and v in Y. Control parameters. The following control parameters can be modified by adding an entry to the dictionary options. options['show_progress'] True/False (default: True) options['maxiters'] positive integer (default: 100) options['refinement'] nonnegative integer (default: 1) options['abstol'] scalar (default: 1e-7) options['reltol'] scalar (default: 1e-6) options['feastol'] scalar (default: 1e-7). """ import math from cvxopt import base, blas, misc from cvxopt.base import matrix, spmatrix STEP = 0.99 BETA = 0.5 ALPHA = 0.01 EXPON = 3 MAX_RELAXED_ITERS = 8 options = kwargs.get('options',globals()['options']) DEBUG = options.get('debug',False) KKTREG = options.get('kktreg',None) if KKTREG is None: pass elif not isinstance(KKTREG,(float,int,long)) or KKTREG < 0.0: raise ValueError("options['kktreg'] must be a nonnegative scalar") MAXITERS = options.get('maxiters',100) if not isinstance(MAXITERS,(int,long)) or MAXITERS < 1: raise ValueError("options['maxiters'] must be a positive integer") ABSTOL = options.get('abstol',1e-7) if not isinstance(ABSTOL,(float,int,long)): raise ValueError("options['abstol'] must be a scalar") RELTOL = options.get('reltol',1e-6) if not isinstance(RELTOL,(float,int,long)): raise ValueError("options['reltol'] must be a scalar") if RELTOL <= 0.0 and ABSTOL <= 0.0 : raise ValueError("at least one of options['reltol'] and " \ "options['abstol'] must be positive") FEASTOL = options.get('feastol',1e-7) if not isinstance(FEASTOL,(float,int,long)) or FEASTOL <= 0.0: raise ValueError("options['feastol'] must be a positive scalar") show_progress = options.get('show_progress', True) refinement = options.get('refinement',1) if not isinstance(refinement,(int,long)) or refinement < 0: raise ValueError("options['refinement'] must be a nonnegative integer") if kktsolver is None: if dims and (dims['q'] or dims['s']): kktsolver = 'chol' else: kktsolver = 'chol2' defaultsolvers = ('ldl', 'ldl2', 'chol', 'chol2') if type(kktsolver) is str and kktsolver not in defaultsolvers: raise ValueError("'%s' is not a valid value for kktsolver" \ %kktsolver) try: mnl, x0 = F() except: raise ValueError("function call 'F()' failed") # Argument error checking depends on level of customization. customkkt = type(kktsolver) is not str operatorG = G is not None and type(G) not in (matrix, spmatrix) operatorA = A is not None and type(A) not in (matrix, spmatrix) if (operatorG or operatorA) and not customkkt: raise ValueError("use of function valued G, A requires a "\ "user-provided kktsolver") customx = (xnewcopy != None or xdot != None or xaxpy != None or xscal != None) if customx and (not operatorG or not operatorA or not customkkt): raise ValueError("use of non-vector type for x requires "\ "function valued G, A and user-provided kktsolver") customy = (ynewcopy != None or ydot != None or yaxpy != None or yscal != None) if customy and (not operatorA or not customkkt): raise ValueError("use of non vector type for y requires "\ "function valued A and user-provided kktsolver") if not customx: if type(x0) is not matrix or x0.typecode != 'd' or x0.size[1] != 1: raise TypeError("'x0' must be a 'd' matrix with one column") if type(c) is not matrix or c.typecode != 'd' or c.size != x0.size: raise TypeError("'c' must be a 'd' matrix of size (%d,%d)"\ %(x0.size[0],1)) if h is None: h = matrix(0.0, (0,1)) if type(h) is not matrix or h.typecode != 'd' or h.size[1] != 1: raise TypeError("'h' must be a 'd' matrix with 1 column") if not dims: dims = {'l': h.size[0], 'q': [], 's': []} # Dimension of the product cone of the linear inequalities. with 's' # components unpacked. cdim = dims['l'] + sum(dims['q']) + sum([ k**2 for k in dims['s'] ]) if h.size[0] != cdim: raise TypeError("'h' must be a 'd' matrix of size (%d,1)" %cdim) if G is None: if customx: def G(x, y, trans = 'N', alpha = 1.0, beta = 0.0): if trans == 'N': pass else: xscal(beta, y) else: G = spmatrix([], [], [], (0, c.size[0])) if not operatorG: if G.typecode != 'd' or G.size != (cdim, c.size[0]): raise TypeError("'G' must be a 'd' matrix with size (%d, %d)"\ %(cdim, c.size[0])) def fG(x, y, trans = 'N', alpha = 1.0, beta = 0.0): misc.sgemv(G, x, y, dims, trans = trans, alpha = alpha, beta = beta) else: fG = G if A is None: if customx or customy: def A(x, y, trans = 'N', alpha = 1.0, beta = 0.0): if trans == 'N': pass else: yscal(beta, y) else: A = spmatrix([], [], [], (0, c.size[0])) if not operatorA: if A.typecode != 'd' or A.size[1] != c.size[0]: raise TypeError("'A' must be a 'd' matrix with %d columns" \ %c.size[0]) def fA(x, y, trans = 'N', alpha = 1.0, beta = 0.0): base.gemv(A, x, y, trans = trans, alpha = alpha, beta = beta) else: fA = A if not customy: if b is None: b = matrix(0.0, (0,1)) if type(b) is not matrix or b.typecode != 'd' or b.size[1] != 1: raise TypeError("'b' must be a 'd' matrix with one column") if not operatorA and b.size[0] != A.size[0]: raise TypeError("'b' must have length %d" %A.size[0]) if b is None and customy: raise ValueEror("use of non vector type for y requires b") # kktsolver(x, z, W) returns a routine for solving # # [ sum_k zk*Hk(x) A' GG'*W^{-1} ] [ ux ] [ bx ] # [ A 0 0 ] [ uy ] = [ by ] # [ GG 0 -W' ] [ uz ] [ bz ] # # where G = [Df(x); G]. if kktsolver in defaultsolvers: if kktsolver == 'ldl': factor = misc.kkt_ldl(G, dims, A, mnl, kktreg = KKTREG) elif kktsolver == 'ldl2': factor = misc.kkt_ldl2(G, dims, A, mnl) elif kktsolver == 'chol': factor = misc.kkt_chol(G, dims, A, mnl) else: factor = misc.kkt_chol2(G, dims, A, mnl) def kktsolver(x, z, W): f, Df, H = F(x, z) return factor(W, H, Df) if xnewcopy is None: xnewcopy = matrix if xdot is None: xdot = blas.dot if xaxpy is None: xaxpy = blas.axpy if xscal is None: xscal = blas.scal def xcopy(x, y): xscal(0.0, y) xaxpy(x, y) if ynewcopy is None: ynewcopy = matrix if ydot is None: ydot = blas.dot if yaxpy is None: yaxpy = blas.axpy if yscal is None: yscal = blas.scal def ycopy(x, y): yscal(0.0, y) yaxpy(x, y) # Initial points x = xnewcopy(x0) y = ynewcopy(b); yscal(0.0, y) z, s = matrix(0.0, (mnl + cdim, 1)), matrix(0.0, (mnl + cdim, 1)) z[: mnl+dims['l']] = 1.0 s[: mnl+dims['l']] = 1.0 ind = mnl + dims['l'] for m in dims['q']: z[ind] = 1.0 s[ind] = 1.0 ind += m for m in dims['s']: z[ind : ind + m*m : m+1] = 1.0 s[ind : ind + m*m : m+1] = 1.0 ind += m**2 rx, ry = xnewcopy(x0), ynewcopy(b) rznl, rzl = matrix(0.0, (mnl, 1)), matrix(0.0, (cdim, 1)), dx, dy = xnewcopy(x), ynewcopy(y) dz, ds = matrix(0.0, (mnl + cdim, 1)), matrix(0.0, (mnl + cdim, 1)) lmbda = matrix(0.0, (mnl + dims['l'] + sum(dims['q']) + sum(dims['s']), 1)) lmbdasq = matrix(0.0, (mnl + dims['l'] + sum(dims['q']) + sum(dims['s']), 1)) sigs = matrix(0.0, (sum(dims['s']), 1)) sigz = matrix(0.0, (sum(dims['s']), 1)) dz2, ds2 = matrix(0.0, (mnl + cdim, 1)), matrix(0.0, (mnl + cdim, 1)) newx, newy = xnewcopy(x), ynewcopy(y) newz, news = matrix(0.0, (mnl + cdim, 1)), matrix(0.0, (mnl + cdim, 1)) newrx = xnewcopy(x0) newrznl = matrix(0.0, (mnl, 1)) rx0, ry0 = xnewcopy(x0), ynewcopy(b) rznl0, rzl0 = matrix(0.0, (mnl, 1)), matrix(0.0, (cdim, 1)), x0, dx0 = xnewcopy(x), xnewcopy(dx) y0, dy0 = ynewcopy(y), ynewcopy(dy) z0 = matrix(0.0, (mnl + cdim, 1)) dz0 = matrix(0.0, (mnl + cdim, 1)) dz20 = matrix(0.0, (mnl + cdim, 1)) s0 = matrix(0.0, (mnl + cdim, 1)) ds0 = matrix(0.0, (mnl + cdim, 1)) ds20 = matrix(0.0, (mnl + cdim, 1)) W0 = {} W0['dnl'] = matrix(0.0, (mnl, 1)) W0['dnli'] = matrix(0.0, (mnl, 1)) W0['d'] = matrix(0.0, (dims['l'], 1)) W0['di'] = matrix(0.0, (dims['l'], 1)) W0['v'] = [ matrix(0.0, (m, 1)) for m in dims['q'] ] W0['beta'] = len(dims['q']) * [ 0.0 ] W0['r'] = [ matrix(0.0, (m, m)) for m in dims['s'] ] W0['rti'] = [ matrix(0.0, (m, m)) for m in dims['s'] ] lmbda0 = matrix(0.0, (mnl + dims['l'] + sum(dims['q']) + sum(dims['s']), 1)) lmbdasq0 = matrix(0.0, (mnl + dims['l'] + sum(dims['q']) + sum(dims['s']), 1)) if show_progress: print("% 10s% 12s% 10s% 8s% 7s" %("pcost", "dcost", "gap", "pres", "dres")) relaxed_iters = 0 for iters in range(MAXITERS + 1): if refinement or DEBUG: # We need H to compute residuals of KKT equations. f, Df, H = F(x, z[:mnl]) else: f, Df = F(x) f = matrix(f, tc='d') if f.typecode != 'd' or f.size != (mnl, 1): raise TypeError("first output argument of F() must be a "\ "'d' matrix of size (%d, %d)" %(mnl, 1)) if type(Df) is matrix or type(Df) is spmatrix: if customx: raise ValueError("use of non-vector type for x "\ "requires function valued Df") if Df.typecode != 'd' or Df.size != (mnl, c.size[0]): raise TypeError("second output argument of F() must "\ "be a 'd' matrix of size (%d,%d)" %(mnl, c.size[0])) def fDf(u, v, alpha = 1.0, beta = 0.0, trans = 'N'): base.gemv(Df, u, v, alpha = alpha, beta = beta, trans = trans) else: if not customkkt: raise ValueError("use of function valued Df requires "\ "a user-provided kktsolver") fDf = Df if refinement or DEBUG: if type(H) is matrix or type(H) is spmatrix: if customx: raise ValueError("use of non-vector type "\ "for x requires function valued H") if H.typecode != 'd' or H.size != (c.size[0], c.size[0]): raise TypeError("third output argument of F() must "\ "be a 'd' matrix of size (%d,%d)" \ %(c.size[0], c.size[0])) def fH(u, v, alpha = 1.0, beta = 0.0): base.symv(H, u, v, alpha = alpha, beta = beta) else: if not customkkt: raise ValueError("use of function valued H requires "\ "a user-provided kktsolver") fH = H gap = misc.sdot(s, z, dims, mnl) # rx = c + A'*y + Df'*z[:mnl] + G'*z[mnl:] xcopy(c, rx) fA(y, rx, beta = 1.0, trans = 'T') fDf(z[:mnl], rx, beta = 1.0, trans = 'T') fG(z[mnl:], rx, beta = 1.0, trans = 'T') resx = math.sqrt(xdot(rx, rx)) # ry = A*x - b ycopy(b, ry) fA(x, ry, alpha = 1.0, beta = -1.0) resy = math.sqrt(ydot(ry, ry)) # rznl = s[:mnl] + f blas.copy(s[:mnl], rznl) blas.axpy(f, rznl) resznl = blas.nrm2(rznl) # rzl = s[mnl:] + G*x - h blas.copy(s[mnl:], rzl) blas.axpy(h, rzl, alpha = -1.0) fG(x, rzl, beta = 1.0) reszl = misc.snrm2(rzl, dims) # Statistics for stopping criteria. # pcost = c'*x # dcost = c'*x + y'*(A*x-b) + znl'*f(x) + zl'*(G*x-h) # = c'*x + y'*(A*x-b) + znl'*(f(x)+snl) + zl'*(G*x-h+sl) # - z'*s # = c'*x + y'*ry + znl'*rznl + zl'*rzl - gap pcost = xdot(c,x) dcost = pcost + ydot(y, ry) + blas.dot(z[:mnl], rznl) + \ misc.sdot(z[mnl:], rzl, dims) - gap if pcost < 0.0: relgap = gap / -pcost elif dcost > 0.0: relgap = gap / dcost else: relgap = None pres = math.sqrt( resy**2 + resznl**2 + reszl**2 ) dres = resx if iters == 0: resx0 = max(1.0, resx) resznl0 = max(1.0, resznl) pres0 = max(1.0, pres) dres0 = max(1.0, dres) gap0 = gap theta1 = 1.0 / gap0 theta2 = 1.0 / resx0 theta3 = 1.0 / resznl0 phi = theta1 * gap + theta2 * resx + theta3 * resznl pres = pres / pres0 dres = dres / dres0 if show_progress: print("%2d: % 8.4e % 8.4e % 4.0e% 7.0e% 7.0e" \ %(iters, pcost, dcost, gap, pres, dres)) # Stopping criteria. if ( pres <= FEASTOL and dres <= FEASTOL and ( gap <= ABSTOL or (relgap is not None and relgap <= RELTOL) )) or \ iters == MAXITERS: sl, zl = s[mnl:], z[mnl:] ind = dims['l'] + sum(dims['q']) for m in dims['s']: misc.symm(sl, m, ind) misc.symm(zl, m, ind) ind += m**2 ts = misc.max_step(s, dims, mnl) tz = misc.max_step(z, dims, mnl) if iters == MAXITERS: if show_progress: print("Terminated (maximum number of iterations "\ "reached).") status = 'unknown' else: if show_progress: print("Optimal solution found.") status = 'optimal' return {'status': status, 'x': x, 'y': y, 'znl': z[:mnl], 'zl': zl, 'snl': s[:mnl], 'sl': sl, 'gap': gap, 'relative gap': relgap, 'primal objective': pcost, 'dual objective': dcost, 'primal slack': -ts, 'dual slack': -tz, 'primal infeasibility': pres, 'dual infeasibility': dres } # Compute initial scaling W: # # W * z = W^{-T} * s = lambda. # # lmbdasq = lambda o lambda if iters == 0: W = misc.compute_scaling(s, z, lmbda, dims, mnl) misc.ssqr(lmbdasq, lmbda, dims, mnl) # f3(x, y, z) solves # # [ H A' GG'*W^{-1} ] [ ux ] [ bx ] # [ A 0 0 ] [ uy ] = [ by ]. # [ GG 0 -W' ] [ uz ] [ bz ] # # On entry, x, y, z contain bx, by, bz. # On exit, they contain ux, uy, uz. try: f3 = kktsolver(x, z[:mnl], W) except ArithmeticError: singular_kkt_matrix = False if iters == 0: raise ValueError("Rank(A) < p or "\ "Rank([H(x); A; Df(x); G]) < n") elif 0 < relaxed_iters < MAX_RELAXED_ITERS > 0: # The arithmetic error may be caused by a relaxed line # search in the previous iteration. Therefore we restore # the last saved state and require a standard line search. phi, gap = phi0, gap0 mu = gap / ( mnl + dims['l'] + len(dims['q']) + sum(dims['s']) ) blas.copy(W0['dnl'], W['dnl']) blas.copy(W0['dnli'], W['dnli']) blas.copy(W0['d'], W['d']) blas.copy(W0['di'], W['di']) for k in range(len(dims['q'])): blas.copy(W0['v'][k], W['v'][k]) W['beta'][k] = W0['beta'][k] for k in range(len(dims['s'])): blas.copy(W0['r'][k], W['r'][k]) blas.copy(W0['rti'][k], W['rti'][k]) xcopy(x0, x); ycopy(y0, y); blas.copy(s0, s); blas.copy(z0, z) blas.copy(lmbda0, lmbda) blas.copy(lmbdasq, lmbdasq0) xcopy(rx0, rx); ycopy(ry0, ry) resx = math.sqrt(xdot(rx, rx)) blas.copy(rznl0, rznl); blas.copy(rzl0, rzl); resznl = blas.nrm2(rznl) relaxed_iters = -1 try: f3 = kktsolver(x, z[:mnl], W) except ArithmeticError: singular_kkt_matrix = True else: singular_kkt_matrix = True if singular_kkt_matrix: sl, zl = s[mnl:], z[mnl:] ind = dims['l'] + sum(dims['q']) for m in dims['s']: misc.symm(sl, m, ind) misc.symm(zl, m, ind) ind += m**2 ts = misc.max_step(s, dims, mnl) tz = misc.max_step(z, dims, mnl) if show_progress: print("Terminated (singular KKT matrix).") status = 'unknown' return {'status': status, 'x': x, 'y': y, 'znl': z[:mnl], 'zl': zl, 'snl': s[:mnl], 'sl': sl, 'gap': gap, 'relative gap': relgap, 'primal objective': pcost, 'dual objective': dcost, 'primal infeasibility': pres, 'dual infeasibility': dres, 'primal slack': -ts, 'dual slack': -tz } # f4_no_ir(x, y, z, s) solves # # [ 0 ] [ H A' GG' ] [ ux ] [ bx ] # [ 0 ] + [ A 0 0 ] [ uy ] = [ by ] # [ W'*us ] [ GG 0 0 ] [ W^{-1}*uz ] [ bz ] # # lmbda o (uz + us) = bs. # # On entry, x, y, z, x, contain bx, by, bz, bs. # On exit, they contain ux, uy, uz, us. if iters == 0: ws3 = matrix(0.0, (mnl + cdim, 1)) wz3 = matrix(0.0, (mnl + cdim, 1)) def f4_no_ir(x, y, z, s): # Solve # # [ H A' GG' ] [ ux ] [ bx ] # [ A 0 0 ] [ uy ] = [ by ] # [ GG 0 -W'*W ] [ W^{-1}*uz ] [ bz - W'*(lmbda o\ bs) ] # # us = lmbda o\ bs - uz. # s := lmbda o\ s # = lmbda o\ bs misc.sinv(s, lmbda, dims, mnl) # z := z - W'*s # = bz - W' * (lambda o\ bs) blas.copy(s, ws3) misc.scale(ws3, W, trans = 'T') blas.axpy(ws3, z, alpha = -1.0) # Solve for ux, uy, uz f3(x, y, z) # s := s - z # = lambda o\ bs - z. blas.axpy(z, s, alpha = -1.0) if iters == 0: wz2nl, wz2l = matrix(0.0, (mnl,1)), matrix(0.0, (cdim, 1)) def res(ux, uy, uz, us, vx, vy, vz, vs): # Evaluates residuals in Newton equations: # # [ vx ] [ 0 ] [ H A' GG' ] [ ux ] # [ vy ] -= [ 0 ] + [ A 0 0 ] [ uy ] # [ vz ] [ W'*us ] [ GG 0 0 ] [ W^{-1}*uz ] # # vs -= lmbda o (uz + us). # vx := vx - H*ux - A'*uy - GG'*W^{-1}*uz fH(ux, vx, alpha = -1.0, beta = 1.0) fA(uy, vx, alpha = -1.0, beta = 1.0, trans = 'T') blas.copy(uz, wz3) misc.scale(wz3, W, inverse = 'I') fDf(wz3[:mnl], vx, alpha = -1.0, beta = 1.0, trans = 'T') fG(wz3[mnl:], vx, alpha = -1.0, beta = 1.0, trans = 'T') # vy := vy - A*ux fA(ux, vy, alpha = -1.0, beta = 1.0) # vz := vz - W'*us - GG*ux fDf(ux, wz2nl) blas.axpy(wz2nl, vz, alpha = -1.0) fG(ux, wz2l) blas.axpy(wz2l, vz, alpha = -1.0, offsety = mnl) blas.copy(us, ws3) misc.scale(ws3, W, trans = 'T') blas.axpy(ws3, vz, alpha = -1.0) # vs -= lmbda o (uz + us) blas.copy(us, ws3) blas.axpy(uz, ws3) misc.sprod(ws3, lmbda, dims, mnl, diag = 'D') blas.axpy(ws3, vs, alpha = -1.0) # f4(x, y, z, s) solves the same system as f4_no_ir, but applies # iterative refinement. if iters == 0: if refinement or DEBUG: wx, wy = xnewcopy(c), ynewcopy(b) wz = matrix(0.0, (mnl + cdim, 1)) ws = matrix(0.0, (mnl + cdim, 1)) if refinement: wx2, wy2 = xnewcopy(c), ynewcopy(b) wz2 = matrix(0.0, (mnl + cdim, 1)) ws2 = matrix(0.0, (mnl + cdim, 1)) def f4(x, y, z, s): if refinement or DEBUG: xcopy(x, wx) ycopy(y, wy) blas.copy(z, wz) blas.copy(s, ws) f4_no_ir(x, y, z, s) for i in range(refinement): xcopy(wx, wx2) ycopy(wy, wy2) blas.copy(wz, wz2) blas.copy(ws, ws2) res(x, y, z, s, wx2, wy2, wz2, ws2) f4_no_ir(wx2, wy2, wz2, ws2) xaxpy(wx2, x) yaxpy(wy2, y) blas.axpy(wz2, z) blas.axpy(ws2, s) if DEBUG: res(x, y, z, s, wx, wy, wz, ws) print("KKT residuals:") print(" 'x': %e" %math.sqrt(xdot(wx, wx))) print(" 'y': %e" %math.sqrt(ydot(wy, wy))) print(" 'z': %e" %misc.snrm2(wz, dims, mnl)) print(" 's': %e" %misc.snrm2(ws, dims, mnl)) sigma, eta = 0.0, 0.0 for i in [0, 1]: # Solve # # [ 0 ] [ H A' GG' ] [ dx ] # [ 0 ] + [ A 0 0 ] [ dy ] = -(1 - eta)*r # [ W'*ds ] [ GG 0 0 ] [ W^{-1}*dz ] # # lmbda o (dz + ds) = -lmbda o lmbda + sigma*mu*e. # mu = gap / (mnl + dims['l'] + len(dims['q']) + sum(dims['s'])) # ds = -lmbdasq + sigma * mu * e blas.scal(0.0, ds) blas.axpy(lmbdasq, ds, n = mnl + dims['l'] + sum(dims['q']), alpha = -1.0) ds[:mnl + dims['l']] += sigma*mu ind = mnl + dims['l'] for m in dims['q']: ds[ind] += sigma*mu ind += m ind2 = ind for m in dims['s']: blas.axpy(lmbdasq, ds, n = m, offsetx = ind2, offsety = ind, incy = m + 1, alpha = -1.0) ds[ind : ind + m*m : m+1] += sigma*mu ind += m*m ind2 += m # (dx, dy, dz) := -(1-eta) * (rx, ry, rz) xscal(0.0, dx); xaxpy(rx, dx, alpha = -1.0 + eta) yscal(0.0, dy); yaxpy(ry, dy, alpha = -1.0 + eta) blas.scal(0.0, dz) blas.axpy(rznl, dz, alpha = -1.0 + eta) blas.axpy(rzl, dz, alpha = -1.0 + eta, offsety = mnl) try: f4(dx, dy, dz, ds) except ArithmeticError: if iters == 0: raise ValueError("Rank(A) < p or "\ "Rank([H(x); A; Df(x); G]) < n") else: sl, zl = s[mnl:], z[mnl:] ind = dims['l'] + sum(dims['q']) for m in dims['s']: misc.symm(sl, m, ind) misc.symm(zl, m, ind) ind += m**2 ts = misc.max_step(s, dims, mnl) tz = misc.max_step(z, dims, mnl) if show_progress: print("Terminated (singular KKT matrix).") return {'status': 'unknown', 'x': x, 'y': y, 'znl': z[:mnl], 'zl': zl, 'snl': s[:mnl], 'sl': sl, 'gap': gap, 'relative gap': relgap, 'primal objective': pcost, 'dual objective': dcost, 'primal infeasibility': pres, 'dual infeasibility': dres, 'primal slack': -ts, 'dual slack': -tz } # Inner product ds'*dz and unscaled steps are needed in the # line search. dsdz = misc.sdot(ds, dz, dims, mnl) blas.copy(dz, dz2) misc.scale(dz2, W, inverse = 'I') blas.copy(ds, ds2) misc.scale(ds2, W, trans = 'T') # Maximum steps to boundary. # # Also compute the eigenvalue decomposition of 's' blocks in # ds, dz. The eigenvectors Qs, Qz are stored in ds, dz. # The eigenvalues are stored in sigs, sigz. misc.scale2(lmbda, ds, dims, mnl) ts = misc.max_step(ds, dims, mnl, sigs) misc.scale2(lmbda, dz, dims, mnl) tz = misc.max_step(dz, dims, mnl, sigz) t = max([ 0.0, ts, tz ]) if t == 0: step = 1.0 else: step = min(1.0, STEP / t) # Backtrack until newx is in domain of f. backtrack = True while backtrack: xcopy(x, newx); xaxpy(dx, newx, alpha = step) t = F(newx) if t is None: newf = None else: newf, newDf = t[0], t[1] if newf is not None: backtrack = False else: step *= BETA # Merit function # # phi = theta1 * gap + theta2 * norm(rx) + # theta3 * norm(rznl) # # and its directional derivative dphi. phi = theta1 * gap + theta2 * resx + theta3 * resznl if i == 0: dphi = -phi else: dphi = -theta1 * (1 - sigma) * gap - \ theta2 * (1 - eta) * resx - \ theta3 * (1 - eta) * resznl # Line search. # # We use two types of line search. In a standard iteration we # use is a normal backtracking line search terminating with a # sufficient decrease of phi. In a "relaxed" iteration the # line search is terminated after one step, regardless of the # value of phi. We make at most MAX_RELAXED_ITERS consecutive # relaxed iterations. When starting a series of relaxed # iteration, we save the state at the end of the first line # search in the series (scaling, primal and dual iterates, # steps, sigma, eta, i.e., all information needed to resume # the line search at some later point). If a series of # MAX_RELAXED_ITERS relaxed iterations does not result in a # sufficient decrease compared to the value of phi at the start # of the series, then we resume the first line search in the # series as a standard line search (using the saved state). # On the other hand, if at some point during the series of # relaxed iterations we obtain a sufficient decrease of phi # compared with the value at the start of the series, then we # start a new series of relaxed line searches. # # To implement this we use a counter relaxed_iters. # # 1. If 0 <= relaxed_iters < MAX_RELAXED_ITERS, we use a # relaxed line search (one full step). If relaxed_iters # is 0, we save the value phi0 of the merit function at the # current point, as well as the state at the end of the # line search. If the relaxed line search results in a # sufficient decrease w.r.t. phi0, we reset relaxed_iters # to 0. Otherwise we increment relaxed_iters. # # 2. If relaxed_iters is MAX_RELAXED_ITERS, we use a standard # line search. If this results in a sufficient decrease # of the merit function compared with phi0, we set # relaxed_iters to 0. If phi decreased compared with phi0, # but not sufficiently, we set relaxed_iters to -1. # If phi increased compared with phi0, we resume the # backtracking at the last saved iteration, and after # completing the line search, set relaxed_iters to 0. # # 3. If relaxed_iters is -1, we use a standard line search # and increment relaxed_iters to 0. backtrack = True while backtrack: xcopy(x, newx); xaxpy(dx, newx, alpha = step) ycopy(y, newy); yaxpy(dy, newy, alpha = step) blas.copy(z, newz); blas.axpy(dz2, newz, alpha = step) blas.copy(s, news); blas.axpy(ds2, news, alpha = step) t = F(newx) newf, newDf = matrix(t[0], tc = 'd'), t[1] if type(newDf) is matrix or type(Df) is spmatrix: if newDf.typecode != 'd' or \ newDf.size != (mnl, c.size[0]): raise TypeError("second output argument "\ "of F() must be a 'd' matrix of size "\ "(%d,%d)" %(mnl, c.size[0])) def newfDf(u, v, alpha = 1.0, beta = 0.0, trans = 'N'): base.gemv(newDf, u, v, alpha = alpha, beta = beta, trans = trans) else: newfDf = newDf # newrx = c + A'*newy + newDf'*newz[:mnl] + G'*newz[mnl:] xcopy(c, newrx) fA(newy, newrx, beta = 1.0, trans = 'T') newfDf(newz[:mnl], newrx, beta = 1.0, trans = 'T') fG(newz[mnl:], newrx, beta = 1.0, trans = 'T') newresx = math.sqrt(xdot(newrx, newrx)) # newrznl = news[:mnl] + newf blas.copy(news[:mnl], newrznl) blas.axpy(newf, newrznl) newresznl = blas.nrm2(newrznl) newgap = (1.0 - (1.0 - sigma) * step) * gap \ + step**2 * dsdz newphi = theta1 * newgap + theta2 * newresx + \ theta3 * newresznl if i == 0: if newgap <= (1.0 - ALPHA * step) * gap and \ ( 0 <= relaxed_iters < MAX_RELAXED_ITERS or \ newphi <= phi + ALPHA * step * dphi ): backtrack = False sigma = min(newgap/gap, (newgap/gap) ** EXPON) eta = 0.0 else: step *= BETA else: if relaxed_iters == -1 or ( relaxed_iters == 0 == MAX_RELAXED_ITERS ): # Do a standard line search. if newphi <= phi + ALPHA * step * dphi: relaxed_iters == 0 backtrack = False else: step *= BETA elif relaxed_iters == 0 < MAX_RELAXED_ITERS: if newphi <= phi + ALPHA * step * dphi: # Relaxed l.s. gives sufficient decrease. relaxed_iters = 0 else: # Save state. phi0, dphi0, gap0 = phi, dphi, gap step0 = step blas.copy(W['dnl'], W0['dnl']) blas.copy(W['dnli'], W0['dnli']) blas.copy(W['d'], W0['d']) blas.copy(W['di'], W0['di']) for k in range(len(dims['q'])): blas.copy(W['v'][k], W0['v'][k]) W0['beta'][k] = W['beta'][k] for k in range(len(dims['s'])): blas.copy(W['r'][k], W0['r'][k]) blas.copy(W['rti'][k], W0['rti'][k]) xcopy(x, x0); xcopy(dx, dx0) ycopy(y, y0); ycopy(dy, dy0) blas.copy(s, s0); blas.copy(z, z0) blas.copy(ds, ds0) blas.copy(dz, dz0) blas.copy(ds2, ds20) blas.copy(dz2, dz20) blas.copy(lmbda, lmbda0) blas.copy(lmbdasq, lmbdasq0) dsdz0 = dsdz sigma0, eta0 = sigma, eta xcopy(rx, rx0); ycopy(ry, ry0) blas.copy(rznl, rznl0); blas.copy(rzl, rzl0) relaxed_iters = 1 backtrack = False elif 0 <= relaxed_iters < MAX_RELAXED_ITERS > 0: if newphi <= phi0 + ALPHA * step0 * dphi0: # Relaxed l.s. gives sufficient decrease. relaxed_iters = 0 else: # Relaxed line search relaxed_iters += 1 backtrack = False elif relaxed_iters == MAX_RELAXED_ITERS > 0: if newphi <= phi0 + ALPHA * step0 * dphi0: # Series of relaxed line searches ends # with sufficient decrease w.r.t. phi0. backtrack = False relaxed_iters = 0 else: # Resume last saved line search. phi, dphi, gap = phi0, dphi0, gap0 step = step0 blas.copy(W0['dnl'], W['dnl']) blas.copy(W0['dnli'], W['dnli']) blas.copy(W0['d'], W['d']) blas.copy(W0['di'], W['di']) for k in range(len(dims['q'])): blas.copy(W0['v'][k], W['v'][k]) W['beta'][k] = W0['beta'][k] for k in range(len(dims['s'])): blas.copy(W0['r'][k], W['r'][k]) blas.copy(W0['rti'][k], W['rti'][k]) xcopy(x0, x); xcopy(dx0, dx); ycopy(y0, y); ycopy(dy0, dy); blas.copy(s0, s); blas.copy(z0, z) blas.copy(ds0, ds) blas.copy(dz0, dz) blas.copy(ds20, ds2) blas.copy(dz20, dz2) blas.copy(lmbda0, lmbda) dsdz = dsdz0 sigma, eta = sigma0, eta0 relaxed_iters = -1 # Update x, y. xaxpy(dx, x, alpha = step) yaxpy(dy, y, alpha = step) # Replace nonlinear, 'l' and 'q' blocks of ds and dz with the # updated variables in the current scaling. # Replace 's' blocks of ds and dz with the factors Ls, Lz in a # factorization Ls*Ls', Lz*Lz' of the updated variables in the # current scaling. # ds := e + step*ds for nonlinear, 'l' and 'q' blocks. # dz := e + step*dz for nonlinear, 'l' and 'q' blocks. blas.scal(step, ds, n = mnl + dims['l'] + sum(dims['q'])) blas.scal(step, dz, n = mnl + dims['l'] + sum(dims['q'])) ind = mnl + dims['l'] ds[:ind] += 1.0 dz[:ind] += 1.0 for m in dims['q']: ds[ind] += 1.0 dz[ind] += 1.0 ind += m # ds := H(lambda)^{-1/2} * ds and dz := H(lambda)^{-1/2} * dz. # # This replaces the nonlinear, 'l' and 'q' components of ds and dz # with the updated variables in the new scaling. # The 's' components of ds and dz are replaced with # # diag(lmbda_k)^{1/2} * Qs * diag(lmbda_k)^{1/2} # diag(lmbda_k)^{1/2} * Qz * diag(lmbda_k)^{1/2} misc.scale2(lmbda, ds, dims, mnl, inverse = 'I') misc.scale2(lmbda, dz, dims, mnl, inverse = 'I') # sigs := ( e + step*sigs ) ./ lambda for 's' blocks. # sigz := ( e + step*sigz ) ./ lambda for 's' blocks. blas.scal(step, sigs) blas.scal(step, sigz) sigs += 1.0 sigz += 1.0 blas.tbsv(lmbda, sigs, n = sum(dims['s']), k = 0, ldA = 1, offsetA = mnl + dims['l'] + sum(dims['q']) ) blas.tbsv(lmbda, sigz, n = sum(dims['s']), k = 0, ldA = 1, offsetA = mnl + dims['l'] + sum(dims['q']) ) # dsk := Ls = dsk * sqrt(sigs). # dzk := Lz = dzk * sqrt(sigz). ind2, ind3 = mnl + dims['l'] + sum(dims['q']), 0 for k in range(len(dims['s'])): m = dims['s'][k] for i in range(m): blas.scal(math.sqrt(sigs[ind3+i]), ds, offset = ind2 + m*i, n = m) blas.scal(math.sqrt(sigz[ind3+i]), dz, offset = ind2 + m*i, n = m) ind2 += m*m ind3 += m # Update lambda and scaling. misc.update_scaling(W, lmbda, ds, dz) # Unscale s, z (unscaled variables are used only to compute # feasibility residuals). blas.copy(lmbda, s, n = mnl + dims['l'] + sum(dims['q'])) ind = mnl + dims['l'] + sum(dims['q']) ind2 = ind for m in dims['s']: blas.scal(0.0, s, offset = ind2) blas.copy(lmbda, s, offsetx = ind, offsety = ind2, n = m, incy = m+1) ind += m ind2 += m*m misc.scale(s, W, trans = 'T') blas.copy(lmbda, z, n = mnl + dims['l'] + sum(dims['q'])) ind = mnl + dims['l'] + sum(dims['q']) ind2 = ind for m in dims['s']: blas.scal(0.0, z, offset = ind2) blas.copy(lmbda, z, offsetx = ind, offsety = ind2, n = m, incy = m+1) ind += m ind2 += m*m misc.scale(z, W, inverse = 'I') gap = blas.dot(lmbda, lmbda) def cp(F, G = None, h = None, dims = None, A = None, b = None, kktsolver = None, xnewcopy = None, xdot = None, xaxpy = None, xscal = None, ynewcopy = None, ydot = None, yaxpy = None, yscal = None, **kwargs): """ Solves a convex optimization problem minimize f0(x) subject to fk(x) <= 0, k = 1, ..., mnl G*x <= h A*x = b. f = (f0, f1, ..., fmnl) is convex and twice differentiable. The linear inequalities are with respect to a cone C defined as the Cartesian product of N + M + 1 cones: C = C_0 x C_1 x .... x C_N x C_{N+1} x ... x C_{N+M}. The first cone C_0 is the nonnegative orthant of dimension ml. The next N cones are second order cones of dimension mq[0], ..., mq[N-1]. The second order cone of dimension m is defined as { (u0, u1) in R x R^{m-1} | u0 >= ||u1||_2 }. The next M cones are positive semidefinite cones of order ms[0], ..., ms[M-1] >= 0. Input arguments (basic usage). F is a function that handles the following arguments. F() returns a tuple (mnl, x0). mnl is the number of nonlinear inequality constraints. x0 is a point in the domain of f. F(x) returns a tuple (f, Df). f is a dense 'd' matrix of size (mnl+1,1) containing f(x). Df is a dense or sparse 'd' matrix of size (mnl+1, n), containing the derivatives of f at x: Df[k,:] is the transpose of the gradient of fk at x. If x is not in dom f, F(x) returns None or (None, None). F(x, z) with z a positive 'd' matrix of size (mnl+1,1), returns a tuple (f, Df, H). f and Df are defined as above. H is a dense or sparse 'd' matrix of size (n,n). The lower triangular part of H contains the lower triangular part of sum_k z[k] * Hk where Hk is the Hessian of fk at x. If F is called with two arguments, it can be assumed that x is dom f. If Df and H are returned as sparse matrices, their sparsity patterns must be the same for each call to F(x) or F(x,z). dims is a dictionary with the dimensions of the components of C. It has three fields. - dims['l'] = ml, the dimension of the nonnegative orthant C_0. (ml >= 0.) - dims['q'] = mq = [ mq[0], mq[1], ..., mq[N-1] ], a list of N integers with the dimensions of the second order cones C_1, ..., C_N. (N >= 0 and mq[k] >= 1.) - dims['s'] = ms = [ ms[0], ms[1], ..., ms[M-1] ], a list of M integers with the orders of the semidefinite cones C_{N+1}, ..., C_{N+M}. (M >= 0 and ms[k] >= 0.) The default value of dims = {'l': G.size[0], 'q': [], 's': []}. G is a dense or sparse 'd' matrix of size (K,n), where K = ml + mq[0] + ... + mq[N-1] + ms[0]**2 + ... + ms[M-1]**2. Each column of G describes a vector v = ( v_0, v_1, ..., v_N, vec(v_{N+1}), ..., vec(v_{N+M}) ) in V = R^ml x R^mq[0] x ... x R^mq[N-1] x S^ms[0] x ... x S^ms[M-1] stored as a column vector [ v_0; v_1; ...; v_N; vec(v_{N+1}); ...; vec(v_{N+M}) ]. Here, if u is a symmetric matrix of order m, then vec(u) is the matrix u stored in column major order as a vector of length m**2. We use BLAS unpacked 'L' storage, i.e., the entries in vec(u) corresponding to the strictly upper triangular entries of u are not referenced. h is a dense 'd' matrix of size (K,1), representing a vector in V, in the same format as the columns of G. A is a dense or sparse 'd' matrix of size (p,n). The default value is a sparse 'd' matrix of size (0,n). b is a dense 'd' matrix of size (p,1). The default value is a dense 'd' matrix of size (0,1). It is assumed that rank(A) = p and rank([H; A; Df; G]) = n at all x in dom f. The other arguments are normally not needed. They make it possible to exploit certain types of structure, as described below. Output arguments. cp() returns a dictionary with keys 'status', 'x', 'snl', 'sl', 'znl', 'zl', 'y', 'primal objective', 'dual objective', 'gap', 'relative gap', 'primal infeasibility', 'dual infeasibility', 'primal slack', 'dual slack'. The 'status' field has values 'optimal' or 'unknown'. If status is 'optimal', x, snl, sl, y, znl, zl are approximate solutions of the primal and dual optimality conditions f(x)[1:] + snl = 0, G*x + sl = h, A*x = b Df(x)'*[1; znl] + G'*zl + A'*y + c = 0 snl >= 0, znl >= 0, sl >= 0, zl >= 0 snl'*znl + sl'* zl = 0. If status is 'unknown', x, snl, sl, y, znl, zl are the last iterates before termination. They satisfy snl > 0, znl > 0, sl > 0, zl > 0, but are not necessarily feasible. The values of the other fields are the values returned by cpl() applied to the epigraph form problem minimize t subjec to f0(x) <= t fk(x) <= 0, k = 1, ..., mnl G*x <= h A*x = b. Termination with status 'unknown' indicates that the algorithm failed to find a solution that satisfies the specified tolerances. In some cases, the returned solution may be fairly accurate. If the primal and dual infeasibilities, the gap, and the relative gap are small, then x, y, snl, sl, znl, zl are close to optimal. Advanced usage. Three mechanisms are provided to express problem structure. 1. The user can provide a customized routine for solving linear equations (`KKT systems') [ sum_k zk*Hk(x) A' GG' ] [ ux ] [ bx ] [ A 0 0 ] [ uy ] = [ by ] [ GG 0 -W'*W ] [ uz ] [ bz ] where GG = [ Df[1:,:]; G ], uz = (uznl, uzl), bz = (bznl, bzl). z is a positive vector of length mnl+1 and x is a point in the domain of f. W is a scaling matrix, a block diagonal mapping W*u = ( Wnl*unl, W0*u_0, ..., W_{N+M}*u_{N+M} ) defined as follows. - For the nonlinear block (Wnl): Wnl = diag(dnl), with dnl a positive vector of length mnl. - For the 'l' block (W_0): W_0 = diag(d), with d a positive vector of length ml. - For the 'q' blocks (W_{k+1}, k = 0, ..., N-1): W_{k+1} = beta_k * ( 2 * v_k * v_k' - J ) where beta_k is a positive scalar, v_k is a vector in R^mq[k] with v_k[0] > 0 and v_k'*J*v_k = 1, and J = [1, 0; 0, -I]. - For the 's' blocks (W_{k+N}, k = 0, ..., M-1): W_k * u = vec(r_k' * mat(u) * r_k) where r_k is a nonsingular matrix of order ms[k], and mat(x) is the inverse of the vec operation. The optional argument kktsolver is a Python function that will be called as g = kktsolver(x, z, W). W is a dictionary that contains the parameters of the scaling: - W['dnl'] is a positive 'd' matrix of size (mnl, 1). - W['dnli'] is a positive 'd' matrix with the elementwise inverse of W['dnl']. - W['d'] is a positive 'd' matrix of size (ml, 1). - W['di'] is a positive 'd' matrix with the elementwise inverse of W['d']. - W['beta'] is a list [ beta_0, ..., beta_{N-1} ] - W['v'] is a list [ v_0, ..., v_{N-1} ] - W['r'] is a list [ r_0, ..., r_{M-1} ] - W['rti'] is a list [ rti_0, ..., rti_{M-1} ], with rti_k the inverse of the transpose of r_k. The call g = kktsolver(x, z, W) should return a function g that solves the KKT system by g(ux, uy, uz). On entry, ux, uy, uz contain the righthand side bx, by, bz. On exit, they contain the solution, with uz scaled: (Wnl*uznl, Wl*uzl) is returned instead of (uznl, uzl). 2. The linear operators Df*u, H*u, G*u and A*u can be specified by providing Python functions instead of matrices. This can only be done in combination with 1. above, i.e., it requires the kktsolver argument. If G is a function, the call G(u, v, alpha, beta, trans) should evaluate the matrix-vector products v := alpha * G * u + beta * v if trans is 'N' v := alpha * G' * u + beta * v if trans is 'T'. The arguments u and v are required. The other arguments have default values alpha = 1.0, beta = 0.0, trans = 'N'. If A is a function, the call A(u, v, alpha, beta, trans) should evaluate the matrix-vectors products v := alpha * A * u + beta * v if trans is 'N' v := alpha * A' * u + beta * v if trans is 'T'. The arguments u and v are required. The other arguments have default values alpha = 1.0, beta = 0.0, trans = 'N'. If Df is a function, the call Df(u, v, alpha, beta, trans) should evaluate the matrix-vectors products v := alpha * Df(x) * u + beta * v if trans is 'N' v := alpha * Df(x)' * u + beta * v if trans is 'T'. If H is a function, the call H(u, v, alpha, beta) should evaluate the matrix-vectors product v := alpha * H * u + beta * v. 3. Instead of using the default representation of the primal variable x and the dual variable y as one-column 'd' matrices, we can represent these variables and the corresponding parameters c and b by arbitrary Python objects (matrices, lists, dictionaries, etc). This can only be done in combination with 1. and 2. above, i.e., it requires a user-provided KKT solver and a function description of the linear mappings. It also requires the arguments xnewcopy, xdot, xscal, xaxpy, ynewcopy, ydot, yscal, yaxpy. These arguments are functions defined as follows. If X is the vector space of primal variables x, then: - xnewcopy(u) creates a new copy of the vector u in X. - xdot(u, v) returns the inner product of two vectors u and v in X. - xscal(alpha, u) computes u := alpha*u, where alpha is a scalar and u is a vector in X. - xaxpy(u, v, alpha = 1.0, beta = 0.0) computes v := alpha*u + v for a scalar alpha and two vectors u and v in X. If Y is the vector space of primal variables y: - ynewcopy(u) creates a new copy of the vector u in Y. - ydot(u, v) returns the inner product of two vectors u and v in Y. - yscal(alpha, u) computes u := alpha*u, where alpha is a scalar and u is a vector in Y. - yaxpy(u, v, alpha = 1.0, beta = 0.0) computes v := alpha*u + v for a scalar alpha and two vectors u and v in Y. Control parameters. The following control parameters can be modified by adding an entry to the dictionary options. options['show_progress'] True/False (default: True) options['maxiters'] positive integer (default: 100) options['refinement'] nonnegative integer (default: 1) options['abstol'] scalar (default: 1e-7) options['reltol'] scalar (default: 1e-6) options['feastol'] scalar (default: 1e-7). """ options = kwargs.get('options',globals()['options']) KKTREG = options.get('kktreg',None) import math from cvxopt import base, blas, misc from cvxopt.base import matrix, spmatrix mnl, x0 = F() # Argument error checking depends on level of customization. customkkt = type(kktsolver) is not str operatorG = G is not None and type(G) not in (matrix, spmatrix) operatorA = A is not None and type(A) not in (matrix, spmatrix) if (operatorG or operatorA) and not customkkt: raise ValueError("use of function valued G, A requires a "\ "user-provided kktsolver") customx = (xnewcopy != None or xdot != None or xaxpy != None or xscal != None) if customx and (not operatorG or not operatorA or not customkkt): raise ValueError("use of non-vector type for x requires "\ "function valued G, A and user-provided kktsolver") customy = (ynewcopy != None or ydot != None or yaxpy != None or yscal != None) if customy and (not operatorA or not customkkt): raise ValueError("use of non vector type for y requires "\ "function valued A and user-provided kktsolver") if not customx: if type(x0) is not matrix or x0.typecode != 'd' or x0.size[1] != 1: raise TypeError("'x0' must be a 'd' matrix with one column") if h is None: h = matrix(0.0, (0,1)) if type(h) is not matrix or h.typecode != 'd' or h.size[1] != 1: raise TypeError("'h' must be a 'd' matrix with one column") if not dims: dims = {'l': h.size[0], 'q': [], 's': []} # Dimension of the product cone of the linear inequalities. with 's' # components unpacked. cdim = dims['l'] + sum(dims['q']) + sum([ k**2 for k in dims['s'] ]) if h.size[0] != cdim: raise TypeError("'h' must be a 'd' matrix of size (%d,1)" %cdim) if G is None: if customx: def G(x, y, trans = 'N', alpha = 1.0, beta = 0.0): if trans == 'N': pass else: xscal(beta, y) else: G = spmatrix([], [], [], (0, x0.size[0])) if type(G) is matrix or type(G) is spmatrix: if G.typecode != 'd' or G.size != (cdim, x0.size[0]): raise TypeError("'G' must be a 'd' matrix with size (%d, %d)"\ %(cdim, x0.size[0])) def fG(x, y, trans = 'N', alpha = 1.0, beta = 0.0): misc.sgemv(G, x, y, dims, trans = trans, alpha = alpha, beta = beta) else: fG = G if A is None: if customy: def A(x, y, trans = 'N', alpha = 1.0, beta = 0.0): if trans == 'N': pass else: xscal(beta, y) else: A = spmatrix([], [], [], (0, x0.size[0])) if type(A) is matrix or type(A) is spmatrix: if A.typecode != 'd' or A.size[1] != x0.size[0]: raise TypeError("'A' must be a 'd' matrix with %d columns" \ %x0.size[0]) def fA(x, y, trans = 'N', alpha = 1.0, beta = 0.0): base.gemv(A, x, y, trans = trans, alpha = alpha, beta = beta) else: fA = A if not customy: if b is None: b = matrix(0.0, (0,1)) if type(b) is not matrix or b.typecode != 'd' or b.size[1] != 1: raise TypeError("'b' must be a 'd' matrix with one column") if not operatorA and b.size[0] != A.size[0]: raise TypeError("'b' must have length %d" %A.size[0]) if b is None and customy: raise ValueEror("use of non vector type for y requires b") if xnewcopy is None: xnewcopy = matrix if xdot is None: xdot = blas.dot if xaxpy is None: xaxpy = blas.axpy if xscal is None: xscal = blas.scal def xcopy(x, y): xscal(0.0, y) xaxpy(x, y) if ynewcopy is None: ynewcopy = matrix if ydot is None: ydot = blas.dot if yaxpy is None: yaxpy = blas.axpy if yscal is None: yscal = blas.scal def ycopy(x, y): yscal(0.0, y) yaxpy(x, y) # The problem is solved by applying cpl() to the epigraph form # # minimize t # subject to f0(x) - t <= 0 # f1(x) <= 0 # ... # fmnl(x) <= 0 # G*x <= h # A*x = b. # # The epigraph form variable is stored as a list [x, t]. # Epigraph form objective c = (0, 1). c = [ xnewcopy(x0), 1 ] xscal(0.0, c[0]) # Nonlinear inequalities for the epigraph problem # # f_e(x,t) = (f0(x) - t, f1(x), ..., fmnl(x)). # def F_e(x = None, z = None): if x is None: return mnl+1, [ x0, 0.0 ] else: if z is None: v = F(x[0]) if v is None or v[0] is None: return None, None val = matrix(v[0], tc = 'd') val[0] -= x[1] Df = v[1] else: val, Df, H = F(x[0], z) val = matrix(val, tc = 'd') val[0] -= x[1] if type(Df) in (matrix, spmatrix): def Df_e(u, v, alpha = 1.0, beta = 0.0, trans = 'N'): if trans == 'N': base.gemv(Df, u[0], v, alpha = alpha, beta = beta, trans = 'N') v[0] -= alpha * u[1] else: base.gemv(Df, u, v[0], alpha = alpha, beta = beta, trans = 'T') v[1] = -alpha * u[0] + beta * v[1] else: def Df_e(u, v, alpha = 1.0, beta = 0.0, trans = 'N'): if trans == 'N': Df(u[0], v, alpha = alpha, beta = beta, trans = 'N') v[0] -= alpha * u[1] else: Df(u, v[0], alpha = alpha, beta = beta, trans = 'T') v[1] = -alpha * u[0] + beta * v[1] if z is None: return val, Df_e else: if type(H) in (matrix, spmatrix): def H_e(u, v, alpha = 1.0, beta = 1.0): base.symv(H, u[0], v[0], alpha = alpha, beta = beta) v[1] += beta*v[1] else: def H_e(u, v, alpha = 1.0, beta = 1.0): H(u[0], v[0], alpha = alpha, beta = beta) v[1] += beta*v[1] return val, Df_e, H_e # Linear inequality constraints. # # G_e = [ G, 0 ] # if type(G) in (matrix, spmatrix): def G_e(u, v, alpha = 1.0, beta = 0.0, trans = 'N'): if trans == 'N': misc.sgemv(G, u[0], v, dims, alpha = alpha, beta = beta) else: misc.sgemv(G, u, v[0], dims, alpha = alpha, beta = beta, trans = 'T') v[1] *= beta else: def G_e(u, v, alpha = 1.0, beta = 0.0, trans = 'N'): if trans == 'N': G(u[0], v, alpha = alpha, beta = beta) else: G(u, v[0], alpha = alpha, beta = beta, trans = 'T') v[1] *= beta # Linear equality constraints. # # A_e = [ A, 0 ] # if type(A) in (matrix, spmatrix): def A_e(u, v, alpha = 1.0, beta = 0.0, trans = 'N'): if trans == 'N': base.gemv(A, u[0], v, alpha = alpha, beta = beta) else: base.gemv(A, u, v[0], alpha = alpha, beta = beta, trans = 'T') v[1] *= beta else: def A_e(u, v, alpha = 1.0, beta = 0.0, trans = 'N'): if trans == 'N': A(u[0], v, alpha = alpha, beta = beta) else: A(u, v[0], alpha = alpha, beta = beta, trans = 'T') v[1] *= beta # kktsolver(x, z, W) returns a routine for solving equations with # coefficient matrix # # [ H A' [Df[1:]; G]' ] # K = [ A 0 0 ]. # [ [Df[1:]; G] 0 -W'*W ] if kktsolver is None: if dims and (dims['q'] or dims['s']): kktsolver = 'chol' else: kktsolver = 'chol2' if kktsolver in ('ldl', 'chol', 'chol2', 'qr'): if kktsolver == 'ldl': factor = misc.kkt_ldl(G, dims, A, mnl, kktreg = KKTREG) elif kktsolver == 'qr': factor = misc.kkt_qr(G, dims, A, mnl) elif kktsolver == 'chol': factor = misc.kkt_chol(G, dims, A, mnl) else: factor = misc.kkt_chol2(G, dims, A, mnl) def kktsolver(x, z, W): f, Df, H = F(x, z) return factor(W, H, Df[1:,:]) ux, uz = xnewcopy(x0), matrix(0.0, (mnl + cdim, 1)) def kktsolver_e(x, znl, W): We = W.copy() We['dnl'] = W['dnl'][1:] We['dnli'] = W['dnli'][1:] g = kktsolver(x[0], znl, We) f, Df = F(x[0]) if type(Df) is matrix: gradf0 = Df[0,:].T elif type(Df) is spmatrix: gradf0 = matrix(Df[0,:].T) else: gradf0 = xnewcopy(x[0]) e0 = matrix(0.0, (mnl + 1, 1)) e0[0] = 1.0 Df(e0, gradf0, trans = 'T') def solve(x, y, z): # Solves # # [ [ H 0 ] [ A' ] [ Df' G'] ] [ ux ] [ bx ] # [ [ 0 0 ] [ 0 ] [ -e0' 0 ] ] [ ] [ ] # [ ] [ ] [ ] # [ [ A 0 ] 0 0 ] [ uy ] = [ by ]. # [ ] [ ] [ ] # [ [ Df -e0 ] 0 -W'*W ] [ uz ] [ bz ] # [ [ G 0 ] ] [ ] [ ] # # The solution is: # # uz[0] = -bx[1] # # [ ux[0] ] [ bx[0] + bx[1] * gradf0 ] # [ uy ] = K^-1 * [ by ]. # [ uz[1:] ] [ bz[1:] ] # # ux[1] = gradf0' * ux[0] - W['dnl'][0]**2 * uz[0] - bz[0] # = gradf0' * ux[0] + W['dnl'][0]**2 * bx[1] - bz[0]. # # Instead of uz we return the scaled solution W*uz. a = z[0] xcopy(x[0], ux) xaxpy(gradf0, ux, alpha = x[1]) blas.copy(z, uz, offsetx = 1) g(ux, y, uz) z[0] = -x[1] * W['dnl'][0] blas.copy(uz, z, offsety = 1) xcopy(ux, x[0]) x[1] = xdot(gradf0, x[0]) + W['dnl'][0]**2 * x[1] - a return solve def xnewcopy_e(x): return [ xnewcopy(x[0]), x[1] ] def xdot_e(x, y): return xdot(x[0], y[0]) + x[1]*y[1] def xaxpy_e(x, y, alpha = 1.0): xaxpy(x[0], y[0], alpha = alpha) y[1] += alpha*x[1] def xscal_e(alpha, x): xscal(alpha, x[0]) x[1] *= alpha sol = cpl(c, F_e, G_e, h, dims, A_e, b, kktsolver_e, xnewcopy_e, xdot_e, xaxpy_e, xscal_e, options = options) sol['x'] = sol['x'][0] sol['znl'], sol['snl'] = sol['znl'][1:], sol['snl'][1:] return sol def gp(K, F, g, G=None, h=None, A=None, b=None, kktsolver=None, **kwargs): """ Solves a geometric program minimize log sum exp (F0*x+g0) subject to log sum exp (Fi*x+gi) <= 0, i=1,...,m G*x <= h A*x = b Input arguments. K is a list of positive integers [K0, K1, K2, ..., Km]. F is a sum(K)xn dense or sparse 'd' matrix with block rows F0, F1, ..., Fm. Each Fi is Kixn. g is a sum(K)x1 dense or sparse 'd' matrix with blocks g0, g1, g2, ..., gm. Each gi is Kix1. G is an mxn dense or sparse 'd' matrix. h is an mx1 dense 'd' matrix. A is a pxn dense or sparse 'd' matrix. b is a px1 dense 'd' matrix. The default values for G, h, A and b are empty matrices with zero rows. Output arguments. Returns a dictionary with keys 'status', 'x', 'snl', 'sl', 'znl', 'zl', 'y', 'primal objective', 'dual objective', 'gap', 'relative gap', 'primal infeasibility', 'dual infeasibility', 'primal slack', 'dual slack'. The 'status' field has values 'optimal' or 'unknown'. If status is 'optimal', x, snl, sl, y, znl, zl are approximate solutions of the primal and dual optimality conditions f(x)[1:] + snl = 0, G*x + sl = h, A*x = b Df(x)'*[1; znl] + G'*zl + A'*y + c = 0 snl >= 0, znl >= 0, sl >= 0, zl >= 0 snl'*znl + sl'* zl = 0, where fk(x) = log sum exp (Fk*x + gk). If status is 'unknown', x, snl, sl, y, znl, zl are the last iterates before termination. They satisfy snl > 0, znl > 0, sl > 0, zl > 0, but are not necessarily feasible. The values of the other fields are the values returned by cpl() applied to the epigraph form problem minimize t subjec to f0(x) <= t fk(x) <= 0, k = 1, ..., mnl G*x <= h A*x = b. Termination with status 'unknown' indicates that the algorithm failed to find a solution that satisfies the specified tolerances. In some cases, the returned solution may be fairly accurate. If the primal and dual infeasibilities, the gap, and the relative gap are small, then x, y, snl, sl, znl, zl are close to optimal. Control parameters. The following control parameters can be modified by adding an entry to the dictionary options. options['show_progress'] True/False (default: True) options['maxiters'] positive integer (default: 100) options['refinement'] nonnegative integer (default: 1) options['abstol'] scalar (default: 1e-7) options['reltol'] scalar (default: 1e-6) options['feastol'] scalar (default: 1e-7). """ options = kwargs.get('options',globals()['options']) import math from cvxopt import base, blas, misc from cvxopt.base import matrix, spmatrix if type(K) is not list or [ k for k in K if type(k) is not int or k <= 0 ]: raise TypeError("'K' must be a list of positive integers") mnl = len(K)-1 l = sum(K) if type(F) not in (matrix, spmatrix) or F.typecode != 'd' or \ F.size[0] != l: raise TypeError("'F' must be a dense or sparse 'd' matrix "\ "with %d rows" %l) if type(g) is not matrix or g.typecode != 'd' or g.size != (l,1): raise TypeError("'g' must be a dene 'd' matrix of "\ "size (%d,1)" %l) n = F.size[1] if G is None: G = spmatrix([], [], [], (0,n)) if h is None: h = matrix(0.0, (0,1)) if type(G) not in (matrix, spmatrix) or G.typecode != 'd' or \ G.size[1] != n: raise TypeError("'G' must be a dense or sparse 'd' matrix "\ "with %d columns" %n) ml = G.size[0] if type(h) is not matrix or h.typecode != 'd' or h.size != (ml,1): raise TypeError("'h' must be a dense 'd' matrix of "\ "size (%d,1)" %ml) dims = {'l': ml, 's': [], 'q': []} if A is None: A = spmatrix([], [], [], (0,n)) if b is None: b = matrix(0.0, (0,1)) if type(A) not in (matrix, spmatrix) or A.typecode != 'd' or \ A.size[1] != n: raise TypeError("'A' must be a dense or sparse 'd' matrix "\ "with %d columns" %n) p = A.size[0] if type(b) is not matrix or b.typecode != 'd' or b.size != (p,1): raise TypeError("'b' must be a dense 'd' matrix of "\ "size (%d,1)" %p) y = matrix(0.0, (l,1)) u = matrix(0.0, (max(K),1)) Fsc = matrix(0.0, (max(K),n)) cs1 = [ sum(K[:i]) for i in range(mnl+1) ] cs2 = [ cs1[i] + K[i] for i in range(mnl+1) ] ind = list(zip(range(mnl+1), cs1, cs2)) def Fgp(x = None, z = None): if x is None: return mnl, matrix(0.0, (n,1)) f = matrix(0.0, (mnl+1,1)) Df = matrix(0.0, (mnl+1,n)) # y = F*x+g blas.copy(g, y) base.gemv(F, x, y, beta=1.0) if z is not None: H = matrix(0.0, (n,n)) for i, start, stop in ind: # yi := exp(yi) = exp(Fi*x+gi) ymax = max(y[start:stop]) y[start:stop] = base.exp(y[start:stop] - ymax) # fi = log sum yi = log sum exp(Fi*x+gi) ysum = blas.asum(y, n=stop-start, offset=start) f[i] = ymax + math.log(ysum) # yi := yi / sum(yi) = exp(Fi*x+gi) / sum(exp(Fi*x+gi)) blas.scal(1.0/ysum, y, n=stop-start, offset=start) # gradfi := Fi' * yi # = Fi' * exp(Fi*x+gi) / sum(exp(Fi*x+gi)) base.gemv(F, y, Df, trans='T', m=stop-start, incy=mnl+1, offsetA=start, offsetx=start, offsety=i) if z is not None: # Hi = Fi' * (diag(yi) - yi*yi') * Fi # = Fisc' * Fisc # where # Fisc = diag(yi)^1/2 * (I - 1*yi') * Fi # = diag(yi)^1/2 * (Fi - 1*gradfi') Fsc[:K[i], :] = F[start:stop, :] for k in range(start,stop): blas.axpy(Df, Fsc, n=n, alpha=-1.0, incx=mnl+1, incy=Fsc.size[0], offsetx=i, offsety=k-start) blas.scal(math.sqrt(y[k]), Fsc, inc=Fsc.size[0], offset=k-start) # H += z[i]*Hi = z[i] * Fisc' * Fisc blas.syrk(Fsc, H, trans='T', k=stop-start, alpha=z[i], beta=1.0) if z is None: return f, Df else: return f, Df, H return cp(Fgp, G, h, dims, A, b, kktsolver = kktsolver, options = options)