A lot of the more technical terms (basis vectors, etc.) are lost on me, so I can only partially speak to the theory in this question.
my code uses linear optimization
Ish. From lp_prob += 0 (which should have instead used setObjective) we can see that there actually isn't an optimization objective at all. This is a linear program, and qualifies as MILP, but is only constraint programming and not minimization.
Since this is for theoretical mathematics, a (shaky) argument could be made that you're better off constructing your LP matrices directly instead of through PuLP. Advantages are typically faster setup time, a more bare-metal representation of the problem construction, and easier graphing of the sparsity pattern; disadvantages are that it takes more code and relies more heavily on knowledge of linear algebra (especially: Kronecker products proliferate).
I want to use sympy, so I cast it to a string and have sympy then parse it
I think you can trim off several of your dependencies. I was able to hack through a solution with only scipy, numpy and sympy, which I consider more standard than symengine and xpress.
each of the y and z terms have unit coefficient (depending on the sign), but when we express it in this level of generality some of them can be the same
This is expressible by just limiting every factor to have unique variables, and allowing multiple factors.
My program expresses it with factors of exactly one y with exactly one z
This isn't strictly necessary. I demonstrate how to allow for anywhere from two variables to all variables on either side of each factor.
is there a way to speed this up, perhaps by trimming down the number of basis vectors?
Ill-defined. You haven't shown any profiling results (and I don't, either, other than measuring setup and solve times).
Included below is a more substantial example that takes a very long time on my machine.
That example is quite different from the first two, since it has cubic terms and so can't be solved with the pair-wise LP construction I showed. Oh well - perhaps that can still serve as inspiration.
Fiddly bits about the code:
Avoid global side-effects.
except; return None is not great - the exception is too broad, and anyway, better to just let the exception fall through and let the caller catch it (or perhaps re-throw a domain-specific exception), which is how Python exception handling was designed to work.
Run your code through a PEP8 linter. For one example, you don't have enough newlines between functions.
Add type hints to your function signatures.
Include a __main__ guard.
Example code
This is more general in one sense, as it accommodates for longer factor expressions, but less general in another sense as it can only handle order-2 polynomials.
import time
from typing import Iterator
import matplotlib.pyplot as plt
import numpy as np
import sympy
from scipy.optimize import milp, Bounds, LinearConstraint
import scipy.sparse as sp
'''
Factorize an expression to the form
a*(yi + yj + ... - zk - zl + ...)(ym + yn + ... - zo - zp + ...) + ...
where:
- the outer product's coefficient a is always a positive integer
- the sign of y terms is always positive
- the sign of z terms is always negative
- within one factor, the number of y is equal to the number of z
- within one factor, there must be at least one y and one z
LP Variables:
n_max_factors: non-negative outer factor integers
n_max_factors * 2 * n_vars: binary selection variables, +1 for y, -1 for z
n_max_factors * n_addends: addends per factor, left variable, and right variable
'''
n_max_factors = 4 # Estimate; will need increasing for larger problems
def parse(expr: str) -> tuple[
sympy.Expr, # The parsed expression object
list[sympy.Symbol], # All free symbols, sorted
int, # number of variables
int, # number of y variables
int, # number of z variables
int, # number of unique addends
]:
# Parse and immediately expand to SOP form
input_expr = sympy.expand(sympy.parse_expr(expr))
# Unique tuple of all symbols sorted by name and index
symbols = sorted(input_expr.free_symbols, key=lambda s: s.name)
n_vars = len(symbols) # Number of unique polynomial variables
n_y = sum(s.name.startswith('y') for s in symbols) # Length of y vector
n_z = sum(s.name.startswith('z') for s in symbols) # Length of z vector
# For each combination of the variables. The number of unique multiplicative terms.
n_addends = n_vars * n_vars
return input_expr, symbols, n_vars, n_y, n_z, n_addends
def flatten_factors(term: sympy.Expr) -> Iterator[sympy.Expr]:
"""
Flatten a multiplicative Sympy expression into an iterator of terms
"""
if isinstance(term, (sympy.Symbol, sympy.Number)):
yield term
elif isinstance(term, sympy.Pow):
yield from tuple(flatten_factors(term.args[0])) * term.args[1]
elif isinstance(term, sympy.Mul):
for arg in term.args:
yield from flatten_factors(arg)
else:
raise NotImplementedError()
def iterate_equation(
input_expr: sympy.Expr,
symbols: list[sympy.Symbol],
n_vars: int,
n_y: int,
n_z: int,
) -> tuple[
np.ndarray, # equation coefficients, square, dense
np.ndarray, # flattened upper triangular equation coefficients
np.ndarray, # signs for each variable, one dimension
np.ndarray, # matrix of signs for each variable pair
np.ndarray, # upper triangular indices, first dimension
np.ndarray, # upper triangular indices, second dimension
]:
# Equation coefficient matrix, dense, square
equation_coefficients = np.zeros((n_vars, n_vars))
for term in input_expr.args:
args = tuple(flatten_factors(term))
coef = next(
(s for s in args if isinstance(s, sympy.Number)),
1,
)
syms = sorted(
(s for s in args if isinstance(s, sympy.Symbol)),
key=lambda s: s.name,
)
if len(syms) != 2:
raise NotImplementedError()
i = symbols.index(syms[0])
j = symbols.index(syms[1])
equation_coefficients[i, j] = coef
# Flattened version of the upper triangular half of the equation matrix
equation_triu = equation_coefficients[np.triu_indices(n_vars)]
# Vector of signs over one dimension of the equation matrix
equation_sign_vector = np.concatenate((
np.ones(n_y), -np.ones(n_z),
))
# Square outer product matrix of signs, dense, square
equation_signs = np.outer(equation_sign_vector, equation_sign_vector)
trii, trij = np.triu_indices(n_vars)
return equation_coefficients, equation_triu, equation_sign_vector, equation_signs, trii, trij
def make_bounds(
equation_signs: np.ndarray,
equation_sign_vector: np.ndarray,
n_addends: int,
) -> Bounds:
# Bounds for selectors
sel_lower = np.tile(equation_sign_vector*0.5-0.5, n_max_factors*2)
sel_upper = np.tile(equation_sign_vector*0.5+0.5, n_max_factors*2)
# Bounds for addends
add_lower = np.zeros((n_max_factors, n_addends))
add_upper = np.zeros((n_max_factors, n_addends))
add_lower[:, equation_signs.ravel() == -1] = -np.inf
add_upper[:, equation_signs.ravel() == +1] = +np.inf
return Bounds(
lb=np.concatenate((
np.zeros(n_max_factors), # outer factor coefficients
sel_lower , # selectors
add_lower.ravel(),
)),
ub=np.concatenate((
np.full(shape=n_max_factors, fill_value=np.inf),
sel_upper,
add_upper.ravel(),
)),
)
def make_constraints(
equation_triu: np.ndarray,
equation_sign_vector: np.ndarray,
equation_signs: np.ndarray,
trii: np.ndarray,
trij: np.ndarray,
n_vars: int,
n_addends: int,
) -> tuple[
tuple[LinearConstraint, ...],
int,
]:
# Per factor, there must be at least two variables selected
constraint_factorcount = LinearConstraint(
A=sp.hstack((
sp.csc_array((n_max_factors*2, n_max_factors)),
sp.kron(
A=sp.eye(n_max_factors*2),
B=equation_sign_vector,
),
sp.csc_array((n_max_factors*2, n_max_factors * n_addends)),
), format='csc'),
lb=2,
)
# Per factor, the number of y and z variables must be equal
constraint_factorparity = LinearConstraint(
A=sp.hstack((
sp.csc_array((n_max_factors*2, n_max_factors)),
sp.kron(
A=sp.eye(n_max_factors*2),
B=np.ones(n_vars),
),
sp.csc_array((n_max_factors*2, n_max_factors * n_addends)),
), format='csc'),
lb=0,
ub=0,
)
# The sum of each addend over all factors is pinned at the equation value
A_equation = sp.lil_array((len(equation_triu), n_addends))
vidx = np.arange(len(equation_triu))
A_equation[vidx, trii*n_vars + trij] = 1
A_equation[vidx, trij*n_vars + trii] = 1
constraint_equation = LinearConstraint(
A=sp.hstack((
sp.csc_array((len(equation_triu), n_max_factors)),
sp.csc_array((len(equation_triu), n_max_factors*2*n_vars)),
sp.kron(
A=np.ones(n_max_factors),
B=A_equation,
),
), format='csc'),
lb=equation_triu,
ub=equation_triu,
)
# Each addend will never exceed its corresponding outer factor
# positive addend < factor -> factor - addend > 0
# negative addend > -factor -> factor + addend > 0
constraint_addends_factor = LinearConstraint(
A=sp.hstack((
sp.kron(
A=sp.eye(n_max_factors),
B=np.ones((n_addends, 1)),
),
sp.csc_array((n_max_factors * n_addends, n_max_factors*2*n_vars)),
sp.kron(
A=-sp.eye(n_max_factors),
B=sp.diags(equation_signs.ravel()),
)
), format='csc'),
lb=0,
)
'''
If either of a given selector in a pair is zero, the corresponding addend must be zero
addend == factor * sela * selb * eqsign
addend*asign <= M*asel*vsign -> M*asel*vsign - addend*asign >= 0
addend*asign <= M*bsel*vsign -> M*bsel*vsign - addend*asign >= 0
'''
M = 10*n_max_factors*np.abs(equation_triu).max() # rough estimate; must exceed all factors
A_add = sp.vstack((
sp.diags(-equation_signs.ravel()),
sp.lil_array((n_addends, n_addends)),
), format='lil')
i, j = np.indices((n_vars, n_vars))
A_add[
n_vars*n_vars + n_vars*i + j,
n_vars*j + i,
] = -equation_signs
constraint_addends_zero = LinearConstraint(
A=sp.hstack((
sp.csc_array((n_max_factors*2*n_addends, n_max_factors)),
sp.kron(
A=sp.kron(
A=sp.eye(n_max_factors * 2),
B=sp.diags(M * equation_sign_vector),
),
B=np.ones((n_vars, 1)),
),
sp.kron(
A=sp.eye(n_max_factors),
B=A_add,
),
), format='csc'),
lb=0,
)
'''
If both selectors in a pair are non-zero, the corresponding addend must equal the factor
addend*asign >= factor - M(1 - sl*sgnl) - M(1 - sr*sgnr)
2M >= factor + Msl*sgnl + Msr*sgnr - asign*addend
'''
constraint_addends_nonzero = LinearConstraint(
A=sp.hstack((
sp.kron(
A=sp.eye(n_max_factors),
B=np.ones((n_addends, 1)),
),
sp.kron(
A=sp.eye(n_max_factors),
B=sp.hstack((
sp.kron(
A=sp.diags(M * equation_sign_vector),
B=np.ones((n_vars, 1)),
),
sp.kron(
A=np.ones((n_vars, 1)),
B=sp.diags(M * equation_sign_vector),
)
)),
),
sp.kron(
A=sp.eye(n_max_factors),
B=sp.diags(-equation_signs.ravel()),
),
), format='csc'),
ub=2*M,
)
constraints = (
constraint_factorcount,
constraint_factorparity,
constraint_equation,
constraint_addends_factor,
constraint_addends_zero,
constraint_addends_nonzero,
)
n_rows = sum(
con.A.shape[0] for con in constraints
)
return constraints, n_rows
def solve(
bounds: Bounds,
constraints: tuple[LinearConstraint, ...],
n_vars: int,
n_addends: int,
) -> np.ndarray:
result = milp(
# This is constraint programming; there is no minimization objective
c=np.zeros(n_max_factors + n_max_factors*2*n_vars + n_max_factors*n_addends),
integrality=np.concatenate((
np.ones(n_max_factors), # Coefficients are integral
np.ones(n_max_factors * 2 * n_vars), # selectors are integral
np.zeros(n_max_factors * n_addends), # addends are continuous
)),
bounds=bounds,
constraints=constraints,
)
if not result.success:
raise ValueError(result.message)
print(result.message)
return result.x
def triu_fold(square: np.ndarray) -> np.ndarray:
"""
Convert a matrix so that its upper triangular is the sum of the transpose and the lower triangular is zero
"""
output = square + square.T
trii, trij = np.tril_indices(n=square.shape[0], k=-1)
output[trii, trij] = 0
output[np.diag_indices_from(output)] = np.diag(square)
return output
def process_results(
result_x: np.ndarray,
symbols: list[sympy.Symbol],
n_vars: int,
) -> tuple[
np.ndarray, # outer factor coefficients
np.ndarray, # binary selectors
np.ndarray, # term addends
np.ndarray, # term addends, sum
np.ndarray, # term addends, upper triangular
sympy.Expr, # final output expression
]:
coeffs, selectors, addends = np.split(
result_x.round().astype(int),
(n_max_factors, n_max_factors + n_max_factors*2*n_vars),
)
selectors = selectors.reshape((n_max_factors, 2, n_vars))
addends = addends.reshape((n_max_factors, n_vars, n_vars))
addends_total = addends.sum(axis=0)
addends_triu = triu_fold(addends_total)
expr = sum(
coeff
* factor_selectors[0, :].dot(symbols)
* factor_selectors[1, :].dot(symbols)
for coeff, factor_selectors in zip(coeffs, selectors)
)
return coeffs, selectors, addends, addends_total, addends_triu, expr
def print_results(
eqn: str,
coeffs: np.ndarray,
selectors: np.ndarray,
addends: np.ndarray,
addends_triu: np.ndarray,
equation_coefficients: np.ndarray,
expr: sympy.Expr,
verbose: bool = False,
) -> None:
if verbose:
print(f'Selectors {selectors.shape}:')
print(selectors.reshape((-1, addends.shape[1])))
print(f'Term addends, output {addends.shape}:')
print(addends)
print(f'Term addends, upper triangular {addends_triu.shape}:')
print(addends_triu)
print(f'Term addends, input {equation_coefficients.shape}:')
print(equation_coefficients.round().astype(int))
print('Outer factor coefficients:', coeffs)
print('Input expression:', eqn)
print('Final expression:', expr)
def test(
equation_coefficients: np.ndarray,
coeffs: np.ndarray,
selectors: np.ndarray,
addends_full: np.ndarray,
addends_triu: np.ndarray,
) -> None:
assert np.allclose(equation_coefficients, addends_triu)
addends_reco = np.zeros_like(addends_triu)
for i_factor, (coef, (lhs, rhs), addend_slice) in enumerate(zip(coeffs, selectors, addends_full)):
product = coef*np.outer(lhs, rhs)
assert np.allclose(addend_slice, product)
addends_reco += product
addends_reco = triu_fold(addends_reco)
assert np.allclose(addends_reco, addends_triu)
def eval_constraints(
constraints: tuple[LinearConstraint, ...],
result_x: np.ndarray,
) -> None:
for constraint in constraints:
evald = constraint.A @ result_x
print('Evaluated constraint vector:')
if evald.size % 10:
print(evald)
else:
print(evald.reshape((-1, 10)))
print('Lower bound:', constraint.lb)
print(evald >= constraint.lb)
print('Upper bound:', constraint.ub)
print(evald <= constraint.ub)
print()
def graph_constraints(
eqn: str,
constraints: tuple[LinearConstraint, ...],
) -> plt.Figure:
all_constraints = sp.vstack([
con.A for con in constraints
]).T.toarray() != 0
fig, ax = plt.subplots()
fig: plt.Figure
ax: plt.Axes
fig.tight_layout()
ax.imshow(all_constraints, cmap='binary')
ax.set_title('Sparsity pattern (transpose) to solve\n' + eqn)
return fig
def run(eqn: str) -> None:
start = time.perf_counter()
input_expr, symbols, n_vars, n_y, n_z, n_addends = parse(eqn)
(
equation_coefficients, equation_triu, equation_sign_vector, equation_signs, trii, trij,
) = iterate_equation(input_expr, symbols, n_vars, n_y, n_z)
bounds = make_bounds(equation_signs, equation_sign_vector, n_addends)
constraints, n_rows = make_constraints(
equation_triu, equation_sign_vector, equation_signs, trii, trij, n_vars, n_addends,
)
setup = time.perf_counter()
result = solve(bounds, constraints, n_vars, n_addends)
solved = time.perf_counter()
coeffs, selectors, addends_full, addends_total, addends_triu, expr = process_results(result, symbols, n_vars)
print(f'Problem size: {n_rows}x{result.size}')
print(f'Setup time: {(setup - start)*1e3:.1f} ms')
print(f'Solve time: {(solved - setup)*1e3:.1f} ms')
print_results(
eqn, coeffs, selectors, addends_total, addends_triu, equation_coefficients, expr,
)
print()
graph_constraints(eqn, constraints)
# eval_constraints(constraints, result)
test(equation_coefficients, coeffs, selectors, addends_full, addends_triu)
def main() -> None:
# In-question first example
run(
'(y_1 - z_1)*(y_3 - z_2) '
'+ (y_1 + y_2 - z_1 - z_2)*(y_1 - z_3) '
'- (y_1 - z_2)*(y_1 - z_3)'
)
# Example hard-coded in OP script
run(
'(y_1-z_5)*(y_5-z_1) '
'+ (y_3 + y_7 - z_1 - z_2)*(y_2 + y_3 - z_7 - z_8) '
'- (y_3 - z_1)*(y_2 - z_7)'
)
# From chat
run(
'2*y_1**2 + 3*y_1*y_3 + 3*y_1*y_4 + 2*y_1*y_5 + 3*y_1*y_6 + 3*y_1*y_7 - 4*y_1*z_1 - 4*y_1*z_2 - 3*y_1*z_3 - 3*y_1*z_4 - 2*y_1*z_5 - 2*y_1*z_6 + y_3**2 + 2*y_3*y_4 + y_3*y_5 + 2*y_3*y_6 + 2*y_3*y_7 - 3*y_3*z_1 - 3*y_3*z_2 - 2*y_3*z_3 - 2*y_3*z_4 - y_3*z_5 - y_3*z_6 + y_4**2 + y_4*y_5 +'
'2*y_4*y_6 + 2*y_4*y_7 - 3*y_4*z_1 - 3*y_4*z_2 - 2*y_4*z_3 - 2*y_4*z_4 - y_4*z_5 - y_4*z_6 + y_5*y_6 + y_5*y_7 - 2*y_5*z_1 - 2*y_5*z_2 - y_5*z_3 - y_5*z_4 + y_6**2 + 2*y_6*y_7 - 3*y_6*z_1 - 3*y_6*z_2 - 2*y_6*z_3 - 2*y_6*z_4 - y_6*z_5 - y_6*z_6 + y_7**2 - 3*y_7*z_1 - 3*y_7*z_2 - 2*y_7*z_3 - 2*y_7*z_4 - y_7*z_5'
'- y_7*z_6 + 2*z_1**2 + 4*z_1*z_2 + 3*z_1*z_3 + 3*z_1*z_4 + 2*z_1*z_5 + 2*z_1*z_6 + 2*z_2**2 + 3*z_2*z_3 + 3*z_2*z_4 + 2*z_2*z_5 + 2*z_2*z_6 + z_3**2 + 2*z_3*z_4 + z_3*z_5 + z_3*z_6 + z_4**2 + z_4*z_5 + z_4*z_6'
'+' # Missing from pasted expression, inferred from context
'(y_1 - z_1)*(y_7 - z_2) + (y_1 - z_1)*(y_1 - z_3) + (y_1 - z_1)*(y_6 - z_4) + 2*(y_1 - z_1)*(y_4 - z_6) + (y_3 - z_1)*(y_1 - z_2) + (y_3 - z_1)*(y_3 - z_4) + (y_4 - z_1)*(y_7 - z_1) + (y_4 - z_1)*(y_1 - z_2) + (y_4 - z_1)*(y_6 - z_2) + (y_4 - z_1)*(y_3 - z_4) + (y_6 - z_1)*(y_3 - z_3)'
'+ (y_6 - z_1)*(y_6 - z_3) + (y_6 - z_1)*(y_5 - z_5) + (y_7 - z_1)*(y_5 - z_1) + (y_7 - z_1)*(y_1 - z_5) + (y_1 - z_2)*(y_7 - z_4) + 2*(y_1 - z_2)*(y_5 - z_5) + (y_3 - z_2)*(y_6 - z_6) + (y_4 - z_2)*(y_4 - z_4) + (y_7 - z_2)*(y_3 - z_2) + (y_7 - z_2)*(y_6 - z_2) + (y_7 - z_2)*(y_6 - z_3) + (y_7 - z_2)*(y_3 - z_4) + 2*(y_1 - z_3)*(y_3 - z_2)'
'+ (y_1 - z_3)*(y_1 - z_4) + (y_4 - z_3)*(y_5 - z_3) + (y_7 - z_3)*(y_7 - z_6) + (y_1 - z_4)*(y_6 - z_3) + (y_4 - z_4)*(y_6 - z_5) + (y_4 - z_4)*(y_7 - z_6) + (y_5 - z_4)*(y_3 - z_4) + (y_4 - z_5)*(y_3 - z_3) + (y_1 - z_6)*(y_6 - z_2)'
)
# Last example in question; has cubics so cannot be solved in pair-factor form
'''
run(
'(y_3 - z_1)*((y_1 + y_3 + y_4 + y_6 - z_1 - z_3 - z_5 - z_6)*(((y_6 - z_1)*(y_4 - z_1)*(y_1 + y_3 - z_1 - z_2) '
'+ (y_3 - z_2)*(y_1 - z_2)*(y_4 + y_6 - z_1 - z_2))*((y_4 - z_3)*(y_6 - z_3) + (y_1 - z_4)*(y_3 - z_4) '
'+ (y_4 + y_6 - z_1 - z_3)*(y_1 + y_3 - z_1 - z_4)) '
'- (y_6 - z_1)*(y_3 - z_1)*(y_1 - z_1)*(y_4 - z_1)*(y_1 + y_3 + y_4 + y_6 - z_1 - z_2 - z_3 - z_4)) '
'- (y_6 - z_1)*(y_3 - z_1)*(y_1 - z_1)*(y_4 - z_1)*((y_4 - z_3)*(y_6 - z_3) + (y_1 - z_4)*(y_3 - z_4) '
'+ (y_1 + y_3 - z_2 - z_4)*(y_4 + y_6 - z_2 - z_3)) - (((y_6 - z_1)*(y_4 - z_1)*(y_1 + y_3 - z_1 - z_3) '
'+ (y_1 - z_3)*(y_3 - z_3)*(y_4 + y_6 - z_1 - z_3))*((y_6 - z_1)*(y_4 - z_1)*(y_1 + y_3 - z_1 - z_2) '
'+ (y_3 - z_2)*(y_1 - z_2)*(y_4 + y_6 - z_1 - z_2)) '
'- (y_6 - z_1)*(y_3 - z_1)*(y_1 - z_1)*(y_4 - z_1)*((y_4 - z_2)*(y_6 - z_2) + (y_1 - z_3)*(y_3 - z_3) '
'+ (y_4 + y_6 - z_1 - z_2)*(y_1 + y_3 - z_1 - z_3)))) '
'- (-(y_6 - z_1)*(y_3 - z_1)*(y_1 - z_1)*(y_4 - z_1)*((y_4 - z_2)*(y_6 - z_2)*(y_1 + y_3 - z_2 - z_3) '
'+ (y_1 - z_3)*(y_3 - z_3)*(y_4 + y_6 - z_2 - z_3)) '
'+ (y_6 - z_1)*(y_3 - z_1)*(y_1 - z_1)*(y_4 - z_1)*(y_1 + y_3 + y_4 + y_6 - z_2 - z_3 - z_5 - z_6)*((y_4 - z_3)*(y_6 - z_3) '
'+ (y_1 - z_4)*(y_3 - z_4) + (y_1 + y_3 - z_2 - z_4)*(y_4 + y_6 - z_2 - z_3)))'
)
'''
plt.show()
if __name__ == '__main__':
main()
Optimization terminated successfully. (HiGHS Status 7: Optimal)
Problem size: 613x196
Setup time: 62.9 ms
Solve time: 33.4 ms
Input expression: (y_1 - z_1)*(y_3 - z_2) + (y_1 + y_2 - z_1 - z_2)*(y_1 - z_3) - (y_1 - z_2)*(y_1 - z_3)
Final expression: (y_1 - z_1)*(y_3 - z_2) + (y_1 - z_3)*(y_2 - z_1)
Optimization terminated successfully. (HiGHS Status 7: Optimal)
Problem size: 1671x484
Setup time: 29.5 ms
Solve time: 21.3 ms
Input expression: (y_1-z_5)*(y_5-z_1) + (y_3 + y_7 - z_1 - z_2)*(y_2 + y_3 - z_7 - z_8) - (y_3 - z_1)*(y_2 - z_7)
Final expression: (y_1 - z_5)*(y_5 - z_1) + (y_2 - z_7)*(y_7 - z_2) + (y_3 - z_1)*(y_3 - z_8) + (y_3 - z_8)*(y_7 - z_2)
Optimization terminated successfully. (HiGHS Status 7: Optimal)
Problem size: 2398x676
Setup time: 142.7 ms
Solve time: 939.3 ms
Input expression: 2*y_1**2 + 3*y_1*y_3 + 3*y_1*y_4 + 2*y_1*y_5 + 3*y_1*y_6 + 3*y_1*y_7 - 4*y_1*z_1 - 4*y_1*z_2 - 3*y_1*z_3 - 3*y_1*z_4 - 2*y_1*z_5 - 2*y_1*z_6 + y_3**2 + 2*y_3*y_4 + y_3*y_5 + 2*y_3*y_6 + 2*y_3*y_7 - 3*y_3*z_1 - 3*y_3*z_2 - 2*y_3*z_3 - 2*y_3*z_4 - y_3*z_5 - y_3*z_6 + y_4**2 + y_4*y_5 +2*y_4*y_6 + 2*y_4*y_7 - 3*y_4*z_1 - 3*y_4*z_2 - 2*y_4*z_3 - 2*y_4*z_4 - y_4*z_5 - y_4*z_6 + y_5*y_6 + y_5*y_7 - 2*y_5*z_1 - 2*y_5*z_2 - y_5*z_3 - y_5*z_4 + y_6**2 + 2*y_6*y_7 - 3*y_6*z_1 - 3*y_6*z_2 - 2*y_6*z_3 - 2*y_6*z_4 - y_6*z_5 - y_6*z_6 + y_7**2 - 3*y_7*z_1 - 3*y_7*z_2 - 2*y_7*z_3 - 2*y_7*z_4 - y_7*z_5- y_7*z_6 + 2*z_1**2 + 4*z_1*z_2 + 3*z_1*z_3 + 3*z_1*z_4 + 2*z_1*z_5 + 2*z_1*z_6 + 2*z_2**2 + 3*z_2*z_3 + 3*z_2*z_4 + 2*z_2*z_5 + 2*z_2*z_6 + z_3**2 + 2*z_3*z_4 + z_3*z_5 + z_3*z_6 + z_4**2 + z_4*z_5 + z_4*z_6+(y_1 - z_1)*(y_7 - z_2) + (y_1 - z_1)*(y_1 - z_3) + (y_1 - z_1)*(y_6 - z_4) + 2*(y_1 - z_1)*(y_4 - z_6) + (y_3 - z_1)*(y_1 - z_2) + (y_3 - z_1)*(y_3 - z_4) + (y_4 - z_1)*(y_7 - z_1) + (y_4 - z_1)*(y_1 - z_2) + (y_4 - z_1)*(y_6 - z_2) + (y_4 - z_1)*(y_3 - z_4) + (y_6 - z_1)*(y_3 - z_3)+ (y_6 - z_1)*(y_6 - z_3) + (y_6 - z_1)*(y_5 - z_5) + (y_7 - z_1)*(y_5 - z_1) + (y_7 - z_1)*(y_1 - z_5) + (y_1 - z_2)*(y_7 - z_4) + 2*(y_1 - z_2)*(y_5 - z_5) + (y_3 - z_2)*(y_6 - z_6) + (y_4 - z_2)*(y_4 - z_4) + (y_7 - z_2)*(y_3 - z_2) + (y_7 - z_2)*(y_6 - z_2) + (y_7 - z_2)*(y_6 - z_3) + (y_7 - z_2)*(y_3 - z_4) + 2*(y_1 - z_3)*(y_3 - z_2)+ (y_1 - z_3)*(y_1 - z_4) + (y_4 - z_3)*(y_5 - z_3) + (y_7 - z_3)*(y_7 - z_6) + (y_1 - z_4)*(y_6 - z_3) + (y_4 - z_4)*(y_6 - z_5) + (y_4 - z_4)*(y_7 - z_6) + (y_5 - z_4)*(y_3 - z_4) + (y_4 - z_5)*(y_3 - z_3) + (y_1 - z_6)*(y_6 - z_2)
Final expression: (2*y_1 + 2*y_7 - 2*z_1 - 2*z_2)*(y_1 + y_3 + y_4 + y_5 + y_6 + y_7 - z_1 - z_2 - z_3 - z_4 - z_5 - z_6) + (y_1 + y_3 + y_4 + y_6 - z_1 - z_2 - z_3 - z_4)*(2*y_1 + 2*y_3 + 2*y_4 + 2*y_5 + 2*y_6 + 2*y_7 - 2*z_1 - 2*z_2 - 2*z_3 - 2*z_4 - 2*z_5 - 2*z_6)
expandwhich will fully expand the result in terms of monomials with integer coefficients. The variables don't have anything substituted, these are abstract polynomials. It arbitrarily assigned indexes to the monomials and populates "monom_to_vec" which is a lookup translating the monomial to the index. Then, given an expanded polynomial, poly_to_vec will fill in those coordinates which are the integer coefficients of the monomials in a numpy array. So the value at the index in the array is the coefficient of the monomial. Is this clearer? \$\endgroup\$(y_5 - z_1)with nonnegative coefficients. So(y_5 - z_1)*(y_1-z_2)+(y_5-z_1)*(y_3-z_5)would be good, and the reason this is so hard is because you could also express the same polynomial as(y_5 - z_1)*(y_1-z_5)+(y_5-z_1)*(y_3-z_2). We're not looking for a unique expression or factorized, just positive in terms of terms like that. The assumption is that this is always possible (if it's not in a certain case, we don't care about that case). \$\endgroup\$(y_5 - z_1)*(y_1-z_2)+(y_5-z_1)*(y_3-z_5)as(y_5 - z_1)*(y_1-z_2+y_3-z_5)? \$\endgroup\$