diff --git a/.vscode/settings.json b/.vscode/settings.json index 34791c3..1549fc7 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -8,6 +8,7 @@ "Gregor's", "gurobi", "gurobipy", + "Heckel", "hendo", "linalg", "loadtxt", diff --git a/Scratch/Timing Test/README.md b/Scratch/Timing Test/README.md new file mode 100644 index 0000000..891f917 --- /dev/null +++ b/Scratch/Timing Test/README.md @@ -0,0 +1,3 @@ +# Timing Test + +Quick timing test of a non-trivial problem to show that it's the same in Python vs Jupyter. \ No newline at end of file diff --git a/Scratch/Timing Test/timing-test.ipynb b/Scratch/Timing Test/timing-test.ipynb new file mode 100644 index 0000000..b5d0f2d --- /dev/null +++ b/Scratch/Timing Test/timing-test.ipynb @@ -0,0 +1,361 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Robust Optimization for Genetic Selection\n", + "\n", + "This notebook provides a didactic for using Python to solve the multi-objective optimization problem which arises in the context of robust genetic selection. It has been tested with Python 3.10 specifically and there are some standard packages which this depends on, imported below." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np # defines matrix structures\n", + "from qpsolvers import solve_qp # used for quadratic optimization\n", + "import gurobipy as gp # Gurobi optimization interface (1)\n", + "from gurobipy import GRB # Gurobi optimization interface (2)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Utility functions and output settings used in this notebook are defined in the two cells below." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "# want to round rather than truncate when printing\n", + "np.set_printoptions(threshold=np.inf)\n", + "\n", + "# only show numpy output to five decimal places\n", + "np.set_printoptions(formatter={'float_kind':\"{:.5f}\".format})" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Real (Simulated) Data\n", + "\n", + "Now we've looked at how these problems can be approached with Gurobi, we can try an example with realistic simulated data. Here we have data for a cohort of 50 candidates split across three data files, each containing space-separated values.\n", + "\n", + "1. In `A50.txt` the matrix $\\Sigma$ is described, with columns for $i$, $j$, and $a_{ij}$, where only the upper triangle is described since the matrix is symmetric.\n", + "2. In `EBV50.txt` the vector $\\bar{\\mu}$ is described, with only one column containing the posterior mean over 1000 samples.\n", + "3. In `S50.txt` the matrix $\\Omega$ is described, with columns for $i$, $j$, and $\\sigma_{ij}$, where only the upper triangle is described since the matrix is symmetric.\n", + "\n", + "Note that this particular problem does not contain a separate file for sex data; odd indexed candidates are sires and even indexed candidates are dams. For now, we also don't have a $l$ or $u$ bounding the possible weights.\n", + "\n", + "\n", + "The first step then is clearly going to be loading the matrices in particular from file. We create the following function to do so.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "def load_problem(A_filename, E_filename, S_filename, dimension=False):\n", + " \"\"\"\n", + " Used to load genetic selection problems into NumPy. It takes three\n", + " string inputs for filenames where Sigma, Mu, and Omega are stored,\n", + " as well as an optional integer input for problem dimension if this\n", + " is known. If it's know know, it's worked out based on E_filename.\n", + "\n", + " As output, it returns (A, E, S, n), where A and S are n-by-n NumPy\n", + " arrays, E is a length n NumPy array, and n is an integer.\n", + " \"\"\"\n", + "\n", + " def load_symmetric_matrix(filename, dimension):\n", + " \"\"\"\n", + " Since NumPy doesn't have a stock way to load matrices\n", + " stored in coordinate format format, this adds one.\n", + " \"\"\"\n", + "\n", + " matrix = np.zeros([dimension, dimension])\n", + "\n", + " with open(filename, 'r') as file:\n", + " for line in file:\n", + " i, j, entry = line.split(\" \")\n", + " # data files indexed from 1, not 0\n", + " matrix[int(i)-1, int(j)-1] = entry\n", + " matrix[int(j)-1, int(i)-1] = entry\n", + "\n", + " return matrix\n", + "\n", + "\n", + " # if dimension wasn't supplied, need to find that\n", + " if not dimension:\n", + " # get dimension from EBV, since it's the smallest file\n", + " with open(E_filename, 'r') as file:\n", + " dimension = sum(1 for _ in file)\n", + "\n", + " # EBV isn't in coordinate format so can be loaded directly\n", + " E = np.loadtxt(E_filename) \n", + " # A and S are stored by coordinates so need special loader\n", + " A = load_symmetric_matrix(A_filename, dimension)\n", + " S = load_symmetric_matrix(S_filename, dimension)\n", + "\n", + " return A, E, S, dimension" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For the given example problem we can now solve it with Gurobi using the exact same methods as before, for both the standard genetic selection problem and the robust version of the problem. The following cell does both alongside each other to accentuate the differences." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Set parameter Username\n", + "Academic license - for non-commercial use only - expires 2025-02-26\n", + "Set parameter TimeLimit to value 300\n", + "Set parameter MIPGap to value 0.009\n", + "Set parameter MIPGap to value 0.009\n", + "Gurobi Optimizer version 11.0.0 build v11.0.0rc2 (linux64 - \"Ubuntu 22.04.4 LTS\")\n", + "\n", + "CPU model: Intel(R) Core(TM) i5-8350U CPU @ 1.70GHz, instruction set [SSE2|AVX|AVX2]\n", + "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n", + "\n", + "Optimize a model with 4 rows, 50 columns and 100 nonzeros\n", + "Model fingerprint: 0x534ff85e\n", + "Model has 1275 quadratic objective terms\n", + "Coefficient statistics:\n", + " Matrix range [1e+00, 1e+00]\n", + " Objective range [7e-01, 1e+00]\n", + " QObjective range [5e-02, 1e+00]\n", + " Bounds range [0e+00, 0e+00]\n", + " RHS range [5e-01, 5e-01]\n", + "Presolve removed 2 rows and 0 columns\n", + "Presolve time: 0.01s\n", + "Presolved: 2 rows, 50 columns, 50 nonzeros\n", + "Presolved model has 1275 quadratic objective terms\n", + "Ordering time: 0.00s\n", + "\n", + "Barrier statistics:\n", + " Free vars : 49\n", + " AA' NZ : 1.274e+03\n", + " Factor NZ : 1.326e+03\n", + " Factor Ops : 4.553e+04 (less than 1 second per iteration)\n", + " Threads : 1\n", + "\n", + " Objective Residual\n", + "Iter Primal Dual Primal Dual Compl Time\n", + " 0 3.76397773e+05 -4.22619805e+05 2.50e+04 1.12e+00 9.99e+05 0s\n", + " 1 -5.01117677e+01 -9.99230833e+02 2.68e+01 1.25e-04 1.09e+03 0s\n", + " 2 -8.81632395e-01 -9.63667643e+02 2.68e-05 1.25e-10 1.93e+01 0s\n", + " 3 -8.81606636e-01 -2.03034027e+00 5.19e-09 2.41e-14 2.30e-02 0s\n", + " 4 -9.04332054e-01 -1.04180295e+00 4.73e-10 2.19e-15 2.75e-03 0s\n", + " 5 -9.59656693e-01 -1.01318293e+00 1.67e-15 8.88e-16 1.07e-03 0s\n", + " 6 -9.73793052e-01 -9.76849441e-01 2.22e-16 1.33e-15 6.11e-05 0s\n", + " 7 -9.76163966e-01 -9.76239391e-01 1.39e-15 6.66e-16 1.51e-06 0s\n", + " 8 -9.76213879e-01 -9.76217890e-01 9.91e-15 5.55e-16 8.02e-08 0s\n", + " 9 -9.76215229e-01 -9.76215592e-01 1.65e-13 5.55e-16 7.25e-09 0s\n", + " 10 -9.76215438e-01 -9.76215482e-01 7.03e-13 6.66e-16 8.86e-10 0s\n", + " 11 -9.76215473e-01 -9.76215476e-01 4.24e-13 6.66e-16 4.91e-11 0s\n", + "\n", + "Barrier solved model in 11 iterations and 0.03 seconds (0.00 work units)\n", + "Optimal objective -9.76215473e-01\n", + "\n", + "Gurobi Optimizer version 11.0.0 build v11.0.0rc2 (linux64 - \"Ubuntu 22.04.4 LTS\")\n", + "\n", + "CPU model: Intel(R) Core(TM) i5-8350U CPU @ 1.70GHz, instruction set [SSE2|AVX|AVX2]\n", + "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n", + "\n", + "Optimize a model with 5 rows, 51 columns and 101 nonzeros\n", + "Model fingerprint: 0x8c48581e\n", + "Model has 1275 quadratic objective terms\n", + "Model has 50 quadratic constraints\n", + "Coefficient statistics:\n", + " Matrix range [1e+00, 1e+00]\n", + " QMatrix range [7e-06, 1e+00]\n", + " Objective range [7e-01, 1e+00]\n", + " QObjective range [5e-02, 1e+00]\n", + " Bounds range [0e+00, 0e+00]\n", + " RHS range [5e-01, 5e-01]\n", + "Presolve removed 3 rows and 0 columns\n", + "\n", + "Continuous model is non-convex -- solving as a MIP\n", + "\n", + "Presolve removed 3 rows and 0 columns\n", + "Presolve time: 0.01s\n", + "Presolved: 2552 rows, 1327 columns, 7600 nonzeros\n", + "Presolved model has 1275 quadratic objective terms\n", + "Presolved model has 1 quadratic constraint(s)\n", + "Presolved model has 1275 bilinear constraint(s)\n", + "Variable types: 1327 continuous, 0 integer (0 binary)\n", + "Found heuristic solution: objective -0.9763584\n", + "\n", + "Root relaxation: objective -1.285097e+00, 31 iterations, 0.00 seconds (0.00 work units)\n", + "\n", + " Nodes | Current Node | Objective Bounds | Work\n", + " Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time\n", + "\n", + " 0 0 -1.28510 0 36 -0.97636 -1.28510 31.6% - 0s\n", + " 0 0 -1.13066 0 69 -0.97636 -1.13066 15.8% - 0s\n", + " 0 0 -1.05344 0 51 -0.97636 -1.05344 7.89% - 0s\n", + " 0 0 -1.01483 0 51 -0.97636 -1.01483 3.94% - 0s\n", + " 0 0 -0.99552 0 37 -0.97636 -0.99552 1.96% - 0s\n", + " 0 0 -0.99119 0 1175 -0.97636 -0.99119 1.52% - 0s\n", + " 0 0 -0.99098 0 1118 -0.97636 -0.99098 1.50% - 0s\n", + " 0 0 -0.98587 0 47 -0.97636 -0.98587 0.97% - 1s\n", + " 0 2 -0.98587 0 47 -0.97636 -0.98587 0.97% - 1s\n", + " 184 201 -0.98413 26 341 -0.97636 -0.98567 0.95% 156 5s\n", + " 434 418 -0.97843 58 533 -0.97636 -0.98567 0.95% 318 10s\n", + "H 609 566 -0.9768443 -0.98567 0.90% 296 12s\n", + "H 699 541 -0.9768479 -0.98567 0.90% 282 14s\n", + " 704 545 -0.97802 28 0 -0.97685 -0.98567 0.90% 280 15s\n", + "\n", + "Cutting planes:\n", + " RLT: 39\n", + " PSD: 143\n", + "\n", + "Explored 704 nodes (208237 simplex iterations) in 15.37 seconds (17.22 work units)\n", + "Thread count was 8 (of 8 available processors)\n", + "\n", + "Solution count 1: -0.976848 \n", + "\n", + "Optimal solution found (tolerance 9.00e-03)\n", + "Best objective -9.768443464472e-01, best bound -9.854709478274e-01, gap 0.8831%\n", + "\n", + "SIRE WEIGHTS\t\t\t DAM WEIGHTS\n", + "--------------------\t\t --------------------\n", + " i w_std w_rbs\t\t i w_std w_rbs\n", + "00 0.00000 0.00002 01 0.08319 0.08339\n", + "02 0.11649 0.11508 03 0.00000 0.00001\n", + "04 0.00000 0.00002 05 0.00000 0.00002\n", + "06 0.00000 0.00001 07 0.06817 0.06830\n", + "08 0.00000 0.00001 09 0.03762 0.03701\n", + "10 0.00000 0.00002 11 0.00000 0.00002\n", + "12 0.00000 0.00001 13 0.00355 0.00564\n", + "14 0.11542 0.11450 15 0.06966 0.06913\n", + "16 0.00000 0.00006 17 0.00000 0.00001\n", + "18 0.00000 0.00002 19 0.00000 0.00003\n", + "20 0.00000 0.00002 21 0.00001 0.00062\n", + "22 0.00000 0.00017 23 0.00000 0.00004\n", + "24 0.03537 0.03493 25 0.00000 0.00001\n", + "26 0.00000 0.00001 27 0.00000 0.00002\n", + "28 0.00000 0.00002 29 0.00000 0.00002\n", + "30 0.00000 0.00001 31 0.03580 0.03612\n", + "32 0.08447 0.08517 33 0.00000 0.00003\n", + "34 0.00000 0.00001 35 0.02849 0.02747\n", + "36 0.13489 0.13487 37 0.07446 0.07264\n", + "38 0.00000 0.00001 39 0.00000 0.00001\n", + "40 0.00000 0.00003 41 0.07312 0.07388\n", + "42 0.00000 0.00004 43 0.00000 0.00003\n", + "44 0.01336 0.01485 45 0.00000 0.00002\n", + "46 0.00000 0.00009 47 0.02593 0.02555\n", + "48 0.00000 0.00003 49 0.00000 0.00000\n" + ] + } + ], + "source": [ + "sigma, mubar, omega, n = load_problem(\n", + " \"../Example/A50.txt\",\n", + " \"../Example/EBV50.txt\",\n", + " \"../Example/S50.txt\",\n", + " 50)\n", + "\n", + "lam = 0.5\n", + "kappa = 2\n", + "\n", + "# define the M so that column i is [1;0] if i is a sire (so even) and [0;1] otherwise \n", + "M = np.zeros((2, n))\n", + "M[0, range(0,50,2)] = 1\n", + "M[1, range(1,50,2)] = 1\n", + "# define the right hand side of the constraint Mx = m\n", + "m = np.array([[0.5], [0.5]])\n", + "\n", + "# create models for standard and robust genetic selection\n", + "model_std = gp.Model(\"n50standard\")\n", + "model_rbs = gp.Model(\"n50robust\")\n", + "\n", + "# initialise w for both models, z for robust model\n", + "w_std = model_std.addMVar(shape=n, vtype=GRB.CONTINUOUS, name=\"w\") \n", + "w_rbs = model_rbs.addMVar(shape=n, vtype=GRB.CONTINUOUS, name=\"w\")\n", + "z_rbs = model_rbs.addVar(name=\"z\")\n", + "\n", + "# define the objective functions for both models\n", + "model_std.setObjective(\n", + " 0.5*w_std@(sigma@w_std) - lam*w_std.transpose()@mubar,\n", + "GRB.MINIMIZE)\n", + "\n", + "model_rbs.setObjective(\n", + " # Gurobi does offer a way to set one objective in terms of another, i.e.\n", + " # we could use `model_std.getObjective() - lam*kappa*z_rbs` to define this\n", + " # robust objective, but it results in a significant slowdown in code.\n", + " 0.5*w_rbs@(sigma@w_rbs) - lam*w_rbs.transpose()@mubar - lam*kappa*z_rbs,\n", + "GRB.MINIMIZE)\n", + "\n", + "# add sum-to-half constraints to both models\n", + "model_std.addConstr(M @ w_std == m, name=\"sum-to-half\")\n", + "model_rbs.addConstr(M @ w_std == m, name=\"sum-to-half\")\n", + "\n", + "# add quadratic uncertainty constraint to the robust model\n", + "model_rbs.addConstr(z_rbs**2 <= np.inner(w_rbs, omega@w_rbs), name=\"uncertainty\")\n", + "model_rbs.addConstr(z_rbs >= 0, name=\"z positive\")\n", + "\n", + "# since working with non-trivial size, set a time limit\n", + "time_limit = 60*5 # 5 minutes\n", + "model_std.setParam(GRB.Param.TimeLimit, time_limit)\n", + "model_std.setParam(GRB.Param.TimeLimit, time_limit)\n", + "\n", + "# for the same reason, also set a duality gap tolerance\n", + "duality_gap = 0.009\n", + "model_std.setParam('MIPGap', duality_gap)\n", + "model_rbs.setParam('MIPGap', duality_gap)\n", + "\n", + "# solve both problems with Gurobi\n", + "model_std.optimize()\n", + "model_rbs.optimize()\n", + "\n", + "# HACK code which prints the results for comparison in a nice format\n", + "print(\"\\nSIRE WEIGHTS\\t\\t\\t DAM WEIGHTS\")\n", + "print(\"-\"*20 + \"\\t\\t \" + \"-\"*20)\n", + "print(\" i w_std w_rbs\\t\\t i w_std w_rbs\")\n", + "for candidate in range(25):\n", + " print(f\"{candidate*2:02d} {w_std.X[candidate*2]:.5f} {w_rbs.X[candidate*2]:.5f} \\\n", + " {candidate*2+1:02d} {w_std.X[candidate*2+1]:.5f} {w_rbs.X[candidate*2+1]:.5f}\")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.12" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/Scratch/Timing Test/timing-test.py b/Scratch/Timing Test/timing-test.py new file mode 100644 index 0000000..053ff32 --- /dev/null +++ b/Scratch/Timing Test/timing-test.py @@ -0,0 +1,122 @@ +import numpy as np # defines matrix structures +from qpsolvers import solve_qp # used for quadratic optimization +import gurobipy as gp # Gurobi optimization interface (1) +from gurobipy import GRB # Gurobi optimization interface (2) + +# want to round rather than truncate when printing +np.set_printoptions(threshold=np.inf) + +# only show numpy output to five decimal places +np.set_printoptions(formatter={'float_kind':"{:.5f}".format}) + + +def load_problem(A_filename, E_filename, S_filename, dimension=False): + """ + Used to load genetic selection problems into NumPy. It takes three + string inputs for filenames where Sigma, Mu, and Omega are stored, + as well as an optional integer input for problem dimension if this + is known. If it's know know, it's worked out based on E_filename. + + As output, it returns (A, E, S, n), where A and S are n-by-n NumPy + arrays, E is a length n NumPy array, and n is an integer. + """ + + def load_symmetric_matrix(filename, dimension): + """ + Since NumPy doesn't have a stock way to load matrices + stored in coordinate format format, this adds one. + """ + + matrix = np.zeros([dimension, dimension]) + + with open(filename, 'r') as file: + for line in file: + i, j, entry = line.split(" ") + # data files indexed from 1, not 0 + matrix[int(i)-1, int(j)-1] = entry + matrix[int(j)-1, int(i)-1] = entry + + return matrix + + + # if dimension wasn't supplied, need to find that + if not dimension: + # get dimension from EBV, since it's the smallest file + with open(E_filename, 'r') as file: + dimension = sum(1 for _ in file) + + # EBV isn't in coordinate format so can be loaded directly + E = np.loadtxt(E_filename) + # A and S are stored by coordinates so need special loader + A = load_symmetric_matrix(A_filename, dimension) + S = load_symmetric_matrix(S_filename, dimension) + + return A, E, S, dimension + + +sigma, mubar, omega, n = load_problem( + "../Example/A50.txt", + "../Example/EBV50.txt", + "../Example/S50.txt", + 50) + +lam = 0.5 +kappa = 2 + +# define the M so that column i is [1;0] if i is a sire (so even) and [0;1] otherwise +M = np.zeros((2, n)) +M[0, range(0,50,2)] = 1 +M[1, range(1,50,2)] = 1 +# define the right hand side of the constraint Mx = m +m = np.array([[0.5], [0.5]]) + +# create models for standard and robust genetic selection +model_std = gp.Model("n50standard") +model_rbs = gp.Model("n50robust") + +# initialise w for both models, z for robust model +w_std = model_std.addMVar(shape=n, vtype=GRB.CONTINUOUS, name="w") +w_rbs = model_rbs.addMVar(shape=n, vtype=GRB.CONTINUOUS, name="w") +z_rbs = model_rbs.addVar(name="z") + +# define the objective functions for both models +model_std.setObjective( + 0.5*w_std@(sigma@w_std) - lam*w_std.transpose()@mubar, +GRB.MINIMIZE) + +model_rbs.setObjective( + # Gurobi does offer a way to set one objective in terms of another, i.e. + # we could use `model_std.getObjective() - lam*kappa*z_rbs` to define this + # robust objective, but it results in a significant slowdown in code. + 0.5*w_rbs@(sigma@w_rbs) - lam*w_rbs.transpose()@mubar - lam*kappa*z_rbs, +GRB.MINIMIZE) + +# add sum-to-half constraints to both models +model_std.addConstr(M @ w_std == m, name="sum-to-half") +model_rbs.addConstr(M @ w_std == m, name="sum-to-half") + +# add quadratic uncertainty constraint to the robust model +model_rbs.addConstr(z_rbs**2 <= np.inner(w_rbs, omega@w_rbs), name="uncertainty") +model_rbs.addConstr(z_rbs >= 0, name="z positive") + +# since working with non-trivial size, set a time limit +time_limit = 60*5 # 5 minutes +model_std.setParam(GRB.Param.TimeLimit, time_limit) +model_std.setParam(GRB.Param.TimeLimit, time_limit) + +# for the same reason, also set a duality gap tolerance +duality_gap = 0.009 +model_std.setParam('MIPGap', duality_gap) +model_rbs.setParam('MIPGap', duality_gap) + +# solve both problems with Gurobi +model_std.optimize() +model_rbs.optimize() + +# HACK code which prints the results for comparison in a nice format +print("\nSIRE WEIGHTS\t\t\t DAM WEIGHTS") +print("-"*20 + "\t\t " + "-"*20) +print(" i w_std w_rbs\t\t i w_std w_rbs") +for candidate in range(25): + print(f"{candidate*2:02d} {w_std.X[candidate*2]:.5f} {w_rbs.X[candidate*2]:.5f} \ + {candidate*2+1:02d} {w_std.X[candidate*2+1]:.5f} {w_rbs.X[candidate*2+1]:.5f}") \ No newline at end of file diff --git a/Scratch/numerical-error.ipynb b/Scratch/numerical-error.ipynb index 73a59c2..63583a3 100644 --- a/Scratch/numerical-error.ipynb +++ b/Scratch/numerical-error.ipynb @@ -282,7 +282,7 @@ "\n", "# setup the robust objective function\n", "model.setObjective(\n", - " 0.5*w@(relationship_matrix@w) - lam*w.transpose()@expected_breeding_values - kappa*z,\n", + " 0.5*w@(relationship_matrix@w) - lam*w.transpose()@expected_breeding_values - lam*kappa*z,\n", "GRB.MINIMIZE)\n", "\n", "# add quadratic uncertainty constraint\n", diff --git a/Scratch/robust-notes.ipynb b/Scratch/robust-notes.ipynb index 31e5b60..5f0fec7 100644 --- a/Scratch/robust-notes.ipynb +++ b/Scratch/robust-notes.ipynb @@ -11,7 +11,7 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 1, "metadata": {}, "outputs": [], "source": [ @@ -32,7 +32,7 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ @@ -48,7 +48,7 @@ }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 3, "metadata": {}, "outputs": [], "source": [ @@ -71,7 +71,7 @@ }, { "cell_type": "code", - "execution_count": 18, + "execution_count": 4, "metadata": {}, "outputs": [], "source": [ @@ -107,7 +107,7 @@ }, { "cell_type": "code", - "execution_count": 19, + "execution_count": 5, "metadata": {}, "outputs": [], "source": [ @@ -146,7 +146,7 @@ }, { "cell_type": "code", - "execution_count": 20, + "execution_count": 6, "metadata": {}, "outputs": [], "source": [ @@ -210,7 +210,7 @@ }, { "cell_type": "code", - "execution_count": 21, + "execution_count": 7, "metadata": {}, "outputs": [], "source": [ @@ -240,7 +240,7 @@ }, { "cell_type": "code", - "execution_count": 22, + "execution_count": 8, "metadata": {}, "outputs": [], "source": [ @@ -308,7 +308,7 @@ }, { "cell_type": "code", - "execution_count": 23, + "execution_count": 9, "metadata": {}, "outputs": [ { @@ -333,7 +333,7 @@ " [ 0. 0. -1. -1. 2.615 0.615 -1.231]\n", " [-1. 0. 0. -1. 0.615 2.615 -1.231]\n", " [ 0. 0. 0. 0. -1.231 -1.231 2.462]]\n", - "Numpy inversion in 0.00021 seconds\n", + "Numpy inversion in 0.00093 seconds\n", "\n", "A^-1 =\t(from pedigree)\n", " [[ 2.333 0.5 -0.667 -0.5 0. -1. 0. ]\n", @@ -343,7 +343,7 @@ " [ 0. 0. -1. -1. 2.615 0.615 -1.231]\n", " [-1. 0. 0. -1. 0.615 2.615 -1.231]\n", " [ 0. 0. 0. 0. -1.231 -1.231 2.462]]\n", - "Henderson inversion in 0.00034 seconds\n", + "Henderson inversion in 0.00033 seconds\n", "\n" ] } @@ -441,7 +441,7 @@ }, { "cell_type": "code", - "execution_count": 24, + "execution_count": 10, "metadata": {}, "outputs": [], "source": [ @@ -472,7 +472,7 @@ }, { "cell_type": "code", - "execution_count": 25, + "execution_count": 11, "metadata": {}, "outputs": [], "source": [ @@ -497,13 +497,15 @@ }, { "cell_type": "code", - "execution_count": 26, + "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ + "Set parameter Username\n", + "Academic license - for non-commercial use only - expires 2025-02-26\n", "QP solution: w = [0.5 0.375 0.125]\n" ] } @@ -544,7 +546,7 @@ }, { "cell_type": "code", - "execution_count": 27, + "execution_count": 13, "metadata": {}, "outputs": [ { @@ -589,7 +591,7 @@ " 7 -8.12500000e-01 -8.12500040e-01 0.00e+00 2.78e-17 1.01e-08 0s\n", " 8 -8.12500000e-01 -8.12500000e-01 2.22e-16 5.55e-17 1.01e-11 0s\n", "\n", - "Barrier solved model in 8 iterations and 0.03 seconds (0.00 work units)\n", + "Barrier solved model in 8 iterations and 0.02 seconds (0.00 work units)\n", "Optimal objective -8.12500000e-01\n", "\n", "[0.5 0.375 0.125]\n", @@ -598,11 +600,11 @@ } ], "source": [ - "# define robust optimization variables\n", - "omega = np.diagflat(np.random.rand(problem_size, 1)) # TODO work out how best to set this param\n", - "kappa = float(np.random.rand(1)) # TODO work out how best to set this param\n", + "# randomly set until received from Gregor \n", + "omega = np.diagflat(np.random.rand(problem_size, 1))\n", + "kappa = float(np.random.rand(1))\n", "\n", - "# TODO a temp(?) workaround to get the square root with SciPy rather than Numpy.\n", + "# a temp(?) workaround to get the square root with SciPy rather than Numpy.\n", "from scipy.linalg import sqrtm\n", "root_omega = sqrtm(omega)\n", "\n", @@ -618,15 +620,19 @@ " # STANDARD OPTIMIZATION OBJECTIVE\n", " 0.5*w@(relationship_matrix@w) - lam*w.transpose()@expected_breeding_values, # standard optimization\n", "\n", - " # ROBUST OPTIMIZATION OBJECTIVE WITH QUADRATIC UNCERTAINTY \n", + " # ROBUST OPTIMIZATION OBJECTIVE WITH QUADRATIC UNCERTAINTY\n", + " # cannot use np.power, get:\n", + " # TypeError: unsupported operand type(s) for ** or pow(): 'MQuadExpr' and 'float'\n", " # 0.5*w@(relationship_matrix@w) - lam*w.transpose()@expected_breeding_values - kappa*np.power(np.inner(w, omega@w), 0.5), # inner product representation\n", - " # BUG TypeError: unsupported operand type(s) for ** or pow(): 'MQuadExpr' and 'float'\n", + " # \n", + " # cannot use scipy.linalg.sqrtm and np.nla.norm(ord=2), get:\n", + " # ValueError: Improper number of dimensions to norm.\n", " # 0.5*w@(relationship_matrix@w) - lam*w.transpose()@expected_breeding_values - kappa*nla.norm(root_omega@w, ord=2), # norm representation\n", - " # BUG ValueError: Improper number of dimensions to norm.\n", "\n", " # ROBUST OPTIMIZATION OBJECTIVE WITH BOX UNCERTAINTY\n", + " # cannot use np.avs, get:\n", + " # TypeError: bad operand type for abs(): 'MQuadExpr'\n", " # 0.5*w@(relationship_matrix@w) - lam*w.transpose()@expected_breeding_values - kappa*np.abs(np.inner(w, omega@w)), # robust with box uncertainty\n", - " # BUG TypeError: bad operand type for abs(): 'MQuadExpr'\n", " GRB.MINIMIZE)\n", "\n", " # add sub-to-half constraints\n", @@ -650,7 +656,7 @@ }, { "cell_type": "code", - "execution_count": 28, + "execution_count": 14, "metadata": {}, "outputs": [ { @@ -663,13 +669,13 @@ "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n", "\n", "Optimize a model with 23 rows, 4 columns and 25 nonzeros\n", - "Model fingerprint: 0xad441817\n", + "Model fingerprint: 0xe8d7361e\n", "Model has 3 quadratic objective terms\n", "Model has 3 quadratic constraints\n", "Coefficient statistics:\n", " Matrix range [1e+00, 1e+00]\n", - " QMatrix range [3e-01, 1e+00]\n", - " Objective range [3e-01, 2e+00]\n", + " QMatrix range [4e-01, 1e+00]\n", + " Objective range [5e-01, 2e+00]\n", " QObjective range [1e+00, 5e+00]\n", " Bounds range [0e+00, 0e+00]\n", " RHS range [5e-01, 1e+00]\n", @@ -687,31 +693,31 @@ "\n", " Objective Residual\n", "Iter Primal Dual Primal Dual Compl Time\n", - " 0 -8.91467504e-01 -8.91467504e-01 1.52e+00 7.59e-01 7.11e-01 0s\n", - " 1 -9.19418849e-01 -1.57169532e+00 1.15e-01 8.34e-07 9.04e-02 0s\n", - " 2 -8.16296575e-01 -9.27971330e-01 1.27e-07 1.17e-08 9.31e-03 0s\n", - " 3 -8.33236487e-01 -8.39356270e-01 1.43e-13 1.28e-14 5.10e-04 0s\n", - " 4 -8.36888293e-01 -8.36915991e-01 4.33e-14 3.33e-16 2.31e-06 0s\n", - " 5 -8.36893763e-01 -8.36893897e-01 9.70e-13 3.94e-14 1.12e-08 0s\n", + " 0 -9.02965386e-01 -9.02965386e-01 1.51e+00 7.38e-01 6.76e-01 0s\n", + " 1 -9.33513427e-01 -1.53929722e+00 1.29e-01 8.12e-07 8.98e-02 0s\n", + " 2 -8.59473155e-01 -1.00199951e+00 1.42e-07 1.99e-08 1.19e-02 0s\n", + " 3 -8.66933672e-01 -8.68757572e-01 7.54e-10 1.48e-10 1.52e-04 0s\n", + " 4 -8.67494911e-01 -8.67498866e-01 9.07e-13 1.31e-13 3.30e-07 0s\n", + " 5 -8.67495578e-01 -8.67495669e-01 1.93e-10 2.16e-13 7.63e-09 0s\n", "\n", "Barrier solved model in 5 iterations and 0.02 seconds (0.00 work units)\n", - "Optimal objective -8.36893763e-01\n", + "Optimal objective -8.67495578e-01\n", "\n", "\n", " Variable X \n", "-------------------------\n", " w[0] 0.5 \n", - " w[1] 0.352596 \n", - " w[2] 0.147404 \n", - " z 0.0813294 \n", - "Obj: -0.836894\n" + " w[1] 0.328648 \n", + " w[2] 0.171352 \n", + " z 0.133483 \n", + "Obj: -0.867496\n" ] } ], "source": [ - "# define robust optimization variables\n", - "omega = np.diagflat(np.random.rand(problem_size, 1)) # TODO work out how best to set this param\n", - "kappa = float(np.random.rand(1)) # TODO work out how best to set this param\n", + "# randomly set until received from Gregor\n", + "omega = np.diagflat(np.random.rand(problem_size, 1))\n", + "kappa = float(np.random.rand(1))\n", "\n", "try:\n", " # create a new model for robust genetic selection\n", diff --git a/robust-genetics.ipynb b/robust-genetics.ipynb index 064fdb5..178193f 100644 --- a/robust-genetics.ipynb +++ b/robust-genetics.ipynb @@ -6,9 +6,7 @@ "source": [ "# Robust Optimization for Genetic Selection\n", "\n", - "This is a working notebook, looking at how the quadratic optimization problem (QP) which arises in the context of robust genetic selection can be solved with Python (tested with 3.10 specifically). There are some standard packages which this depends on, imported below.\n", - "\n", - "" + "This notebook provides a didactic for using Python to solve the multi-objective optimization problem which arises in the context of robust genetic selection. It has been tested with Python 3.10 specifically and there are some standard packages which this depends on, imported below." ] }, { @@ -240,19 +238,65 @@ "\n", "A wrinkle to this is that we don't typically our problem variables with certainty. There _are_ ways to define $\\Sigma$ based on relationships within the cohort (which are known and discussed more later) but $\\mu$ is an estimated variable. In particular we say $\\mu$ has a univariate normal distribution, $\\mu\\sim N(\\bar{\\mu}, \\Omega)$. This means we must turn to optimization tools which can address this uncertainty.¹\n", "\n", - "Robust optimization is one such tool in which we adjust the objective function to model the inherent uncertainty in the problem. We may either do this with a quadratic uncertainty set, in which case our objective has an additional square-root term as in\n", + "Robust optimization is one such tool in which we adjust the objective function to model the inherent uncertainty in the problem. Within this framework, the problem would be the bilevel optimization problem\n", "$$\n", - " \\min_w \\frac{1}{2}w^{T}\\Sigma w - \\lambda w^{T}\\mu - \\kappa\\sqrt{w^{T}\\Omega w}\\ \\text{ subject to }\\ Mw = \\begin{bmatrix} 0.5 \\\\ 0.5\\end{bmatrix},\\ l\\leq w\\leq u,\n", + " \\min_w \\frac{1}{2}w^{T}\\Sigma w - \\lambda\\left( \\max_{\\mu\\in U_{\\mu}} w^{T}\\mu \\right)\\ \\text{ subject to }\\ Mw = \\begin{bmatrix} 0.5 \\\\ 0.5\\end{bmatrix},\\ l\\leq w\\leq u.\n", "$$\n", - "or with a box uncertainty set, in which case our objective has an additional absolute value term as in\n", + "The result of the inner problem determines the outer problem in a follower-leader setup. However, the inner problem is relatively simple and can be solved explicitly based on our choice of $U_{\\mu}$, representing an uncertainty set around the value of $\\mu$.\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Quadratic Uncertainty Sets\n", + "\n", + "There are many ways to define the uncertainty set $U_{\\mu}$, but for practical reasons relating to continuity and differentiability² it's beneficial to define it using a quadratic uncertainty set. We say that\n", "$$\n", - " \\min_w \\frac{1}{2}w^{T}\\Sigma w - \\lambda w^{T}\\mu - \\kappa\\|\\Omega^{\\frac{1}{2}} w\\|\\ \\text{ subject to }\\ Mw = \\begin{bmatrix} 0.5 \\\\ 0.5\\end{bmatrix},\\ l\\leq w\\leq u,\n", + " U_{\\mu} := \\left\\lbrace \\mu :\\ {(\\mu-\\bar{\\mu})}^{T}\\Omega^{-1}(\\mu-\\bar{\\mu}) \\leq \\kappa^2 \\right\\rbrace\n", "$$\n", - "where $\\kappa\\in\\mathbb{R}$ is our robust optimization parameters. For practical reasons relating to continuity and differentiability, the quadratic uncertainty set is far more favourable to work with.\n", - "\n", - "This is obviously no longer a quadratic problem, so `qpsolvers` is no longer a viable tool. We will instead now need to work with the Gurobi API directly.\n", - "\n", - "" + "which in effect bounds the uncertainty in a ball about $\\bar{\\mu}$ with parameter $\\kappa$.³ This means that our lower level problem $\\max_{\\mu\\in U_{\\mu}} w^{T}\\mu$ becomes\n", + "$$\n", + " \\max_{\\mu} w^{T}\\mu\\quad \\text{subject to}\\quad {(\\mu-\\bar{\\mu})}^{T}\\Omega^{-1}(\\mu-\\bar{\\mu}) \\leq \\kappa^2,\n", + "$$\n", + "or, in standard form,\n", + "$$\n", + " \\min_{\\mu} \\left(-w^{T}\\mu\\right) \\quad\\text{subject to}\\quad \\kappa^2 - {(\\mu-\\bar{\\mu})}^{T}\\Omega^{-1}(\\mu-\\bar{\\mu}) \\geq 0.\n", + "$$\n", + "This is convex, since $\\Omega$ is a positive definite matrix and $-w^{T}\\mu$ is an affine function. If we define a Lagrangian multiplier $\\rho\\in\\mathbb{R}$, the KKT conditions of this problem are:\n", + "\\begin{align}\n", + " \\nabla_{\\mu}L(\\mu, \\rho) = 0 \\quad&\\Rightarrow\\quad \\nabla_{\\mu}\\left( -w^{T}\\mu - \\rho\\big( \\kappa^2 - {(\\mu-\\bar{\\mu})}^{T}\\Omega^{-1}(\\mu-\\bar{\\mu}) \\big) \\right) = 0, \\\\\n", + " c(\\mu) \\geq 0 \\quad&\\Rightarrow\\quad \\kappa^2 - {(\\mu-\\bar{\\mu})}^{T}\\Omega^{-1}(\\mu-\\bar{\\mu}) \\geq 0, \\\\\n", + " \\rho \\geq 0 \\quad&\\Rightarrow\\quad \\rho\\geq0, \\\\\n", + " \\rho c(\\mu) = 0 \\quad&\\Rightarrow\\quad \\rho\\left(\\kappa^2 - {(\\mu-\\bar{\\mu})}^{T}\\Omega^{-1}(\\mu-\\bar{\\mu}) \\right) = 0.\n", + "\\end{align}\n", + "From (1) we have that\n", + "$$\n", + " -w + \\rho 2\\Omega^{-1}(\\mu-\\bar{\\mu}) = 0 \\quad\\Rightarrow\\quad \\mu - \\bar{\\mu} = \\frac{1}{2\\rho}\\Omega w \\qquad (5)\n", + "$$\n", + "\\end{align*}\n", + "which when substituted into (4) gives\n", + "\\begin{align*}\n", + " &\\phantom{\\Rightarrow}\\quad \\rho\\left( \\kappa^2 - {(\\mu-\\bar{\\mu})}^{T}\\Omega^{-1}(\\mu-\\bar{\\mu}) \\right) = 0 \\\\\n", + " &\\Rightarrow\\quad \\rho\\left( \\kappa^2 - \\big(\\frac{1}{2\\rho}\\Omega w\\big)^{T}\\Omega^{-1}\\big(\\frac{1}{2\\rho}\\Omega w\\big) \\right) = 0 \\\\\n", + " &\\Rightarrow\\quad \\rho\\kappa^2 - \\frac{1}{4\\rho}w^{T}\\Omega^{T}\\Omega^{-1}\\Omega w = 0 \\\\\n", + " &\\Rightarrow\\quad \\rho^2\\kappa^2 = \\frac{1}{4}w^{T}\\Omega w\\quad \\text{(since $\\Omega$ symmetric)} \\\\\n", + " &\\Rightarrow\\quad \\rho\\kappa = \\frac{1}{2}\\sqrt{w^{T}\\Omega w}\\quad \\text{(since $\\rho,\\kappa\\geq0$)} \\\\\n", + " &\\Rightarrow\\quad \\frac{1}{2\\rho} = \\frac{\\kappa}{\\sqrt{w^{T}\\Omega w}}.\n", + "\\end{align*}\n", + "Substituting this back into (5), we find that\n", + "$$\n", + " \\mu - \\bar{\\mu} = \\frac{\\kappa}{\\sqrt{w^{T}\\Omega w}}\\Omega w \\quad\\Rightarrow\\quad \\mu = \\bar{\\mu} + \\frac{\\kappa}{\\sqrt{w^{T}\\Omega w}}\\Omega w\n", + "$$\n", + "as the solution for the inner problem. Substituting this back into the outer problem gives\n", + "$$\n", + " \\min_w \\frac{1}{2}w^{T}\\Sigma w - \\lambda w^{T}\\left( \\bar{\\mu} + \\frac{\\kappa}{\\sqrt{w^{T}\\Omega w}}\\Omega w \\right)\\ \\text{ subject to }\\ Mw = \\begin{bmatrix} 0.5 \\\\ 0.5\\end{bmatrix},\\ l\\leq w\\leq u.\n", + "$$\n", + "which after simplifying down and rationalizing the denominator becomes\n", + "$$\n", + " \\min_w \\frac{1}{2}w^{T}\\Sigma w - \\lambda w^{T}\\bar{\\mu} - \\lambda\\kappa\\sqrt{w^{T}\\Omega w} \\ \\text{ subject to }\\ Mw = \\begin{bmatrix} 0.5 \\\\ 0.5\\end{bmatrix},\\ l\\leq w\\leq u\n", + "$$\n", + "where $\\kappa\\in\\mathbb{R}$ is our robust optimization parameter. Since our objective has gained an additional square root term, this is obviously no longer a quadratic problem and `qpsolvers` is no longer a viable tool. We will instead now need to work with the Gurobi API directly." ] }, { @@ -347,16 +391,17 @@ "source": [ "Unfortunately Gurobi cannot handle our problem with its objective function in the form\n", "$$\n", - " \\min_w \\frac{1}{2}w^{T}\\Sigma w - \\lambda w^{T}\\mu - \\kappa\\sqrt{w^{T}\\Omega w}\\ \\text{ subject to }\\ Mw = \\begin{bmatrix} 0.5 \\\\ 0.5\\end{bmatrix},\\ l\\leq w\\leq u,\n", + " \\min_w \\frac{1}{2}w^{T}\\Sigma w - \\lambda w^{T}\\mu - \\lambda\\kappa\\sqrt{w^{T}\\Omega w}\\ \\text{ subject to }\\ Mw = \\begin{bmatrix} 0.5 \\\\ 0.5\\end{bmatrix},\\ l\\leq w\\leq u,\n", "$$\n", "so some further adjustments are needed first. If we define a real auxillary variable $z\\geq0$ such that $z\\leq\\sqrt{w^{T}\\Omega w}$, then our problem becomes\n", "$$\n", - " \\min_w \\frac{1}{2}w^{T}\\Sigma w - \\lambda w^{T}\\mu - \\kappa z\\ \\text{ s.t. }\\ z\\leq\\sqrt{w^{T}\\Omega w},\\ Mw = \\begin{bmatrix} 0.5 \\\\ 0.5\\end{bmatrix},\\ l\\leq w\\leq u.\n", + " \\min_w \\frac{1}{2}w^{T}\\Sigma w - \\lambda w^{T}\\mu - \\lambda\\kappa z\\ \\text{ s.t. }\\ z\\leq\\sqrt{w^{T}\\Omega w},\\ Mw = \\begin{bmatrix} 0.5 \\\\ 0.5\\end{bmatrix},\\ l\\leq w\\leq u.\n", "$$\n", - "\n", + "TODO: _better explain the idea of $z$ pushing up so that this switch doesn't have any practical difference._\n", + "\n", "However, Gurobi _still_ can't handle this due to the presence of the square root, so we further make use of both $z$ and $\\sqrt{w^{T}\\Omega w}$ being positive to note that $z\\leq\\sqrt{w^{T}\\Omega w}$ can be squared on both sides:\n", "$$\n", - " \\min_w \\frac{1}{2}w^{T}\\Sigma w - \\lambda w^{T}\\mu - \\kappa z\\ \\text{ s.t. }\\ z^2\\leq w^{T}\\Omega w,\\ Mw = \\begin{bmatrix} 0.5 \\\\ 0.5\\end{bmatrix},\\ l\\leq w\\leq u.\n", + " \\min_w \\frac{1}{2}w^{T}\\Sigma w - \\lambda w^{T}\\mu - \\lambda\\kappa z\\ \\text{ s.t. }\\ z^2\\leq w^{T}\\Omega w,\\ Mw = \\begin{bmatrix} 0.5 \\\\ 0.5\\end{bmatrix},\\ l\\leq w\\leq u.\n", "$$" ] }, @@ -364,12 +409,13 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We will define\n", + "To explore this with our toy problem from before, we will define\n", "$$\n", " \\bar{\\mu} = \\begin{bmatrix} 1 \\\\ 5 \\\\ 2 \\end{bmatrix},\\quad\n", - " \\Omega = \\begin{bmatrix} 1 & 0 & 0 \\\\ 0 & 4 & 0 \\\\ 0 & 0 & \\frac{1}{8} \\end{bmatrix}, \\kappa = 0.5\n", + " \\Omega = \\begin{bmatrix} 1 & 0 & 0 \\\\ 0 & 4 & 0 \\\\ 0 & 0 & \\frac{1}{8} \\end{bmatrix},\\quad\n", + " \\kappa = 0.5\n", "$$\n", - "and our other problem variables as before." + "and retain the other variables as before." ] }, { @@ -409,13 +455,13 @@ "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n", "\n", "Optimize a model with 23 rows, 4 columns and 25 nonzeros\n", - "Model fingerprint: 0x541bef89\n", + "Model fingerprint: 0xb99edb8a\n", "Model has 3 quadratic objective terms\n", "Model has 3 quadratic constraints\n", "Coefficient statistics:\n", " Matrix range [1e+00, 1e+00]\n", " QMatrix range [1e-01, 4e+00]\n", - " Objective range [5e-01, 2e+00]\n", + " Objective range [2e-01, 2e+00]\n", " QObjective range [1e+00, 5e+00]\n", " Bounds range [0e+00, 0e+00]\n", " RHS range [5e-01, 1e+00]\n", @@ -433,19 +479,19 @@ "\n", " Objective Residual\n", "Iter Primal Dual Primal Dual Compl Time\n", - " 0 -8.68156934e-01 -8.68156934e-01 1.57e+00 7.38e-01 6.17e-01 0s\n", - " 1 -9.99130241e-01 -1.52016137e+00 1.66e-01 8.12e-07 8.63e-02 0s\n", - " 2 -8.04900486e-01 -9.40830018e-01 1.82e-07 1.07e-08 1.13e-02 0s\n", - " 3 -8.32538184e-01 -8.51116759e-01 3.42e-13 1.18e-14 1.55e-03 0s\n", - " 4 -8.36522462e-01 -8.36788831e-01 1.55e-14 4.44e-16 2.22e-05 0s\n", - " 5 -8.36548538e-01 -8.36556228e-01 6.67e-13 4.47e-14 6.41e-07 0s\n", - " 6 -8.36549617e-01 -8.36550259e-01 1.09e-11 5.45e-14 5.35e-08 0s\n", + " 0 -8.36906934e-01 -8.36906934e-01 1.57e+00 7.34e-01 6.19e-01 0s\n", + " 1 -9.44994564e-01 -1.48469043e+00 1.36e-01 8.08e-07 8.01e-02 0s\n", + " 2 -7.95781338e-01 -8.91732825e-01 1.50e-07 1.21e-08 8.00e-03 0s\n", + " 3 -8.15171524e-01 -8.29882203e-01 1.65e-13 1.33e-14 1.23e-03 0s\n", + " 4 -8.24009582e-01 -8.24184428e-01 3.18e-14 3.20e-16 1.46e-05 0s\n", + " 5 -8.24036160e-01 -8.24039018e-01 1.54e-13 7.22e-15 2.38e-07 0s\n", + " 6 -8.24036608e-01 -8.24036837e-01 1.10e-10 1.94e-13 1.90e-08 0s\n", "\n", - "Barrier solved model in 6 iterations and 0.02 seconds (0.00 work units)\n", - "Optimal objective -8.36549617e-01\n", + "Barrier solved model in 6 iterations and 0.03 seconds (0.00 work units)\n", + "Optimal objective -8.24036608e-01\n", "\n", - "w = [0.50000 0.35282 0.14718],\n", - "z = 0.052036865323801265.\n" + "w = [0.50000 0.36395 0.13605],\n", + "z = 0.04810269528692041.\n" ] } ], @@ -459,7 +505,7 @@ "\n", "# setup the robust objective function\n", "model.setObjective(\n", - " 0.5*w@(relationship_matrix@w) - lam*w.transpose()@expected_breeding_values - kappa*z,\n", + " 0.5*w@(relationship_matrix@w) - lam*w.transpose()@expected_breeding_values - lam*kappa*z,\n", "GRB.MINIMIZE)\n", "\n", "# add quadratic uncertainty constraint\n", @@ -485,17 +531,28 @@ "| $\\kappa$ | $w$ | $z$ | $f(w,z)$ |\n", "| ---: | ------------------------: | ------: | -------: |\n", "| 0.0 | [0.50000 0.37500 0.12500] | 0.02756 | -0.81250 |\n", - "| 0.5 | [0.50000 0.35282 0.14718] | 0.05204 | -0.83655 |\n", - "| 1.0 | [0.50000 0.33073 0.16927] | 0.05985 | -0.86451 |\n", - "| 2.0 | [0.50000 0.28663 0.21337] | 0.07544 | -0.93214 |\n", - "| 4.0 | [0.50000 0.19816 0.30184] | 0.10671 | -1.11428 |\n", - "| 8.0 | [0.50000 0.07511 0.42489] | 0.15022 | -1.65453 |\n", - "| 16.0 | [0.50000 0.07511 0.42489] | 0.15022 | -2.85630 |\n", + "| 0.5 | [0.50000 0.36395 0.13605] | 0.04810 | -0.82404 |\n", + "| 1.0 | [0.50000 0.35282 0.14718] | 0.05204 | -0.83655 |\n", + "| 2.0 | [0.50000 0.33073 0.16927] | 0.05985 | -0.86451 |\n", + "| 4.0 | [0.50000 0.28663 0.21337] | 0.07544 | -0.93214 |\n", + "| 8.0 | [0.50000 0.19816 0.30184] | 0.10671 | -1.11428 |\n", + "| 16.0 | [0.50000 0.07511 0.42489] | 0.15022 | -1.65453 |\n", + "| 32.0 | [0.50000 0.07511 0.42489] | 0.15022 | -2.85630 |\n", "\n", "We can see that in the case $\\kappa = 0$ we return the standard optimization solution, as expected. However, we do have $z\\neq0$ which suggests that a small amount of numerical error has been introduced as a result of asking Gurobi to optimize $z$ as well as $w$.\n", - "\n", "\n", - "\n" + "We also note that there are not further changes once $\\kappa$ falls above a certain threshold." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Verifying with KKTs\n", + "\n", + "As with other(?) optimization problems, we can verify that we have produced an optimum solution by checking that it satisfies the KKT conditions.\n", + "\n", + "TODO: _check that these values are correct, talk about subbing into the KKT conditions._" ] }, { @@ -524,7 +581,13 @@ "source": [ "def load_problem(A_filename, E_filename, S_filename, dimension=False):\n", " \"\"\"\n", - " TODO write a docstring\n", + " Used to load genetic selection problems into NumPy. It takes three\n", + " string inputs for filenames where Sigma, Mu, and Omega are stored,\n", + " as well as an optional integer input for problem dimension if this\n", + " is known. If it's know know, it's worked out based on E_filename.\n", + "\n", + " As output, it returns (A, E, S, n), where A and S are n-by-n NumPy\n", + " arrays, E is a length n NumPy array, and n is an integer.\n", " \"\"\"\n", "\n", " def load_symmetric_matrix(filename, dimension):\n", @@ -621,7 +684,7 @@ " 10 -9.76215438e-01 -9.76215482e-01 7.03e-13 6.66e-16 8.86e-10 0s\n", " 11 -9.76215473e-01 -9.76215476e-01 4.24e-13 6.66e-16 4.91e-11 0s\n", "\n", - "Barrier solved model in 11 iterations and 0.03 seconds (0.00 work units)\n", + "Barrier solved model in 11 iterations and 0.04 seconds (0.00 work units)\n", "Optimal objective -9.76215473e-01\n", "\n", "Gurobi Optimizer version 11.0.0 build v11.0.0rc2 (linux64 - \"Ubuntu 22.04.4 LTS\")\n", @@ -666,39 +729,42 @@ "Cutting planes:\n", " RLT: 6\n", "\n", - "Explored 1 nodes (1828 simplex iterations) in 0.34 seconds (0.10 work units)\n", + "Explored 1 nodes (1828 simplex iterations) in 0.29 seconds (0.10 work units)\n", "Thread count was 8 (of 8 available processors)\n", "\n", "Solution count 1: -0.976358 \n", "\n", "Optimal solution found (tolerance 5.00e-02)\n", "Best objective -9.763583540500e-01, best bound -1.025012483088e+00, gap 4.9832%\n", - " i w_std w_rbs \t\t i w_std w_rbs\n", - "00 0.00000 0.00002\t 01 0.08319 0.08339\n", - "02 0.11649 0.11508\t 03 0.00000 0.00001\n", - "04 0.00000 0.00002\t 05 0.00000 0.00002\n", - "06 0.00000 0.00001\t 07 0.06817 0.06830\n", - "08 0.00000 0.00001\t 09 0.03762 0.03701\n", - "10 0.00000 0.00002\t 11 0.00000 0.00002\n", - "12 0.00000 0.00001\t 13 0.00355 0.00564\n", - "14 0.11542 0.11450\t 15 0.06966 0.06913\n", - "16 0.00000 0.00006\t 17 0.00000 0.00001\n", - "18 0.00000 0.00002\t 19 0.00000 0.00003\n", - "20 0.00000 0.00002\t 21 0.00001 0.00062\n", - "22 0.00000 0.00017\t 23 0.00000 0.00004\n", - "24 0.03537 0.03493\t 25 0.00000 0.00001\n", - "26 0.00000 0.00001\t 27 0.00000 0.00002\n", - "28 0.00000 0.00002\t 29 0.00000 0.00002\n", - "30 0.00000 0.00001\t 31 0.03580 0.03612\n", - "32 0.08447 0.08517\t 33 0.00000 0.00003\n", - "34 0.00000 0.00001\t 35 0.02849 0.02747\n", - "36 0.13489 0.13487\t 37 0.07446 0.07264\n", - "38 0.00000 0.00001\t 39 0.00000 0.00001\n", - "40 0.00000 0.00003\t 41 0.07312 0.07388\n", - "42 0.00000 0.00004\t 43 0.00000 0.00003\n", - "44 0.01336 0.01485\t 45 0.00000 0.00002\n", - "46 0.00000 0.00009\t 47 0.02593 0.02555\n", - "48 0.00000 0.00003\t 49 0.00000 0.00000\n" + "\n", + "SIRE WEIGHTS\t\t\t DAM WEIGHTS\n", + "--------------------\t\t --------------------\n", + " i w_std w_rbs\t\t i w_std w_rbs\n", + "00 0.00000 0.00002 01 0.08319 0.08339\n", + "02 0.11649 0.11508 03 0.00000 0.00001\n", + "04 0.00000 0.00002 05 0.00000 0.00002\n", + "06 0.00000 0.00001 07 0.06817 0.06830\n", + "08 0.00000 0.00001 09 0.03762 0.03701\n", + "10 0.00000 0.00002 11 0.00000 0.00002\n", + "12 0.00000 0.00001 13 0.00355 0.00564\n", + "14 0.11542 0.11450 15 0.06966 0.06913\n", + "16 0.00000 0.00006 17 0.00000 0.00001\n", + "18 0.00000 0.00002 19 0.00000 0.00003\n", + "20 0.00000 0.00002 21 0.00001 0.00062\n", + "22 0.00000 0.00017 23 0.00000 0.00004\n", + "24 0.03537 0.03493 25 0.00000 0.00001\n", + "26 0.00000 0.00001 27 0.00000 0.00002\n", + "28 0.00000 0.00002 29 0.00000 0.00002\n", + "30 0.00000 0.00001 31 0.03580 0.03612\n", + "32 0.08447 0.08517 33 0.00000 0.00003\n", + "34 0.00000 0.00001 35 0.02849 0.02747\n", + "36 0.13489 0.13487 37 0.07446 0.07264\n", + "38 0.00000 0.00001 39 0.00000 0.00001\n", + "40 0.00000 0.00003 41 0.07312 0.07388\n", + "42 0.00000 0.00004 43 0.00000 0.00003\n", + "44 0.01336 0.01485 45 0.00000 0.00002\n", + "46 0.00000 0.00009 47 0.02593 0.02555\n", + "48 0.00000 0.00003 49 0.00000 0.00000\n" ] } ], @@ -710,7 +776,7 @@ " 50)\n", "\n", "lam = 0.5\n", - "kappa = 1\n", + "kappa = 2\n", "\n", "# define the M so that column i is [1;0] if i is a sire (so even) and [0;1] otherwise \n", "M = np.zeros((2, n))\n", @@ -734,8 +800,11 @@ "GRB.MINIMIZE)\n", "\n", "model_rbs.setObjective(\n", - " 0.5*w_rbs@(sigma@w_rbs) - lam*w_rbs.transpose()@mubar - kappa*z_rbs,\n", - "GRB.MINIMIZE) # TODO check if there's a way to define this objective in terms of the other\n", + " # Gurobi does offer a way to set one objective in terms of another, i.e.\n", + " # we could use `model_std.getObjective() - lam*kappa*z_rbs` to define this\n", + " # robust objective, but it results in a significant slowdown in code.\n", + " 0.5*w_rbs@(sigma@w_rbs) - lam*w_rbs.transpose()@mubar - lam*kappa*z_rbs,\n", + "GRB.MINIMIZE)\n", "\n", "# add sum-to-half constraints to both models\n", "model_std.addConstr(M @ w_std == m, name=\"sum-to-half\")\n", @@ -760,17 +829,19 @@ "model_rbs.optimize()\n", "\n", "# HACK code which prints the results for comparison in a nice format\n", - "print(\" i w_std w_rbs \\t\\t i w_std w_rbs\")\n", - "for candidate in range(0,25):\n", - " print(f\"{candidate*2:02d} {w_std.X[candidate*2]:.5f} {w_rbs.X[candidate*2]:.5f}\\t \\\n", - " {candidate*2+1:02d} {w_std.X[candidate*2+1]:.5f} {w_rbs.X[candidate*2+1]:.5f}\")" + "print(\"\\nSIRE WEIGHTS\\t\\t\\t DAM WEIGHTS\")\n", + "print(\"-\"*20 + \"\\t\\t \" + \"-\"*20)\n", + "print(\" i w_std w_rbs\\t\\t i w_std w_rbs\")\n", + "for candidate in range(25):\n", + " print(f\"{candidate*2:02d} {w_std.X[candidate*2]:.5f} {w_rbs.X[candidate*2]:.5f} \\\n", + " {candidate*2+1:02d} {w_std.X[candidate*2+1]:.5f} {w_rbs.X[candidate*2+1]:.5f}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "We can see that with parameters $\\lambda = 0.5, \\kappa = 1$, we go from having most candidates having zero contribution in the standard solution to most candidates having a contribution of order 0.001 in the robust solution." + "We can see that with parameters $\\lambda = 0.5, \\kappa = 2$, we go from having most candidates having zero contribution in the standard solution to most candidates having a contribution of order 0.00001 in the robust solution." ] }, { @@ -779,7 +850,30 @@ "source": [ "## Footnotes\n", "\n", - "1. In reality, $\\Omega$ is also an estimated variable but we can ignore this for now to avoid going down a rabbit hole of uncertainties." + "1. In reality, $\\Omega$ is also an estimated variable but we can ignore this for now to avoid going down a rabbit hole of uncertainties.\n", + "\n", + "2. For example, we could alternatively define $U_{\\mu}$ using a box uncertainty set\n", + " $$\n", + " U_{\\mu} := \\left\\lbrace \\mu :\\ | \\mu_i - \\bar{\\mu}_i| \\leq\\xi_i,\\ \\forall i \\right\\rbrace\n", + " $$\n", + " where $\\xi$ is some upper bound on the uncertainty value. In this case, work by [Heckel et al (2016)](https://www.risk.net/journal-of-investment-strategies/2475975/insights-into-robust-optimization-decomposing-into-mean-variance-and-risk-based-portfolios) shows that our objective function would gain an additional absolute value term, becoming\n", + " $$\n", + " \\min_w \\frac{1}{2}w^{T}\\Sigma w - \\lambda w^{T}\\mu - \\kappa\\|\\Omega^{\\frac{1}{2}} w\\|\\ \\text{ subject to }\\ Mw = \\begin{bmatrix} 0.5 \\\\ 0.5\\end{bmatrix},\\ l\\leq w\\leq u.\n", + " $$\n", + " This absolute value term means we lose many of the properties that we would like to work with in our objective functions.\n", + "\n", + "3. By \"ball\" we do not just mean a sphere in the typical geometric sense, but a ball in the sense of metric spaces. If we recognise that in our definition we have the transpose of a vector, multiplied by a positive definite matrix, multiplied by the same vector, we can see that\n", + " \\begin{align*}\n", + " {(\\mu-\\bar{\\mu})}^{T}\\Omega^{-1}(\\mu-\\bar{\\mu}) &=: x^{T}Ax \\\\\n", + " &= x^{T}L^{T}Lx \\quad \\text{(using Cholesky factorisation)} \\\\\n", + " &= {(Lx)}^{T}{(Lx)} \\\\\n", + " &= {\\|Lx\\|}_2^2 \\\\\n", + " &= {\\|x\\|}_{L^{-1}}^2 \\leq \\kappa^2. \n", + " \\end{align*}\n", + " \n", + " In other words, in the norm of $\\Omega$'s Cholesky factor our uncertainty set is a ball of radius $\\kappa$ centred on $\\bar{\\mu}$.\n", + "\n", + "3. Some fourth footnote which has yet to be written. When is, increment and retain. " ] } ],