parameterized vector with PySCIPopt

Question

I'm attempting to use PySCIPopt to solve a traditional Ax-b + constraints type of problem. I have many values of b, and I need to run the optimizer on each of them. How can I reuse the setup? Second question, what is the equivalent of norm in PySCIPopt? Or what is the right way to drive Ax-b as close to zero as possible? See the ??? marks below

import numpy as np
from pyscipopt import Model, quicksum

def make_program():
    A = ... load constant master matrix ...
    model = Model('Match_to_Master')
    x = []
    y = []
    for i in range(A.shape[1]):
        x.append(model.addVar(vtype='C', lb=0.0, ub=4.0, name='x(%s)' % i))
        y.append(model.addVar(vtype='B', name='y(%s)' % i))
        model.addCons(x[i] <= y[i]*4)
    for i in range(0, A.shape[1] - 20, 20):
       model.addCons(quicksum(y[i:i+20]) <= 1)

    #b = Parameter(A.shape[0], nonneg=True) ???
    model.setObjective(norm(A*x - b), sense='minimize') ???
    return b, x, model


def run_program(data, thresh=0.2):
    b, x, model = make_program()
    B = ... from data load matrix for analysis ...
    c = 0
    for column in B.T:
        b.value = column   ???
        model.optimize()  # reuse previous x values as starting point
        x.value[x.value < thresh] = 0.0
        for i in range(0, x.value.size - 20, 20):
            sum = np.sum(x.value[i:i+20])
            if sum > 0.2:
                print('  hit at ', str(i//20), ' for column ', str(c))
        c += 1

Answer 1

Well, I don't think that what you are solving is a traditional Ax - b type problem. The traditional type is min |Ax - b|^2. You don't need SCIP for that.

If you want to solve the program for multiple values of b, you should adapt your make_program function accordingly. You seem to be interested in integral solutions. While LPs have some warm starting capabilities for different right hand sides (using the dual simplex), integer programs don't have that capability.

Still, you can use a solution from a previous run. Use

solution = model.createSol(None)

# for each nonzero variable
model.setSolVal(solution, var, val)

model.trySol(solution)

to create a solution, set its values and add it to the model.

Regarding norm minimization, as pointed out by sascha, you can only minimize polyhedral norms, most likely the 1 / inf norm. In the former case, you can use

min y_1 + ... + y_m + z_1 + ... + z_m

s.t. A x + y - z = b

y, z >= 0

In the latter one, you can use

min d

s.t. A x + y - z = b

y, z >= 0

d >= y_i, i=1,...,m

d >= z_i i=1,...,m

In either case, this is not supported by scip out of the box, so you will have to add the corresponding variables/constraints yourself.

Edit:

Regarding the 2-norm: First, note that the minimization of the 2-norm turns the problem into a mixed-integer nonlinear program (MINLP), specifically, a Mixed Integer Quadratic Problem (MIQP). Thus, the problem becomes more difficult to solve :) I don't know if SCIP is the best fit in this case (I have heard good things about Pajarito). Still, SCIP can solve MINLPs.

To model the norm, you should note that SCIP does not support nonlinear objectives, only nonlinear constraints. Thus, the model should be

min y_1 + ... + y_m

s.t. (A x - b)^2 - y <= 0

You should be able to add arbitrary nonlinear constraints via model.addConst(...) Based on the degree of the expression, PySCIPOpt will do the right thing. In case of a degree-2 constraint, i.e.,

model.addCons(x[i]*y[i] <= 0)

this means adding a quadratic constraint.

Note that this will minimize the squared norm instead of the norm. You will get the correct solution, but the objective value will be off.

Also, I think it will help or maybe be required to compile SCIP with Ipopt support enabled. Ipopt is a solver for nonlinear programs (NLPs), which are the relaxations of MINLPs.