Skip to content

Commit

Permalink
removed avoid_slacks as a compilation option
Browse files Browse the repository at this point in the history
  • Loading branch information
rileyjmurray committed Oct 5, 2019
1 parent f95ce1e commit cc8421c
Show file tree
Hide file tree
Showing 4 changed files with 13 additions and 440 deletions.
32 changes: 5 additions & 27 deletions sageopt/coniclifts/problems/solvers/mosek.py
View file
Expand Up @@ -54,18 +54,9 @@ def apply(c, A, b, K, destructive, compilation_options):
b = b.copy()
K = copy.deepcopy(K)
c = c.copy()
if 'avoid_slacks' in compilation_options:
avoid_slacks = compilation_options['avoid_slacks']
else:
avoid_slacks = False
A, b, K, sep_K, scale, trans = separate_cone_constraints(A, b, K,
destructive=True,
dont_sep={'0', '+'},
avoid_slacks=avoid_slacks)
# A,b,K now reflect a conic system where "x" has been replaced by "np.diag(scale) @ (x + trans)"
inv_data = {'n': A.shape[1]}
A, b, K, sep_K = separate_cone_constraints(A, b, K, destructive=True, dont_sep={'0', '+'})
c = np.hstack([c, np.zeros(shape=(A.shape[1] - len(c)))])
objective_offset = c @ trans
c = c * scale
type_selectors = build_cone_type_selectors(K)
# Below: inequality constraints that the "user" intended to give to MOSEK.
A_ineq = A[type_selectors['+'], :]
Expand All @@ -81,7 +72,6 @@ def apply(c, A, b, K, destructive, compilation_options):
K = [Cone('+', A_ineq.shape[0]), Cone('0', A_z.shape[0])]
# Return values
data = {'A': A, 'b': b, 'K': K, 'sep_K': sep_K, 'c': c}
inv_data = {'scaling': scale, 'translation': trans, 'objective_offset': objective_offset, 'n': A.shape[1]}
return data, inv_data

@staticmethod
Expand Down Expand Up @@ -141,11 +131,7 @@ def parse_result(solver_output, inv_data, var_mapping):
:param solver_output: a dictionary containing the mosek Task and Environment that
were constructed at solve-time. Also contains any parameters that were passed at solved-time.
:param inv_data: a dictionary with keys 'scaling' and 'translation'. The primal variable
values returned by Mosek should be ...
(1) scaled by the inverse of inv_data['scaling'], and
(2) translated by the negative of inv_data['translation'].
Then the objective value should be reduced by inv_data['objective_offset'].
:param inv_data: a dictionary with key 'n'.
:param var_mapping: a dictionary mapping names of coniclifts Variables to arrays
of indices. The array var_mapping['my_var'] contains the indices in
Expand All @@ -164,7 +150,7 @@ def parse_result(solver_output, inv_data, var_mapping):
if solution_status == mosek.solsta.optimal:
# optimal
problem_status = CL_CONSTANTS.solved
problem_value = task.getprimalobj(sol) - inv_data['objective_offset']
problem_value = task.getprimalobj(sol)
x0 = [0.] * inv_data['n']
task.getxxslice(sol, 0, len(x0), x0)
x0 = np.array(x0)
Expand Down Expand Up @@ -237,15 +223,7 @@ def set_default_solver_settings(task):

@staticmethod
def load_variable_values(x0, inv_data, var_mapping):
x = inv_data['scaling'] * x0 + inv_data['translation']
# ^ That produces the correct result for a test, but isn't consistent with
# comments written in the reformulators.py file.
#
# x = np.power(inv_data['scaling'], -1) * x0 - inv_data['translation']
# ^ That's consistent with comments in reformulators.py, but doesn't produce
# correct results for tests.
#TODO: figure out what's going on, and update in reformulators.py.
x = np.hstack([x, 0])
x = np.hstack([x0[:inv_data['n']], 0])
# The final coordinate is a dummy value, which is loaded into ScalarVariables
# which (1) did not participate in the problem, but (2) whose parent Variable
# object did participate in the problem.
Expand Down
182 changes: 8 additions & 174 deletions sageopt/coniclifts/reformulators.py
View file
Expand Up @@ -36,156 +36,6 @@ def build_cone_type_selectors(K):
return type_selectors


def find_cones_with_mutually_disjoint_scope(cones):
# The runtime of this function scales quadratically in the length of the list "cones".
# If Mosek is going to be viable for compiling cone programs of nontrivial size, then
# this function needs to be orders of magnitude faster. In addition to changing how it's
# written, the final function might be simple enough to be accelerated by numba or cython.
disjoint_cone_idxs = set()
remaining_cone_idxs = {i for i in range(len(cones))}
for co_idx, co in enumerate(cones):
if 'isolated' in co.annotations and co.annotations['isolated']:
disjoint_cone_idxs.add(co_idx)
remaining_cone_idxs.remove(co_idx)
else:
co_cols = co.annotations['scope']
is_isolated = True
for other_idx in remaining_cone_idxs:
if other_idx == co_idx:
continue
other_cols = cones[other_idx].annotations['scope']
if len(co_cols.intersection(other_cols)) > 0:
is_isolated = False
break
if is_isolated:
disjoint_cone_idxs.add(co_idx)
remaining_cone_idxs.remove(co_idx)
disjoint_cones = []
for i in disjoint_cone_idxs:
disjoint_cones.append(cones[i])
return disjoint_cones


def find_diagonal_cone_constraints(A, K):
A = sp.csr_matrix(A)
diag_cones = []
for co in K:
co_start, co_stop = co.annotations['ss']
# Check if this cone is diagonal.
curr_A = A[co_start:co_stop, :].tocoo()
if 'diagonal' in co.annotations and co.annotations['diagonal']:
co.annotations['A'] = curr_A
diag_cones.append(co)
elif util.sparse_matrix_is_diag(curr_A):
co.annotations['A'] = curr_A
diag_cones.append(co)
return diag_cones


def prep_sep_dis_diag_cone_cons(A, b, K, dont_sep=None):
"""
"Prepare to separate disjoint diagonal cone constraints (from the conic system (A, b, K))."
K is a list a Cone objects.
"""
if dont_sep is None:
dont_sep = set()
dont_sep = {'0'}.union(dont_sep)
#
# Before doing anything we annotate the cones with the information necessary
# to place themselves relative to (A,b). Then we find the disjoint diagonal cones.
#
Cone.annotate_cone_positions(K) # gives cones 'ss' and 'position' annotations
relevant_cones = [co for co in K if co.type not in dont_sep]
diag_cones = find_diagonal_cone_constraints(A, relevant_cones) # gives cones an 'A' annotation
Cone.annotate_cone_scopes(A, diag_cones) # gives cones a 'scope' annotation.
disjoint_diag_cones = find_cones_with_mutually_disjoint_scope(diag_cones)
disjoint_diagonal_idxs = {co.annotations['position'] for co in disjoint_diag_cones}
#
# Build the return values
#
K0 = [co for i, co in enumerate(K) if i not in disjoint_diagonal_idxs]
for i in range(len(K0)):
K0[i].annotations = dict()
sep_K0 = []
type_selectors = build_cone_type_selectors(K)
scalings = np.ones(A.shape[1])
translates = np.zeros(A.shape[1])
for co in disjoint_diag_cones:
co_start, co_stop = co.annotations['ss']
# record normalization and translation necessary to homogenize these cones
curr_A = co.annotations['A'] # curr_A is COO format.
nz_vals, nz_rows, nz_cols = curr_A.data, curr_A.row, curr_A.col
for i, col_idx in enumerate(nz_cols):
row_idx = co_start + nz_rows[i]
translates[col_idx] -= b[row_idx]
scalings[col_idx] /= nz_vals[i]
# record a map from coordinate indices in this this cone to the columns
# (variable indices) that must belong to the cone at the coordinate index.
# i.e., sort the entries of nz_cols according to the entries of nz_rows
type_selectors[co.type][co_start:co_stop] = False
col_mapping = nz_cols[np.argsort(nz_rows)]
new_co = Cone(co.type, co.len, {'col mapping': col_mapping})
sep_K0.append(new_co)
# Build the final selector
all_selectors = [sel.reshape(-1, 1) for sel in type_selectors.values()]
all_selectors = np.hstack(all_selectors)
row_selector = np.any(all_selectors, axis=1)
return row_selector, K0, sep_K0, scalings, translates


def separate_disjoint_diagonal_cones(A, b, K, destructive=False, dont_sep=None):
"""
This function's purpose ...
coniclifts' standard form for a feasible set is {x : A @ x + b \in K}. Some solvers
require alternative standard forms, such as {x : A @ x == b, x \in K },
or {x : A @ x <= b, x \in K }. Mosek requires an even more peculiar form, namely
{(x,X) : A @ x + F(X) <= b, x \in K, X >> 0 }, where F is a particular linear map.
We can trivially convert from the coniclifts standard form to the above forms by
introduction of slack variables. I.e. represent {x : A @ x + b \in K} as
{ x : y == A @ x + b, y \in K }. The drawback is that doing this substitution
can substantially worsen the conditioning of the problem. This function exists
to facilitate similar reformulation, without introducing unnecessary slack
variables.
:param A: a scipy csc sparse matrix (m rows, n columns)
:param b: a numpy 1darray of length m
:param K: a list of cones, specified by (str, int) pairs, where "str" gives the
cones type, and "int" gives the cone's length.
:param destructive: a boolean. Indicates whether or not inputs (A,b) are modified,
or copied.
:param dont_sep: a set of cone types that are not to be modified while reformulating
the system. Common values are dont_sep={'0'} and dont_sep={'0','+'}.
:return: See in-line comments.
"""
if not destructive:
A = A.copy()
b = b.copy()
sel, K0, sep_K0, scale, trans = prep_sep_dis_diag_cone_cons(A, b, K, dont_sep=dont_sep)
# In order to separate disjoint diagonal cones from the system (A, b, K) as suggested by
# the function above, we will need to apply a diagonal affine change of coordinates to
# {x : A @ x + b \in K}. That diagonal change of coordinates maps "x" to
# x0 == np.diag(scale) @ (x + trans).
#
# Since "x" and "x0" don't yet exist on the computer, we transform (A, b) to
# (A0, b0) in such a way that the set of feasible "x" maps to {x0 : A0 @ x0 + b0 \in K}.
b0 = b + A.dot(trans)
A0 = A @ sp.diags(scale)
# Now that this map has been applied, we can drop a subset of rows from (A, b)
# --- they are now accounted for by "sep_K0". The returned value "K0" has
# already been updated to reflect the fact that we intend to drop these rows.
A0 = A0[sel, :]
b0 = b0[sel]
# Our feasible set has now been mapped to { x0 : A0 @ x0 + b0 \in K0 and x0 \in sep_K0 }.
#
# In order to invert this map later, we return both (A0, b0, K0, sep_K0) and the
# transformation data given by (scale, trans).
return A0, b0, K0, sep_K0, scale, trans


def replace_nonzero_cones_with_zero_cones_and_slacks(A, K, destructive=False, dont_rep=None):
# It only really makes sense to call this function after any isolated diagonal cone
# constraints have been separated / factored out of (A,K).
Expand Down Expand Up @@ -220,7 +70,7 @@ def replace_nonzero_cones_with_zero_cones_and_slacks(A, K, destructive=False, do
return A, K, slacks_K


def separate_cone_constraints(A, b, K, destructive=False, dont_sep=None, avoid_slacks=False):
def separate_cone_constraints(A, b, K, destructive=False, dont_sep=None):
"""
Replace the conic system {x : A @ x + b \in K } by {x0 : A1 @ x1 + b1 \in K}
Expand All @@ -232,8 +82,6 @@ def separate_cone_constraints(A, b, K, destructive=False, dont_sep=None, avoid_s
or copied.
:param dont_sep: a set of cone types that are not to be modified while reformulating
the system. Common values are dont_sep={'0'} and dont_sep={'0','+'}.
:param avoid_slacks: boolean. Indicates whether or not to perform pre-processing to minimize
the number of slack variables.
:return: See in-line comments.
Expand Down Expand Up @@ -280,32 +128,18 @@ def separate_cone_constraints(A, b, K, destructive=False, dont_sep=None, avoid_s
K = K.copy()
if dont_sep is None:
dont_sep = set('0')
if avoid_slacks:
# (1) identify "simple" nonlinear cones, and factor them out from the matrix part of the conic system.
A0, b0, K0, sep_K0, scale, trans = separate_disjoint_diagonal_cones(A, b, K, True, dont_sep)
# A0 has the same number of columns as A, but likely fewer rows.
# K0 is no longer than K.
# sep_K0 contains information for block elementwise constraints on "x".
# scale and trans contain the data for a diagonal affine change of coordinates between our
# original variables {x : A @ x + b \in K} and the new variables {x0 : A0 @ x0 + b0 \in K0, x0 \in sep_K0 }.
else:
A0 = A
K0 = K
b0 = b
sep_K0 = []
scale = np.ones(A.shape[1])
trans = np.zeros(A.shape[1])
# (2) introduce slack variables (and separately record constraints) to factor remaining constraints
# of type "dont_sep" out of the matrix and into disjoint diagonal homogeneous cone constraints.
A0 = A
K0 = K
b0 = b
sep_K0 = []
# introduce slack variables (and separately record constraints) to factor remaining constraints
# of type "dont_sep" out of the matrix and into disjoint diagonal homogeneous cone constraints.
A1, K1, sep_K1 = replace_nonzero_cones_with_zero_cones_and_slacks(A0, K0, destructive=True, dont_rep=dont_sep)
scale = np.hstack((scale, np.ones(A1.shape[1] - len(scale))))
trans = np.hstack((trans, np.zeros(A1.shape[1] - len(trans))))
# A1 has the same number of rows of A0, but likely has more columns.
# K1 is of the same length of K0, but only contains the cones in "dont_sep".
# sep_K1 defines the constraints on any slack variables which allowed us to replace K0 by K1 as above.
# The introduction of new variables required that we appropriately pad "scale" and "trans".
sep_K = sep_K0 + sep_K1
return A1, b0, K1, sep_K, scale, trans
return A1, b0, K1, sep_K


def separated_cones_to_matrix_cones(sep_K, num_cols, destructive=False):
Expand Down

0 comments on commit cc8421c

Please sign in to comment.