r"""
This module is responsible for conversion of a given symbolic optimization problem into semidefinite optimization
problems.
The main classes included in this module are:
+ `SDPRelaxations`
+ `SDRelaxSol`
+ `Mom`
"""
# from __future__ import print_function
from .base import base
from .sdp import sdp
from numpy import array, float64, ndarray, sqrt, zeros, abs, linalg, trim_zeros, where, random, dot
from numpy import zeros as npzeros
from numpy.random import uniform
from numpy.linalg import cholesky, LinAlgError
from sympy import Function, Symbol, QQ, groebner, Poly, zeros, reduced, sympify, Matrix, expand, latex, lambdify, Abs
from sympy.core.relational import Equality, GreaterThan, LessThan, StrictGreaterThan, StrictLessThan
from sympy.polys.polyerrors import PolynomialError
from sympy.polys.matrices import DomainMatrix
from scipy import optimize as opt
from scipy.linalg import eigvals
from scipy import linalg as spla
from math import ceil
from functools import reduce
from Irene.program import OptimizationProblem
from itertools import product
from operator import mul
from time import time
import multiprocessing as mp
from copy import copy
from pickle import load, loads, dump, dumps
[docs]
def Calpha_(expn, Mmnt):
r"""
Given an exponent `expn`, this function finds the corresponding
:math:`C_{expn}` matrix which can be used for parallel processing.
"""
r = Mmnt.shape[0]
C = zeros(r, r)
for i in range(r):
for j in range(i, r):
entity = Mmnt[i, j]
if expn in entity:
C[i, j] = entity[expn]
C[j, i] = C[i, j]
return array(C.tolist()).astype(float64)
[docs]
def Calpha__(expn, Mmnt, ii, q):
r"""
Given an exponent `expn`, this function finds the corresponding
:math:`C_{expn}` matrix which can be used for parallel processing.
"""
r = Mmnt.shape[0]
C = zeros(r, r)
for i in range(r):
for j in range(i, r):
entity = Mmnt[i, j].as_dict()
if expn in entity:
C[i, j] = entity[expn]
C[j, i] = C[i, j]
q.put([ii, array(C.tolist()).astype(float64)])
[docs]
class SDPRelaxations(base):
r"""
This class defines a function space by taking a family of sympy
symbolic functions and relations among them.
Simply, it initiates a commutative free real algebra on the symbolic
functions and defines the function space as the quotient of the free
algebra by the ideal generated by the given relations.
It takes three arguments:
- `gens` which is a list of ``sympy`` symbols and function symbols,
- `relations` which is a set of ``sympy`` expressions in terms of `gens` that defines an ideal.
- `name` is a given name which is used to save the state of the instant at break.
"""
GensError = r"""The `gens` must be a list of sympy functions or symbols"""
RelsError = r"""The `relations` must be a list of relation among generators"""
MonoOrdError = r"""`ord` must be one of 'lex', 'grlex', 'grevlex', 'ilex', 'igrlex', 'igrevlex'"""
MmntOrdError = r"""The order of moments must be a positive integer"""
SDPInpTypError = r"""The input of the SDP solver must be either a numpy matrix or ndarray"""
# Monomial order: "lex", "grlex", "grevlex", "ilex", "igrlex", "igrevlex"
MonomialOrder = 'lex'
SDPSolver = 'cvxopt'
Info = {}
ErrorTolerance = 1e-6
AvailableSolvers = []
PSDMoment = True
Probability = True
Parallel = True
def __init__(self, gens, relations=(), name="SDPRlx"):
assert type(gens) is list, self.GensError
assert type(gens) is list, self.RelsError
super(SDPRelaxations, self).__init__()
self.NumCores = mp.cpu_count()
self.EQ = Equality
self.GEQ = GreaterThan
self.LEQ = LessThan
self.GT = StrictGreaterThan
self.LT = StrictLessThan
self.ExpTypes = [Equality, GreaterThan,
LessThan, StrictGreaterThan, StrictLessThan]
self.Field = QQ
self.Generators = []
self.SymDict = {}
self.RevSymDict = {}
self.AuxSyms = []
self.NumGenerators = 0
self.FreeRelations = []
self.Groebner = []
self.MmntOrd = 0
self.ReducedBases = {}
#
self.Constraints = []
self.OrgConst = []
self.MomConst = []
self.OrgMomConst = []
self.ObjDeg = 0
self.ObjHalfDeg = 0
self.CnsDegs = []
self.CnsHalfDegs = []
self.MmntCnsDeg = 0
self.Blck = []
self.C_ = []
self.InitIdx = 0
self.LastIdxVal = 0
self.Stage = None
self.PrevStage = None
self.Name = name
self.SDP = sdp()
self.MatSize = []
self.InitTime = 0
self.Solution = None
self.f_min = 0
# check generators
for f in gens:
if isinstance(f, Function) or isinstance(f, Symbol):
self.Generators.append(f)
self.NumGenerators += 1
t_sym = Symbol('X%d' % self.NumGenerators)
self.SymDict[f] = t_sym
self.RevSymDict[t_sym] = f
self.AuxSyms.append(t_sym)
else:
raise TypeError(self.GensError)
self.Objective = Poly(0, *self.Generators)
self.RedObjective = Poly(0, *self.AuxSyms)
# check the relations
# TBI
for r in relations:
t_rel = r.subs(self.SymDict)
self.FreeRelations.append(t_rel)
if self.FreeRelations:
self.Groebner = groebner(
self.FreeRelations, domain=self.Field, order=self.MonomialOrder)
self.AvailableSolvers = self.AvailableSDPSolvers()
[docs]
def SetMonoOrd(self, ordr):
r"""
Changes the default monomial order to `ord` which mustbe among
`lex`, `grlex`, `grevlex`, `ilex`, `igrlex`, `igrevlex`.
"""
assert ordr in ['lex', 'grlex', 'grevlex', 'ilex', 'igrlex', 'igrevlex'], self.MonoOrdError
self.MonomialOrder = ordr
if self.FreeRelations:
self.Groebner = groebner(
self.FreeRelations, domain=self.Field, order=self.MonomialOrder)
[docs]
@classmethod
def from_problem(cls, optim_prob: OptimizationProblem, name="SDPRlx"):
"""
Creates an SDPRelaxations instance from an OptimizationProblem.
This method acts as an alternative constructor to bridge compatibility
with problems defined using SemigroupAlgebra.
Args:
optim_prob (OptimizationProblem): The optimization problem defined
with SemigroupAlgebra.
name (str): A name for the relaxation instance.
Returns:
An instance of SDPRelaxations.
"""
sga = optim_prob.sga
gen_names = sga.gens
sympy_gens = [Symbol(g) for g in gen_names]
sym_map = {name: sym for name, sym in zip(gen_names, sympy_gens)}
# Convert relations if they exist
relations = [optim_prob.to_sympy(rel, sym_map) for rel in optim_prob.relations] if optim_prob.relations else []
rlx = cls(sympy_gens, relations, name)
rlx.SetObjective(optim_prob.to_sympy(optim_prob.objective, sym_map))
for const in optim_prob.constraints:
rlx.AddConstraint(optim_prob.to_sympy(const, sym_map) >= 0)
return rlx
[docs]
def SetNumCores(self, num):
r"""
Sets the maximum number of workers which cannot be bigger than
number of available cores.
"""
assert (num > 0) and type(
num) is int, "Number of cores must be a positive integer."
self.NumCores = min(self.NumCores, num)
[docs]
def SetSDPSolver(self, solver):
r"""
Sets the default SDP solver. The followings are currently supported:
- CVXOPT
- DSDP
- SDPA
- CSDP
The selected solver must be installed otherwise it cannot be called.
The default solver is `CVXOPT` which has an interface for Python.
`DSDP` is called through the CVXOPT's interface. `SDPA` and `CSDP`
are called independently.
"""
assert solver.upper() in ['CVXOPT', 'DSDP', 'SDPA',
'CSDP'], "'%s' sdp solver is not supported" % solver
self.SDPSolver = solver
[docs]
def ReduceExp(self, expr):
r"""
Takes an expression `expr`, either in terms of internal free symbolic
variables or generating functions and returns the reduced expression
in terms of internal symbolic variables, if a relation among generators
is present, otherwise it just substitutes generating functions with
their corresponding internal symbols.
"""
try:
T = expr.subs(self.SymDict)
except:
T = Poly(expr, *self.AuxSyms)
if self.Groebner:
return reduced(T, self.Groebner)[1]
else:
return T
[docs]
def SetObjective(self, obj):
r"""
Takes the objective function `obj` as an algebraic combination
of the generating symbolic functions, replace the symbolic
functions with corresponding auxiliary symbols and reduce them
according to the given relations.
"""
self.Objective = sympify(obj)
self.RedObjective = self.ReduceExp(sympify(obj))
# self.CheckVars(obj)
tot_deg = Poly(self.RedObjective, *self.AuxSyms).total_degree()
self.ObjDeg = tot_deg
self.ObjHalfDeg = int(ceil(tot_deg / 2.))
[docs]
def AddConstraint(self, cnst):
r"""
Takes an (in)equality as an algebraic combination of the
generating functions that defines the feasibility region.
It reduces the defining (in)equalities according to the
given relations.
"""
self.OrgConst.append(cnst)
if isinstance(cnst, (self.GEQ, self.GT)):
non_red_exp = cnst.lhs - cnst.rhs
expr = self.ReduceExp(non_red_exp)
self.Constraints.append(expr)
tot_deg = Poly(expr, *self.AuxSyms).total_degree()
self.CnsDegs.append(tot_deg)
self.CnsHalfDegs.append(int(ceil(tot_deg / 2.)))
elif isinstance(cnst, (self.LEQ, self.LT)):
non_red_exp = cnst.rhs - cnst.lhs
expr = self.ReduceExp(non_red_exp)
self.Constraints.append(expr)
tot_deg = Poly(expr, *self.AuxSyms).total_degree()
self.CnsDegs.append(tot_deg)
self.CnsHalfDegs.append(int(ceil(tot_deg / 2.)))
elif isinstance(cnst, self.EQ):
non_red_exp = cnst.lhs - cnst.rhs
expr = self.ReduceExp(non_red_exp)
self.Constraints.append(self.ErrorTolerance + expr)
self.Constraints.append(self.ErrorTolerance - expr)
tot_deg = Poly(expr, *self.AuxSyms).total_degree()
# add twice
self.CnsDegs.append(tot_deg)
self.CnsDegs.append(tot_deg)
self.CnsHalfDegs.append(int(ceil(tot_deg / 2.)))
self.CnsHalfDegs.append(int(ceil(tot_deg / 2.)))
[docs]
def MomentConstraint(self, cnst):
r"""
Takes constraints on the moments. The input must be an instance of
`Mom` class.
"""
assert isinstance(
cnst, Mom), "The argument must be of moment type 'Mom'"
self.OrgMomConst.append(cnst)
CnsTyp = cnst.TYPE
if CnsTyp in ['ge', 'gt']:
expr = self.ReduceExp(cnst.Content)
tot_deg = Poly(expr, *self.AuxSyms).total_degree()
self.MmntCnsDeg = max(int(ceil(tot_deg / 2.)), self.MmntCnsDeg)
self.MomConst.append([expr, cnst.rhs])
elif CnsTyp in ['le', 'lt']:
expr = self.ReduceExp(-cnst.Content)
tot_deg = Poly(expr, *self.AuxSyms).total_degree()
self.MmntCnsDeg = max(int(ceil(tot_deg / 2.)), self.MmntCnsDeg)
self.MomConst.append([expr, -cnst.rhs])
elif CnsTyp == 'eq':
non_red_exp = cnst.Content - cnst.rhs
expr = self.ReduceExp(cnst.Content)
tot_deg = Poly(expr, *self.AuxSyms).total_degree()
self.MmntCnsDeg = max(int(ceil(tot_deg / 2.)), self.MmntCnsDeg)
self.MomConst.append([expr, cnst.rhs - self.ErrorTolerance])
self.MomConst.append([-expr, -cnst.rhs - self.ErrorTolerance])
[docs]
def ReducedMonomialBase(self, deg):
r"""
Returns a reduce monomial basis up to degree `d`.
"""
if deg in self.ReducedBases:
return self.ReducedBases[deg]
all_monos = product(range(deg + 1), repeat=self.NumGenerators)
req_monos = filter(lambda x: sum(x) <= deg, all_monos)
monos = [reduce(mul, [self.AuxSyms[i] ** expn[i]
for i in range(self.NumGenerators)], 1) for expn in req_monos]
RBase = []
for expr in monos:
rexpr = self.ReduceExp(expr)
expr_monos = Poly(rexpr, *self.AuxSyms).as_dict()
for mono_exp in expr_monos:
t_mono = reduce(mul, [self.AuxSyms[i] ** mono_exp[i]
for i in range(self.NumGenerators)], 1)
if t_mono not in RBase:
RBase.append(t_mono)
self.ReducedBases[deg] = RBase
return RBase
[docs]
def ExponentsVec(self, deg):
r"""
Returns all the exponents that appear in the reduced basis of all
monomials of the auxiliary symbols of degree at most `deg`.
"""
basis = self.ReducedMonomialBase(deg)
exponents = []
for elmnt in basis:
rbp = Poly(elmnt, *self.AuxSyms).as_dict()
for expnt in rbp:
if expnt not in exponents:
exponents.append(expnt)
return exponents
[docs]
def MomentsOrd(self, ordr):
r"""
Sets the order of moments to be considered.
"""
# from types import IntType
# assert (type(ordr) is IntType) and (ordr > 0), self.MmntOrdError
assert (type(ordr) is int) and (ordr > 0), self.MmntOrdError
self.MmntOrd = ordr
[docs]
def RelaxationDeg(self):
r"""
Finds the minimum required order of moments according to user's
request, objective function and constraints.
"""
if not self.CnsHalfDegs:
CHD = 0
else:
CHD = max(self.CnsHalfDegs)
RlxDeg = max([CHD, self.ObjHalfDeg, self.MmntOrd, self.MmntCnsDeg])
self.MmntOrd = RlxDeg
return RlxDeg
[docs]
def PolyCoefFullVec(self):
r"""
return the vector of coefficient of the reduced objective function
as an element of the vector space of elements of degree up to the
order of moments.
"""
c = []
fmono = Poly(self.RedObjective, *self.AuxSyms).as_dict()
exponents = self.ExponentsVec(2 * self.MmntOrd)
for expn in exponents:
if expn in fmono:
c.append(fmono[expn])
else:
c.append(0)
return c
def _poly_total_degree_or_raise(self, expr, context):
try:
return Poly(expr, *self.AuxSyms).total_degree()
except PolynomialError as exc:
raise ValueError("Unable to determine polynomial degree for %s" % context) from exc
[docs]
def LocalizedMoment(self, p):
r"""
Computes the reduced symbolic moment generating matrix localized
at `p`.
"""
tot_deg = self._poly_total_degree_or_raise(p, 'localized moment')
half_deg = int(ceil(tot_deg / 2.))
mmntord = self.MmntOrd - half_deg
m = Matrix(self.ReducedMonomialBase(mmntord))
LMmnt = expand(p * m * m.T)
LrMmnt = zeros(*LMmnt.shape)
for i in range(LMmnt.shape[0]):
for j in range(i, LMmnt.shape[1]):
LrMmnt[i, j] = self.ReduceExp(LMmnt[i, j])
LrMmnt[j, i] = LrMmnt[i, j]
return LrMmnt
[docs]
def LocalizedMoment_(self, p):
r"""
Computes the reduced symbolic moment generating matrix localized
at `p`.
"""
from sympy.polys.polymatrix import PolyMatrix
tot_deg = self._poly_total_degree_or_raise(p, 'localized moment')
half_deg = int(ceil(tot_deg / 2.))
mmntord = self.MmntOrd - half_deg
m = Matrix(self.ReducedMonomialBase(mmntord))
LMmnt = expand(p * m * m.T)
# LrMmnt = zeros(*LMmnt.shape)
tmp = [[0.*self.AuxSyms[0] for _ in range(LMmnt.shape[1])] for __ in range(LMmnt.shape[0])]
LrMmnt = tmp # PolyMatrix(tmp)
for i in range(LMmnt.shape[0]):
for j in range(i, LMmnt.shape[1]):
#LrMmnt[i, j] = Poly(self.ReduceExp(
#LMmnt[i, j]), *self.AuxSyms).as_dict()
#LrMmnt[j, i] = LrMmnt[i, j]
LrMmnt[i][j] = Poly(self.ReduceExp(
LMmnt[i, j]), *self.AuxSyms)
LrMmnt[j][i] = LrMmnt[i][j]
return PolyMatrix(LrMmnt)
[docs]
def MomentMat(self):
r"""
Returns the numerical moment matrix resulted from solving the SDP.
"""
assert 'moments' in self.Info, "The sdp has not been (successfully) solved (yet)."
Mmnt = self.LocalizedMoment(1.)
for i in range(Mmnt.shape[0]):
for j in range(Mmnt.shape[1]):
t_monos = Poly(Mmnt[i, j], *self.AuxSyms).as_dict()
t_mmnt = 0
for expn in t_monos:
mono = reduce(mul, [self.AuxSyms[k] ** expn[k]
for k in range(self.NumGenerators)], 1)
t_mmnt += t_monos[expn] * self.Info['moments'][mono]
Mmnt[i, j] = t_mmnt
Mmnt[j, i] = Mmnt[i, j]
return array(Mmnt.tolist()).astype(float64)
[docs]
def Calpha(self, expn, Mmnt):
r"""
Given an exponent `expn`, this method finds the corresponding
:math:`C_{expn}` matrix.
"""
r = Mmnt.shape[0]
C = zeros(r, r)
for i in range(r):
for j in range(i, r):
entity = Mmnt[i, j]
entity_monos = Poly(entity, *self.AuxSyms).as_dict()
if expn in entity_monos:
C[i, j] = entity_monos[expn]
C[j, i] = C[i, j]
return array(C.tolist()).astype(float64)
[docs]
def sInitSDP(self):
r"""
Initializes the semidefinite program (SDP), in serial mode, whose
solution is a lower bound for the minimum of the program.
"""
start = time()
self.SDP = sdp(self.SDPSolver)
self.RelaxationDeg()
N = len(self.ReducedMonomialBase(2 * self.MmntOrd))
self.MatSize = [len(self.ReducedMonomialBase(self.MmntOrd)), N]
Blck = [[] for _ in range(N)]
C = []
# Number of constraints
NumCns = len(self.CnsDegs)
# Number of moment constraints
NumMomCns = len(self.MomConst)
# Reduced vector of monomials of the given order
ExpVec = self.ExponentsVec(2 * self.MmntOrd)
# The localized moment matrices should be psd ##
for idx in range(NumCns):
d = len(self.ReducedMonomialBase(
self.MmntOrd - self.CnsHalfDegs[idx]))
# Corresponding C block is 0
h = zeros(d, d)
C.append(array(h.tolist()).astype(float64))
Mmnt = self.LocalizedMoment(self.Constraints[idx])
for i in range(N):
Blck[i].append(self.Calpha(ExpVec[i], Mmnt))
# Moment matrix should be psd ##
if self.PSDMoment:
d = len(self.ReducedMonomialBase(self.MmntOrd))
# Corresponding C block is 0
h = zeros(d, d)
C.append(array(h.tolist()).astype(float64))
Mmnt = self.LocalizedMoment(1.)
for i in range(N):
Blck[i].append(self.Calpha(ExpVec[i], Mmnt))
# L(1) = 1 #
if self.Probability:
for i in range(N):
Blck[i].append(array(
zeros(1, 1).tolist()).astype(float64))
Blck[i].append(array(
zeros(1, 1).tolist()).astype(float64))
# Blck[0][NumCns + 1][0] = 1
# Blck[0][NumCns + 2][0] = -1
Blck[0][-2][0] = 1
Blck[0][-1][0] = -1
C.append(array(Matrix([1]).tolist()).astype(float64))
C.append(array(Matrix([-1]).tolist()).astype(float64))
# Moment constraints
for idx in range(NumMomCns):
MomCns = Matrix([self.MomConst[idx][0]])
for i in range(N):
Blck[i].append(self.Calpha(ExpVec[i], MomCns))
C.append(array(
Matrix([self.MomConst[idx][1]]).tolist()).astype(float64))
self.SDP.C = C
self.SDP.b = self.PolyCoefFullVec()
self.SDP.A = Blck
elapsed = (time() - start)
self.InitTime = elapsed
[docs]
def Commit(self, blk, c, idx):
r"""
Sets the latest computed values for the final SDP
and saves the current state.
"""
self.Blck = copy(blk)
self.C_ = copy(c)
self.InitIdx = idx
def _commit_stage_state(self, blk, c, idx):
r"""
Commit stage state, preserving current interrupt semantics.
"""
try:
self.Commit(blk, c, idx)
except:
self.Commit(blk, c, idx)
raise KeyboardInterrupt
def _parallel_calpha_results(self, expvec, mmnt):
r"""
Computes ``Calpha__`` results in parallel and cleans up worker
processes deterministically.
"""
queue = mp.Queue(self.NumCores)
started_procs = []
results = [None for _ in range(len(expvec))]
try:
for idx, expn in enumerate(expvec):
proc = mp.Process(target=Calpha__,
args=(expn, mmnt, idx, queue))
proc.start()
started_procs.append(proc)
for _ in range(len(expvec)):
tmp = queue.get()
results[tmp[0]] = tmp[1]
return results
except BaseException:
for proc in started_procs:
if proc.is_alive():
proc.terminate()
raise
finally:
for proc in started_procs:
proc.join()
if hasattr(queue, 'close'):
queue.close()
if hasattr(queue, 'join_thread'):
queue.join_thread()
[docs]
def pInitSDP(self):
r"""
Initializes the semidefinite program (SDP), in parallel, whose
solution is a lower bound for the minimum of the program.
"""
start = time()
self.SDP = sdp(self.SDPSolver, solver_path=self.Path)
self.RelaxationDeg()
N = len(self.ReducedMonomialBase(2 * self.MmntOrd))
self.MatSize = [len(self.ReducedMonomialBase(self.MmntOrd)), N]
if not self.Blck:
self.Blck = [[] for _ in range(N)]
# Number of constraints
NumCns = len(self.CnsDegs)
# Number of moment constraints
NumMomCns = len(self.MomConst)
# Reduced vector of monomials of the given order
ExpVec = self.ExponentsVec(2 * self.MmntOrd)
# The localized moment matrices should be psd ##
if (self.PrevStage is None) or (self.PrevStage == "PSDLocMom"):
self.Stage = "PSDLocMom"
self.PrevStage = None
idx = self.LastIdxVal
while idx < NumCns:
d = len(self.ReducedMonomialBase(
self.MmntOrd - self.CnsHalfDegs[idx]))
Mmnt = self.LocalizedMoment_(self.Constraints[idx])
results = self._parallel_calpha_results(ExpVec, Mmnt)
# stash changes
tBlck = copy(self.Blck)
tC_ = copy(self.C_)
for i in range(N):
tBlck[i].append(results[i])
# Corresponding self.C_ block is 0
h = zeros(d, d)
tC_.append(array(h.tolist()).astype(float64))
# increase loop counter
idx += 1
# commit changes
self._commit_stage_state(tBlck, tC_, idx)
# self.InitIdx = idx
self.LastIdxVal = 0
# Moment matrix should be psd ##
if (self.PrevStage is None) or (self.PrevStage == "PSDMom"):
self.Stage = "PSDMom"
self.PrevStage = None
if self.PSDMoment:
d = len(self.ReducedMonomialBase(self.MmntOrd))
Mmnt = self.LocalizedMoment_(1.)
results = self._parallel_calpha_results(ExpVec, Mmnt)
# stash changes
tBlck = copy(self.Blck)
tC_ = copy(self.C_)
for i in range(N):
tBlck[i].append(results[i])
# Corresponding self.C_ block is 0
h = zeros(d, d)
tC_.append(array(h.tolist()).astype(float64))
# commit changes
self._commit_stage_state(tBlck, tC_, 0)
# self.Blck = copy(tBlck)
# self.C_ = copy(tC_)
# L(1) = 1 ##
if (self.PrevStage is None) or (self.PrevStage == "L(1)=1"):
self.Stage = "L(1)=1"
self.PrevStage = None
if self.Probability:
# stash changes
tBlck = copy(self.Blck)
tC_ = copy(self.C_)
for i in range(N):
tBlck[i].append(array(
zeros(1, 1).tolist()).astype(float64))
tBlck[i].append(array(
zeros(1, 1).tolist()).astype(float64))
# Blck[0][NumCns + 1][0] = 1
# Blck[0][NumCns + 2][0] = -1
tBlck[0][-2][0] = 1
tBlck[0][-1][0] = -1
tC_.append(array(Matrix([1]).tolist()).astype(float64))
tC_.append(
array(Matrix([-1]).tolist()).astype(float64))
# commit changes
self._commit_stage_state(tBlck, tC_, 0)
# self.Blck = copy(tBlck)
# self.C_ = copy(tC_)
# Moment constraints
if (self.PrevStage is None) or (self.PrevStage == "MomConst"):
self.Stage = "MomConst"
self.PrevStage = None
idx = self.LastIdxVal
while idx < NumMomCns:
MomCns = Matrix([self.MomConst[idx][0]])
# stash changes
tBlck = copy(self.Blck)
tC_ = copy(self.C_)
for i in range(N):
tBlck[i].append(self.Calpha(ExpVec[i], MomCns))
tC_.append(array(
Matrix([self.MomConst[idx][1]]).tolist()).astype(float64))
# increase loop counter
idx += 1
# commit changes
self._commit_stage_state(tBlck, tC_, idx)
# self.Blck = copy(tBlck)
# self.C_ = copy(tC_)
# self.InitIdx = idx
self.SDP.C = self.C_
self.SDP.b = self.PolyCoefFullVec()
self.SDP.A = self.Blck
elapsed = (time() - start)
self.InitTime = elapsed
[docs]
def InitSDP(self):
r"""
Initializes the SDP based on the value of ``self.Parallel``.
If it is ``True``, it runs in parallel mode, otherwise
in serial.
"""
if self.Parallel:
try:
self.pInitSDP()
except KeyboardInterrupt:
with open(self.Name + '.rlx', 'wb') as obj_file:
dump(self, obj_file)
print("\n...::: The program is saved in '" +
self.Name + ".rlx' :::...")
raise KeyboardInterrupt
else:
self.sInitSDP()
[docs]
def Minimize(self):
r"""
Finds the minimum of the truncated moment problem which provides
a lower bound for the actual minimum.
"""
self.SDP.solve()
self.Solution = SDRelaxSol(
self.AuxSyms, symdict=self.SymDict, err_tol=self.ErrorTolerance)
self.Info = {}
self.Solution.Status = self.SDP.Info['Status']
if self.SDP.Info['Status'] == 'Optimal':
self.f_min = min(self.SDP.Info['PObj'], self.SDP.Info['DObj'])
self.Solution.Primal = self.SDP.Info['PObj']
self.Solution.Dual = self.SDP.Info['DObj']
self.Info = {"min": self.f_min, "CPU": self.SDP.Info[
'CPU'], 'InitTime': self.InitTime}
self.Solution.RunTime = self.SDP.Info['CPU']
self.Solution.InitTime = self.InitTime
self.Info['status'] = 'Optimal'
self.Info[
'Message'] = 'Feasible solution for moments of order ' + str(self.MmntOrd)
self.Solution.Message = self.Info['Message']
self.Info['tms'] = self.SDP.Info['y']
FullMonVec = self.ReducedMonomialBase(2 * self.MmntOrd)
self.Info['moments'] = {FullMonVec[i]: self.Info[
'tms'][i] for i in range(len(FullMonVec))}
self.Info['solver'] = self.SDP.solver
for idx in self.Info['moments']:
self.Solution.TruncatedMmntSeq[idx.subs(self.RevSymDict)] = self.Info[
'moments'][idx]
self.Solution.MomentMatrix = self.MomentMat()
self.Solution.MonoBase = self.ReducedMonomialBase(self.MmntOrd)
self.Solution.Solver = self.SDP.solver
self.Solution.NumGenerators = self.NumGenerators
else:
self.f_min = None
self.Info['min'] = self.f_min
self.Info['status'] = 'Infeasible'
self.Info['Message'] = 'No feasible solution for moments of order ' + \
str(self.MmntOrd) + ' were found'
self.Solution.Status = 'Infeasible'
self.Solution.Message = self.Info['Message']
self.Solution.Solver = self.SDP.solver
self.Info["Size"] = self.MatSize
return self.f_min
[docs]
def Decompose(self):
r"""
Returns a dictionary that associates a list to every constraint,
:math:`g_i\ge0` for :math:`i=0,\dots,m`, where :math:`g_0=1`.
Each list consists of elements of algebra whose sums of squares
is equal to :math:`\sigma_i` and :math:`f-f_*=\sum_{i=0}^m\sigma_ig_i`.
Here, :math:`f_*` is the output of the ``SDPRelaxation.Minimize()``.
"""
SOSCoefs = {}
blks = []
NumCns = len(self.CnsDegs)
for M in self.SDP.Info['X']:
blks.append(Matrix(cholesky(M)))
for idx in range(NumCns):
SOSCoefs[idx + 1] = []
v = Matrix(self.ReducedMonomialBase(
self.MmntOrd - self.CnsHalfDegs[idx])).T
decomp = v * blks[idx]
for p in decomp:
SOSCoefs[idx + 1].append(p.subs(self.RevSymDict))
v = Matrix(self.ReducedMonomialBase(self.MmntOrd)).T
SOSCoefs[0] = []
decomp = v * blks[NumCns]
for p in decomp:
SOSCoefs[0].append(p.subs(self.RevSymDict))
return SOSCoefs
[docs]
def getObjective(self):
r"""
Returns the objective function of the problem after reduction modulo the relations, if given.
"""
return self.RedObjective.subs(self.RevSymDict)
[docs]
def getConstraint(self, idx):
r"""
Returns the constraint number `idx` of the problem after reduction modulo the relations, if given.
"""
assert idx < len(self.Constraints), "Index out of range."
return self.Constraints[idx].subs(self.RevSymDict) >= 0
[docs]
def getMomentConstraint(self, idx):
r"""
Returns the moment constraint number `idx` of the problem after reduction modulo the relations, if given.
"""
assert idx < len(self.MomConst), "Index out of range."
return self.MomConst[idx][0].subs(self.RevSymDict) >= sympify(self.MomConst[idx][1]).subs(self.RevSymDict)
[docs]
def Resume(self):
r"""
Resumes the process of a previously saved program.
"""
with open(self.Name + '.rlx', 'rb') as obj_file:
return load(obj_file)
[docs]
def SaveState(self):
r"""
Saves the current state of the relaxation object to the file `self.Name+'.rlx'`.
"""
with open(self.Name + '.rlx', 'wb') as obj_file:
dump(self, obj_file)
[docs]
def State(self):
r"""
Returns the latest state of the object at last break and save.
"""
with open(self.Name + '.rlx', 'rb') as obj_file:
ser_inst = load(obj_file)
return ser_inst.PrevStage, ser_inst.LastIdxVal
def __str__(self):
r"""
String output.
"""
out_txt = "=" * 70 + "\n"
out_txt += "Minimize\t"
out_txt += str(self.RedObjective.subs(self.RevSymDict)) + "\n"
out_txt += "Subject to\n"
for cns in self.Constraints:
out_txt += "\t\t" + str(cns.subs(self.RevSymDict) >= 0) + "\n"
out_txt += "And\n"
for cns in self.MomConst:
out_txt += "\t\tMoment " + \
str(cns[0].subs(self.RevSymDict) >= sympify(
cns[1]).subs(self.RevSymDict)) + "\n"
out_txt += "=" * 70 + "\n"
return out_txt
def __getstate__(self):
r"""
Pickling process.
"""
self.PrevStage = self.Stage
self.LastIdxVal = self.InitIdx
# self.SDP = None
exceptions = ['RevSymDict', 'Generators', 'Objective',
'SymDict', 'Solution', 'OrgConst']
cur_inst = self.__dict__
ser_inst = {}
for kw in cur_inst:
if kw in exceptions:
ser_inst[kw] = str(cur_inst[kw])
else:
ser_inst[kw] = dumps(cur_inst[kw])
return dumps(ser_inst)
def __setstate__(self, state):
r"""
Loading pickles
"""
exceptions = ['RevSymDict', 'Generators', 'Objective',
'SymDict', 'Solution', 'OrgConst']
ser_inst = loads(state)
for kw in ser_inst:
if kw in exceptions:
if kw not in ['Solution']:
self.__dict__[kw] = sympify(ser_inst[kw])
else:
self.__dict__[kw] = loads(ser_inst[kw])
def __latex__(self):
r"""
Generates LaTeX code of the optimization problem.
"""
latexcode = "\\left\\lbrace\n"
latexcode += "\\begin{array}{ll}\n"
latexcode += "\t\\min & " + latex(self.Objective) + "\\\\\n"
latexcode += "\t\\textrm{subject to} & \\\\\n"
for cns in self.OrgConst:
latexcode += "\t\t & " + latex(cns) + "\\\\\n"
latexcode += "\t\\textrm{where} & \\\\\n"
for cns in self.OrgMomConst:
latexcode += "\t\t" + cns.__latex__(True) + "\\\\\n"
latexcode += "\\end{array}"
latexcode += "\\right."
return latexcode
#######################################################################
# Solution of the Semidefinite Relaxation
[docs]
class SDRelaxSol(object):
r"""
Instances of this class carry information on the solution of the
semidefinite relaxation associated to an optimization problem.
It includes various pieces of information:
- ``SDRelaxSol.TruncatedMmntSeq`` a dictionary of resulted moments
- ``SDRelaxSol.MomentMatrix`` the resulted moment matrix
- ``SDRelaxSol.Primal`` the value of the SDP in primal form
- ``SDRelaxSol.Dual`` the value of the SDP in dual form
- ``SDRelaxSol.RunTime`` the run time of the sdp solver
- ``SDRelaxSol.InitTime`` the total time consumed for initialization of the sdp
- ``SDRelaxSol.Solver`` the name of sdp solver
- ``SDRelaxSol.Status`` final status of the sdp solver
- ``SDRelaxSol.RelaxationOrd`` order of relaxation
- ``SDRelaxSol.Message`` the message that maybe returned by the sdp solver
- ``SDRelaxSol.ScipySolver`` the scipy solver to extract solutions
- ``SDRelaxSol.err_tol`` the minimum value which is considered to be nonzero
- ``SDRelaxSol.Support`` the support of discrete measure resulted from ``SDPRelaxation.Minimize()``
- ``SDRelaxSol.Weights`` corresponding weights for the Dirac measures
"""
def __init__(self, X, symdict={}, err_tol=10e-6):
self.TruncatedMmntSeq = {}
self.MomentMatrix = None
self.Primal = None
self.Dual = None
self.RunTime = None
self.InitTime = None
self.Solver = None
self.Status = None
self.RelaxationOrd = None
self.Message = None
self.MonoBase = None
self.NumGenerators = None
# SDPRelaxations auxiliary symbols
self.X = X
self.Xij = None
self.SymDict = symdict
self.RevSymDict = {}
for v in self.SymDict:
self.RevSymDict[self.SymDict[v]] = v
self.err_tol = err_tol
self.ScipySolver = 'lm'
self.Support = None
self.Weights = None
self.weight = []
def __str__(self):
r"""
Generate the output for print.
"""
out_str = "Solution of a Semidefinite Program:\n"
out_str += " Solver: " + self.Solver + "\n"
out_str += " Status: " + self.Status + "\n"
out_str += " Initialization Time: " + \
str(self.InitTime) + " seconds\n"
out_str += " Run Time: " + \
str(self.RunTime) + " seconds\n"
out_str += "Primal Objective Value: " + str(self.Primal) + "\n"
out_str += " Dual Objective Value: " + str(self.Dual) + "\n"
if self.Support is not None:
out_str += " Support:\n"
for p in self.Support:
out_str += "\t\t" + str(p) + "\n"
out_str += " Support solver: " + self.ScipySolver + "\n"
out_str += self.Message + "\n"
return out_str
def __getitem__(self, idx):
r"""
Returns the moment corresponding to the index ``idx`` if exists,
otherwise, returns ``None``.
"""
if idx in self.TruncatedMmntSeq:
return self.TruncatedMmntSeq[idx]
else:
return None
def __len__(self):
r"""
Returns the length of the moment sequence.
"""
return len(self.TruncatedMmntSeq)
def __iter__(self):
return iter(self.TruncatedMmntSeq)
[docs]
def SetScipySolver(self, solver):
r"""
Sets the ``scipy.optimize.root`` solver to `solver`.
"""
assert solver.lower() in ['hybr', 'lm', 'broyden1', 'broyden2', 'anderson', 'linearmixing', 'diagbroyden',
'excitingmixing', 'krylov',
'df-sane'], "Unrecognized solver. The solver must be among 'hybr', 'lm', 'broyden1', 'broyden2', 'anderson', 'linearmixing', 'diagbroyden', 'excitingmixing', 'krylov', 'df-sane'"
self.ScipySolver = solver.lower()
[docs]
def Pivot(self, arr):
r"""
Get the leading term of each column.
"""
if arr.ndim == 1:
idxs, = arr.nonzero()
elif arr.ndim == 2:
assert arr.shape[0] == 1
idxs = zip(*arr.nonzero())
else:
raise Exception("Array of unexpected size: " + arr.ndim)
for idx in idxs:
elem = arr[idx]
if abs(elem) > self.err_tol:
if arr.ndim == 1:
return idx, elem
elif arr.ndim == 2:
return idx[1], elem
return 0, arr[0]
[docs]
def StblRedEch(self, A):
r"""
Compute the stabilized row reduced echelon form.
"""
A = array(A)
m, n = A.shape
Q = []
R = npzeros((min(m, n), n)) # Rows
for i, ai in enumerate(A.T):
# Remove any contribution from previous rows
for j, qj in enumerate(Q):
R[j, i] = ai.dot(qj)
ai -= ai.dot(qj) * qj
li = sqrt((ai ** 2).sum())
if li > self.err_tol:
assert len(Q) < min(m, n)
# Add a new column to Q
Q.append(ai / li)
# And write down the contribution
R[len(Q) - 1, i] = li
# Convert to reduced row echelon form
nrows, _ = R.shape
for i in range(nrows - 1, 0, -1):
k, v = self.Pivot(R[i, :])
if v > self.err_tol:
for j in range(i):
R[j, :] -= R[i, :] * R[j, k] / R[i, k]
else:
R[i, :] = 0
# row_normalize
for r in R:
li = sqrt((r ** 2).sum())
if li < self.err_tol:
r[:] = 0
else:
r /= li
return array(Q).T, R
[docs]
def NumericalRank(self):
r"""
Finds the rank of the moment matrix based on the size of its
eigenvalues. It considers those with absolute value less than
``self.err_tol`` to be zero.
"""
num_rnk = 0
eignvls = eigvals(self.MomentMatrix)
for ev in eignvls:
if abs(ev) >= self.err_tol:
num_rnk += 1
return num_rnk
[docs]
def Term2Mmnt(self, trm, rnk, X):
r"""
Converts a moment object into an algebraic equation.
"""
num_vars = len(X)
expr = 0
for i in range(rnk):
expr += self.weight[i] * \
trm.subs({X[j]: self.Xij[i][j] for j in range(num_vars)})
return expr
def __latex__(self):
r"""
Generates LaTeX code for the moment matrix.
"""
a = self.MomentMatrix
lines = str(a).replace('[', '').replace(']', '').splitlines()
rv = [r'\begin{bmatrix}']
rv += [' ' + ' & '.join(l.split()) + r'\\' for l in lines]
rv += [r'\end{bmatrix}']
return '\n'.join(rv)
#######################################################################
# A Symbolic object to handle moment constraints
[docs]
class Mom(object):
r"""
This is a simple interface to define moment constraints to be
used via `SDPRelaxations.MomentConstraint`.
It takes a sympy expression as input and initiates an object
which can be used to force particular constraints on the moment
sequence.
**Example:** Force the moment of :math:`x^2f(x) + f(x)^2` to be at least `.5`::
Mom(x**2 * f + f**2) >= .5
# OR
Mom(x**2 * f) + Mom(f**2) >= .5
"""
def __init__(self, expr):
# from types import IntType, LongType, FloatType
# self.NumericTypes = [IntType, LongType, FloatType]
self.NumericTypes = [int, float]
self.Content = sympify(expr)
self.rhs = 0
self.TYPE = None
def __add__(self, x):
if isinstance(x, Mom):
self.Content += x.Content
else:
self.Content += x
return self
def __sub__(self, x):
if isinstance(x, Mom):
self.Content -= x.Content
else:
self.Content -= x
return self
def __neg__(self):
self.Content = -self.Content
return self
def __mul__(self, x):
if type(x) in self.NumericTypes:
self.Content = x * self.Content
else:
raise Exception("Operation not supported")
return self
def __rmul__(self, x):
if type(x) in self.NumericTypes:
self.Content = x * self.Content
else:
raise Exception("Operation not supported")
return self
def __ge__(self, x):
if isinstance(x, Mom):
self.rhs = 0
self.Content -= x.Content
elif type(x) in self.NumericTypes:
self.rhs = x
self.TYPE = 'ge'
return self
def __gt__(self, x):
if isinstance(x, Mom):
self.rhs = 0
self.Content -= x.Content
elif type(x) in self.NumericTypes:
self.rhs = x
self.TYPE = 'gt'
return self
def __le__(self, x):
if isinstance(x, Mom):
self.rhs = 0
self.Content -= x.Content
elif type(x) in self.NumericTypes:
self.rhs = x
self.TYPE = 'le'
return self
def __lt__(self, x):
if isinstance(x, Mom):
self.rhs = 0
self.Content -= x.Content
elif type(x) in self.NumericTypes:
self.rhs = x
self.TYPE = 'lt'
return self
def __eq__(self, x):
if isinstance(x, Mom):
self.rhs = 0
self.Content -= x.Content
elif type(x) in self.NumericTypes:
self.rhs = x
else:
self.rhs += x
self.TYPE = 'eq'
return self
def __str__(self):
symbs = {'lt': '<', 'le': '<=', 'gt': '>', 'ge': '>=', 'eq': '=='}
strng = str(self.Content)
if self.TYPE is not None:
strng += " " + symbs[self.TYPE]
strng += " " + str(self.rhs)
return strng
def __getstate__(self):
r"""
Pickling process.
"""
ser_inst = {'NumericTypes': dumps(self.NumericTypes), 'Content': str(self.Content), 'rhs': dumps(self.rhs),
'TYPE': dumps(self.TYPE)}
return dumps(ser_inst)
def __setstate__(self, state):
r"""
Loading pickles
"""
ser_inst = loads(state)
self.__dict__['NumericTypes'] = loads(ser_inst['NumericTypes'])
self.__dict__['Content'] = sympify(ser_inst['Content'])
self.__dict__['rhs'] = loads(ser_inst['rhs'])
self.__dict__['TYPE'] = loads(ser_inst['TYPE'])
def __latex__(self, external=False):
r"""
Generates LaTeX code for the moment term.
"""
symbs = {'lt': '<', 'le': '\\leq', 'gt': '>', 'ge': '\\geq', 'eq': '='}
latexcode = "\\textrm{Moment of }"
if external:
latexcode += " & "
latexcode += latex(self.Content)
latexcode += symbs[self.TYPE] + latex(self.rhs)
return latexcode
