please expose the extreme ray code

[BrannonKing]

I've been studying nlhdlr_quadratic.c to understand how to iterate the extreme rays when at an optimum of the LP relaxation. It's more complex than expected. (It might be less complex if the code for handling quadratics was not part of it.) Is it possible to expose some public methods for iterating the extreme rays of an optimal corner? This is surely a common need for modern cutting techniques, no?

In the meantime, can you give some guidance on iterating extreme rays of the corner in a custom separator for MILP problems?

[svigerske]

CC @antoniach @fserra

[BrannonKing]

I tried this again today. Calling getLPBInvRow I get the correct magnitudes. However, I don't get the correct sign on those numbers. I don't know how to get the correct sign on those numbers. They are inconsistent in some way that I cannot detect from inspection. As I understood it from the documentation, I need to use getLPBInvARow to get the extreme rays of the cone. That data does not have the right magnitude or sign; it's all zeros generally.

When I look at the GMI code, I see that it changes the sign on rows that have status "lower". However, of the numbers I'm using, none of them come from "lower" rows. From my data, it seems that I often need to negate the second component (column). However, I can't see anything that makes that column different from any others, looking at all the column's properties -- they're no different than the first column.

[antoniach]

Hi Brannon,

please excuse the late reply and thank you for your patience!

In nlhdlr_quadratic.c, the rays are created by calling createAndStoreSparseRays(). I think the documentation of that function (lines 809-844) gives you a pretty good overview of what you need to do to get the rays. There you also see that, unfortunately, it is not enough to just call SCIPgetLPBInvARow(). What the code creating the rays basically does is first build dense rays and then store them in a sparser format so only nonzero entries are stored. To make sure that we get the right entries for all variables, you will see that there are some helper-functions that map variable indices to their position in the LP, etc. I think to understand the rough idea of the code, this can be ignored for now. Furthermore, let me draw your attention to some key aspects that will hopefully help you understand the code:

1. When we build the rays, we only consider the variables that appear in the constraint we want to enforce. Since we consider quadratic constraints in nlhdlr_quadratic.c, you will see that we often iterate over the variables appearing there (e.g., in line 520). Depending on which rays you want to access, this needs to be adapted.
2. In SCIP, a variable x can have upper and lower bounds (and not just x >= 0). Since we want to store the rays as if every nonbasic variable was at its lower bound (so that all rays move to infinity), we need to change the sign of the entries corresponding to variables at its upper bound (you can see that in lines 910-916).
3. Since SCIP introduces slack variables, we also need to consider those when building the extreme rays. For this, we have to call SCIPgetLPBInvRow() to get the right entries. The entries of the basic non-slack and slack variables are set in function storeDenseTableauRow().

I hope this helps for now. I agree with you that SCIP could benefit from a more accessible method to get the rays. That is something on our ToDo-list. Let me know if you have any other questions.

Best, Antonia

[BrannonKing]

I think I have it working with the code below (translated from nlhdlr), in some simple tests. Please tell if anything looks incorrect. The nested loops and frequent calls to getLPBInvARow get expensive, so I will be more selective as to when I call them in my real program.

        # in sepaexeclp():

        def scale_by_status(st, val):
            if st == 'upper':
                return -1.0*val
            if st == 'lower':
                return val
            if st == 'zero':
                return 0
            raise NotImplemented()

        tab_to_col = []
        tableau = np.zeros((sum(1 for x in bases if x >= 0), len(cols)))
        bases = m.getLPBasisInd()
        cols = m.getLPColsData()
        rows = m.getLPRowsData()

        for r1, row in enumerate(rows):
            c = bases[r1]
            if c < 0:
                continue
            invRow = m.getLPBInvRow(r1)
            invARow = m.getLPBInvARow(r1)
            ray = 0
            for c2, col2 in enumerate(cols):
                status = col2.getBasisStatus()
                if status != 'basic':
                    tableau[len(tab_to_col), ray] = scale_by_status(status, invARow[c2])
                    ray += 1
            for r2, row2 in enumerate(rows):
                status = row2.getBasisStatus()
                if status != 'basic':
                    tableau[len(tab_to_col), ray] = scale_by_status(status, invRow[r2])
                    ray += 1
            tab_to_col.append(c)

[BrannonKing]

I'm reopening this. I'm not convinced that the above code returns the correct tableau. Also, it shouldn't be necessary to do nested enumeration on the rows to get this right. Maybe I need to ignore the "upper" status when both r1 and r2 have it?

[BrannonKing]

There's the example in the Conforti "Integer Programming" book, page 238. Also, as a simple example: max {y : .9x - .9y >= 1, .9x + .6y <= 2.5}. That should give rays (.44, -.66), (-.66, -.66). Turn off presolve, heuristics, and propagation for testing these.

[BrannonKing]

I realized that the documentation on corner rays requires a standard form LP. By having mixed inequality directions, I mess that up. Does SCIP have a way to force or convert the model to standard form before running? Does SCIP give me a way to get the sense of inequality direction for a slack-variable column? Also, does SCIP give me a way to verify that a point is within the LP relaxation, for purposes of testing my rays?

[antoniach]

Hi Brannon,

you are right, SCIP does not transform the LP into standard form. Every row and column has lower and upper bounds (which might be set to infinity). You can access those with the functions SCIPcolGetLb(), SCIPcolGetUb(), SCIProwGetLb() and SCIProwGetUb(). Remember that slack variables are not explicitly introduced in SCIP (so there is no column attached to a slack variable), but only implicitly introduced (so every slack variable corresponds to a row). So if you want to find out the inequality direction of a slack variable, you have to check the bounds of the corresponding row by accessing SCIProwGetLb() and SCIProwGetUb(). To check whether a non-basic slack variable is at its upper or lower bound, you need to call SCIProwGetBasisStatus(). Note though that the basic status of rows are flipped compared to columns (see lines 956 - 963 of nlhdlr_quadratic.c).

Regarding the code you sent earlier, there are at least two other things you need to watch out for:

1. you flip the value wrong when checking the basis status of a column (see lines 910-918 of nlhdlr_quadratic.c and lines 809-844 for an explanation why you need to do it like that).
2. you need to consider that columns and rows have different types of indices. For example, a column has a lppos (column position number in current LP) which does not necessarily correspond to the index the column has when getting it through SCIPgetLPCols(). So you need to be aware of that. That is why the code in nlhdlr_quadratic.c defines so many different maps between different types of indices.

I hope this helps.

Best, Antonia

[BrannonKing]

@antoniach , can you further explain your thinking on this statement quoted below? Rays are scalable by definition, no? You just multiply by a scalar to change their size. Are you implying something about the feasibility of a scalar? E.g., if I scale my ray by some small value, I should still be feasible? Or how would I go about testing that my rays are correct? I was thinking to just multiply them by a small value (0.0001) and add them to my current primsol and then test them with the original constraints. This usually works but it doesn't always work: I struggle with constraints that had negative b to begin with. So is there a transformed constraint model that I also need to account for in testing my rays?

In SCIP, a variable x can have upper and lower bounds (and not just x >= 0). Since we want to store the rays as if every nonbasic variable was at its lower bound (so that all rays move to infinity), we need to change the sign of the entries corresponding to variables at its upper bound (you can see that in lines 910-916)

[antoniach]

I am not sure, what you mean with "Rays are scalable by definition". The goal of getting the rays in nlhdlr_quadratic.c is to get the corner polyhedron at the current basic solution. I.e., we want to find the rays r_1, ..., r_n such that we can represent the variables x in terms of the non-basic variables s_1, ..., s_n: x = f + r_1 s_1 + ... + r_n s_n. So you cannot just scale the rays without changing this constraint. The statement you quoted refers to making sure that we get the rays in a "standardized" form, meaning that we can always move in a positive direction along the rays without immediately leaving the feasible region of the LP relaxation. (So we don't really "need" to change the sign of the rays, as I mentioned in the quote. We just choose to, to make things easier when working with the rays later on).

"I was thinking to just multiply them by a small value (0.0001) and add them to my current primsol and then test them with the original constraints" -- If the resulting point is feasible for the LP relaxation, depends on the problem.

I hope this helps.