search

Foggalong / RobustOCS

Robust optimal contirbution selection problems for genetics with Python
https://pypi.org/project/robustocs
MIT License
1 stars 0 forks source link

SQP Method with HiGHS #10

Closed Foggalong closed 2 months ago

Foggalong commented 3 months ago

An SQP method using Gurobi was added as a proof of concept while I was working out the HiGHS documentation. Now clarified that this is the HiGHS' reference QP example (with this using passModel as an alternative), we will seek to implement this method in HiGHS instead.

Initial draft of the PR just templates the functions, nothing more. Aside from implementing the method, it'll be necessary to write some new loaders to get HiGHS what it needs. The key steps to completion (in no particular order) are:

Foggalong commented 3 months ago

For checking the two methods within Gurobi, aside from verifying output from gurobi_robust_genetics_sqp was consistent with Julian's output (it was), I used the following internal comparison:

import alphargs

# key problem variables loaded from standard format txt files
sigma, mubar, omega, n = alphargs.load_problem(
    "examples/04/A04.txt",
    "examples/04/EBV04.txt",
    "examples/04/S04.txt"
)

sires = range(0, n, 2)
dams = range(1, n, 2)

lam = 0.5
kap = 1

# robust direct optimization using Gurobi
w_rbs, z_rbs, obj_rbs = alphargs.gurobi_robust_genetics(
    sigma, mubar, omega, sires, dams, lam, kap, n)
# robust SQP optimization using Gurobi
w_sqp, z_sqp, obj_sqp = alphargs.gurobi_robust_genetics_sqp(
    sigma, mubar, omega, sires, dams, lam, kap, n)

# compare robust using direct and SQP
alphargs.print_compare_solutions(
    w_rbs, w_sqp, obj_rbs, obj_sqp, z1=z_rbs, z2=z_sqp,
    name1="w_rbs", name2="w_sqp")

print("\nDirect:")
if not alphargs.check_uncertainty_constraint(z_rbs, w_rbs, omega, debug=True):
    raise ValueError

print("\nSQP:")
if not alphargs.check_uncertainty_constraint(z_sqp, w_sqp, omega, debug=True):
    raise ValueError

print("\nDone!")

which produces the output

i  w_rbs    w_sqp
1  0.38200  0.38214
2  0.38200  0.38214
3  0.11800  0.11786
4  0.11800  0.11786

w_rbs objective: 0.77684 (z = 0.37924)
w_sqp objective: 0.77684 (z = 0.37892)
Maximum change: 0.00014
Average change: 0.00014
Minimum change: 0.00014

Direct:
     z: 0.37923871642022844
w'*Ω*w: 0.3792386953366983
  Diff: 2.1083530143961582e-08

SQP:
     z: 0.37892405482420255
w'*Ω*w: 0.37892405801748014
  Diff: -3.1932775867993257e-09

I'm not sure why the slight (but not insignificant) difference between direct robust optimization and using SQP, though I haven't yet checked the MPS files to see if there's a model error or similar.

Foggalong commented 3 months ago

Using a similar test to above but just looking at the standard non-robust problem, 

import alphargs

sigma, mubar, omega, n = alphargs.load_problem(
    "examples/50/A50.txt",
    "examples/50/EBV50.txt",
    "examples/50/S50.txt",
    issparse=True  # loads into SciPy, not NumPy
)

sires = range(0, n, 2)
dams = range(1, n, 2)

lam = 0.5
kap = 1

w_gurobi, obj_gurobi = alphargs.gurobi_standard_genetics(
    sigma, mubar, sires, dams, lam, n)

w_highs, obj_highs = alphargs.highs_standard_genetics(
    sigma, mubar, sires, dams, lam, n)

alphargs.print_compare_solutions(
    w_gurobi, w_highs, obj_gurobi, obj_highs, name1="Gurobi", name2="HiGHS")

gives weights that differ by order $10^{-5}$ between Gurobi and HiGHS, the largest being a difference of $0.00014$. Tried to output model files to compare, but output through highspy didn't seem to be working. Unsure if this is an error somewhere in my order of operations or a bug elsewhere.

jajhall commented 3 months ago

Tried to output model files to compare, but output through highspy didn't seem to be working. Unsure if this is an error somewhere in my order of operations or a bug elsewhere.

model file can be used externally for verification # BUG not working

if debug: h.setOptionValue('write_model_to_file', True) h.setOptionValue('write_model_file', 'standard-opt.mps')

These option values force writing of the model to a file only when running HiGHS from the command line - allows conversion from .mps to .lp. When using the callable library (or highspy) you need to call h.writeModel("foo.mps")

Foggalong commented 2 months ago

I say "working" in ba6d088 because the HiGHS SQP function gives the correct answer for the $n = 4$ example, but is 33 times slower than Gurobi. It also core dumps for the $n = 50$ on the first call of h.run(), never even making it to the first new constraint being added.

Between the two, it points to some issue with how I've defined the $z$ variable as an extension of $\mathbf{w}$. There may be issues with the constraint definition too, but that $n=50$ doesn't even get there shows it's not the immediate problem.

Foggalong commented 2 months ago

@jajhall Once h.passModel(...) has been used, is it still possible to do something like

h.addVar(0, highspy.kHighsInf)
h.changeColCost(dimension, kappa)

to add an additional variable? I tried doing this to incorporate $z$ as an alternative to what was pushed in ba6d088, but that just resulted in 'Model error' without any of the usual hints about what it was that caused the error.

jajhall commented 2 months ago

I've checked out https://github.com/Foggalong/robust-genetics/commit/ba6d0883f7f85c945756bf43f807106480dd25aa, but python3 example.py seems just to run gurobi.

What version of highspy are you running? There was a bug with the naming of addvar - in particular following it with changeColCost that's fixed in v1.7.1

Foggalong commented 2 months ago

What version of highspy are you running?

According to pip show highspy it's version 1.7.1.dev2.

I've checked out https://github.com/Foggalong/robust-genetics/commit/ba6d0883f7f85c945756bf43f807106480dd25aa, but python3 example.py seems just to run gurobi.

Yeah, of the back of what you said earlier I've been trying to avoid committing changes to the example file to avoid flooding the history with minor changes while debugging. I've tried to make usage of the two similar, so it just replaces gurobi_ with highs_ and adds issparse=True to the loader if it wasn't already there.

This is what I'm using to compare HiGHS-HiGHS:

sigma, mubar, omega, n = alphargs.load_problem(
    "examples/04/A04.txt",
    "examples/04/EBV04.txt",
    "examples/04/S04.txt",
    issparse=True
)

sires = range(0, n, 2)
dams = range(1, n, 2)

lam = 0.5
kap = 1

# computes the standard and robust genetic selection solutions
w_std, obj_std = alphargs.highs_standard_genetics(
    sigma, mubar, sires, dams, lam, n)
w_rbs, z_rbs, obj_rbs = alphargs.highs_robust_genetics_sqp(
    sigma, mubar, omega, sires, dams, lam, kap, n)

alphargs.print_compare_solutions(
    w_std, w_rbs, obj_std, obj_rbs, z2=z_rbs, name1="w_std", name2="w_rbs")

And this is what I'm using to compare Gurobi-HiGHS:

sigma, mubar, omega, n = alphargs.load_problem(
    "examples/50/A50.txt",
    "examples/50/EBV50.txt",
    "examples/50/S50.txt",
    issparse=True
)

sires = range(0, n, 2)
dams = range(1, n, 2)

lam = 0.5
kap = 1

t0 = time.time()
w_grb, z_grb, obj_grb = alphargs.gurobi_robust_genetics_sqp(
    sigma, mubar, omega, sires, dams, lam, kap, n, debug=False)
t1 = time.time()
print(f"Gurobi took {t1-t0:.5f} seconds")

t0 = time.time()
w_his, z_his, obj_his = alphargs.highs_robust_genetics_sqp(
    sigma, mubar, omega, sires, dams, lam, kap, n, debug=False)
t1 = time.time()
print(f"HiGHS took  {t1-t0:.5f} seconds")

alphargs.print_compare_solutions(
    w_grb, w_his, obj_grb, obj_his, z1=z_grb, z2=z_his,
    name1="Gurobi", name2="HiGHS", tol=1e-6
)
Foggalong commented 2 months ago

I'm not certain yet, but I think the issue for the pushed method comes from not extending sigma correctly when z is incorporated. Think I may have jumbled the indexes and pointers in those definitions.

jajhall commented 2 months ago

Commenting out the

if not debug: h.setOptionValue('output_flag', False) h.setOptionValue('log_to_console', False)

I immediately see reported

ERROR: Matrix dimension validation fails on index size = 50 < 51 = number of nonzeros ERROR: LP dimension validation (passModel) fails on a_matrix dimensions ERROR: LP dimension validation (passModel) fails h.passModel(model) returns status = HighsStatus.kError

Your model.lp_.a_matrix_.start_ contains values [0, 25, 51] so HiGHS expects model.lp_.a_matrix_.index_ to be of size 51, but it's of size 50

jajhall commented 2 months ago

Correcting the model.lp_.a_matrix_.start_ error and running again, I get

Running HiGHS 1.7.1 (git hash: 0c240d8): Copyright (c) 2024 HiGHS under MIT licence terms ERROR: Hessian matrix packed vector 51 has illegal start of 2500 > 51 = number of nonzeros h.passModel(model) returns status = HighsStatus.kError

That's because model.hessian_.start_ should be

= np.append(sigma.indptr, dimension*dimension)

not

= np.append(sigma.indptr, dimension + 1)

Correcting that, I see that HiGHS runs for max_iterations because reduction of the gap stalls at

gap = 5.2256630966862616e-08 ; robust_gap_tol = 1e-08

If I test against 10*robust_gap_tol then the output is

python3 example.py Set parameter Username Academic license - for non-commercial use only - expires 2024-12-21 Gurobi took 0.07078 seconds

dams = 25

sires = 25

dimension = 50

start = 0 25 50

hessian_start = [0, 50, 100, 150, 200, 250, 300, 350, 400, 450, 500, 550, 600, 650, 700, 750, 800, 850, 900, 950, 1000, 1050, 1100, 1150, 1200, 1250, 1300, 1350, 1400, 1450, 1500, 1550, 1600, 1650, 1700, 1750, 1800, 1850, 1900, 1950, 2000, 2050, 2100, 2150, 2200, 2250, 2300, 2350, 2400, 2450, 2500, 2500]

h.passModel(model) returns status = HighsStatus.kOk Iter 0 : gap = 0.18082919748787718 ; robust_gap_tol = 1e-08 Iter 1 : gap = 0.12806381581161647 ; robust_gap_tol = 1e-08 Iter 2 : gap = 0.032786471129276046 ; robust_gap_tol = 1e-08 Iter 3 : gap = 0.021426794862761944 ; robust_gap_tol = 1e-08 Iter 4 : gap = 0.00900572532136007 ; robust_gap_tol = 1e-08 Iter 5 : gap = 0.004463685289424507 ; robust_gap_tol = 1e-08 Iter 6 : gap = 0.0011589246184516289 ; robust_gap_tol = 1e-08 Iter 7 : gap = 0.00030844113417907715 ; robust_gap_tol = 1e-08 Iter 8 : gap = 0.00023256718986991443 ; robust_gap_tol = 1e-08 Iter 9 : gap = 5.8201795532508704e-05 ; robust_gap_tol = 1e-08 Iter 10 : gap = 1.4559719441559205e-05 ; robust_gap_tol = 1e-08 Iter 11 : gap = 3.8792773184304075e-06 ; robust_gap_tol = 1e-08 Iter 12 : gap = 1.036192015219095e-06 ; robust_gap_tol = 1e-08 Iter 13 : gap = 4.596713923388229e-07 ; robust_gap_tol = 1e-08 Iter 14 : gap = 2.650064735709723e-07 ; robust_gap_tol = 1e-08 Iter 15 : gap = 1.215385057318219e-07 ; robust_gap_tol = 1e-08 Iter 16 : gap = 6.85891470286748e-08 ; robust_gap_tol = 1e-08 HiGHS took 0.01694 seconds i Gurobi HiGHS 02 0.07811 0.07812 03 0.12466 0.12470 08 0.08075 0.08069 10 0.03534 0.03530 15 0.10591 0.10595 16 0.06712 0.06712 25 0.03404 0.03394 32 0.04460 0.04455 33 0.09334 0.09337 36 0.01040 0.01043 37 0.14205 0.14204 38 0.07912 0.07918 42 0.08135 0.08136 48 0.02321 0.02324

Gurobi objective: 1.92338 (z = 0.15064) HiGHS objective: 1.92338 (z = 0.15065) Maximum change: 0.00010 Average change: 0.00001 Minimum change: 0.00000

Done!

My instinct is that, once the gap gets to 5.2256630966862616e-08, the QP solver returns the same point, so each new linearization is identical. Hard to tell, as I can't debug when running from Python.

It would be good to know why you get a segfault in h.run() after h.passModel(model) returns an error. It shouldn't happen. Again, it's hard to debug using Python, and I don't have time to try to reproduce it in C++.

The good news is that Gurobi took 0.07078 seconds; HiGHS took 0.01694 seconds :-)

jajhall commented 2 months ago

A take-home from this is that, when you're developing the coding of model-building, keep the logging on, and test the return code from h.passModel(model), although I know that the documentation of highspy is very limited!

Foggalong commented 2 months ago

Your model.lp_.a_matrix_.start_ contains values [0, 25, 51] so HiGHS expects model.lp_.a_matrix_.index_ to be of size 51, but it's of size 50. model.hessian_.start_ should be np.append(sigma.indptr, dimension*dimension) not np.append(sigma.indptr, dimension + 1)

Ach damn! I'd suspected an issue with how I'd extended it but couldn't find which values were wrong before leaving. Thanks!

when you're developing the coding of model-building, keep the logging on

Is there a way to tell HiGHS what level of logging to provide? It would be useful to print warnings and errors, but skip everything else. Since it's solving a problem at each iteration there's a lot of terminal output which can obscure the important information

Foggalong commented 2 months ago

test the return code from h.passModel(model)

Useful to know that provides a return code, will incorporate that! Are these the values it returns? Tried to find the docs for passModel itself but couldn't find them (for any language, not just Python)

Foggalong commented 2 months ago

It would be good to know why you get a segfault in h.run() after h.passModel(model) returns an error. It shouldn't happen. Again, it's hard to debug using Python, and I don't have time to try to reproduce it in C++.

I've kept a minimum example of this error to one side in case you want it later. Specifically, in that situation where the matrix pointers/indices are setup wrong and it's non-diagonal, there's a segfault if h.writeModel(...) is called

Running HiGHS 1.7.0 (git hash: 1ddc6c347): Copyright (c) 2024 HiGHS under MIT licence terms
ERROR:   Matrix dimension validation fails on index size = 4 < 5 = number of nonzeros
ERROR:   LP dimension validation (passModel) fails on a_matrix dimensions
ERROR:   LP dimension validation (passModel) fails
Segmentation fault (core dumped)

and an aborted core dump if h.run() is called

Running HiGHS 1.7.0 (git hash: 1ddc6c347): Copyright (c) 2024 HiGHS under MIT licence terms
ERROR:   Matrix dimension validation fails on index size = 4 < 5 = number of nonzeros
ERROR:   LP dimension validation (passModel) fails on a_matrix dimensions
ERROR:   LP dimension validation (passModel) fails
Coefficient ranges:
  Matrix [1e+00, 1e+00]
  Cost   [1e+00, 2e+00]
  Bound  [1e+00, 1e+00]
  RHS    [5e-01, 5e-01]
double free or corruption (out)
Aborted (core dumped)
jajhall commented 2 months ago

I've kept a minimum example of this error to one side in case you want it later. Specifically, in that situation where the matrix pointers/indices are setup wrong and it's non-diagonal, there's a segfault if h.writeModel(...) is called

Thanks. If passModel gives an error return, then I thought that the incumbent model remained empty, but there's some evidence that it isn't. I should be able to reproduce what you're observing without the specific example.

jajhall commented 2 months ago

Useful to know that provides a return code, will incorporate that! Are these the values it returns? Tried to find the docs for passModel itself but couldn't find them (for any language, not just Python)

The most complete documentation is for the C API

https://ergo-code.github.io/HiGHS/stable/interfaces/c_api/

(since it's the one most used by developers of interfaces to HiGHS). Unfortunately, since you can't pass Python parameters by reference, some of the highspy parameter lists are different, and the values passed by reference are returned.

The rule of thumb with HiGHS is that if the method doesn't return a structure like HighsLp or HighsBasis (in which case nothing can go wrong, as it's just returning a const reference to a well-defined internal data structure) it will return a status

https://ergo-code.github.io/HiGHS/dev/structures/enums/#HighsStatus

jajhall commented 2 months ago

Is there a way to tell HiGHS what level of logging to provide? It would be useful to print warnings and errors, but skip everything else. Since it's solving a problem at each iteration there's a lot of terminal output which can obscure the important information

No, the INFO logging can't be switched off independently - as with Gurobi. What you can do it switch logging off before run() and then switch it on afterwards. By default, logging also goes to a file, but this might not work with Python.

Foggalong commented 2 months ago

logging also goes to a file, but this might not work with Python.

Yeah that doesn't seem to be working, at least when I first tried playing around with it. These are all good to know though, thanks!

Foggalong commented 2 months ago

Regression introduced, $n = 4$ HiGHS SQP now segfaults. Using h.run() gives

free(): invalid pointer
Aborted (core dumped)

while h.writeModel(...) gives a simple Segmentation fault (core dumped). Associated error log is

WARNING: LP matrix packed vector contains 1 |values| in [0, 0] less than or equal to 1e-09: ignored
ERROR:   Matrix dimension validation fails on index size = 4 < 16 = number of nonzeros
ERROR:   Matrix dimension validation fails on value size = 4 < 16 = number of nonzeros

so likely the same problem as last time, I've just oops a CSR index somewhere. Think it's that model.hessian_.start_ should be

np.append(sigma.indptr, sigma.indptr[-1])

not dimension*dimension or dimension + 1, but just double checking now.

Update: right place, wrong correction. Hopefully about to fix