From c82f19e5df5e80cd03408a80b7a11d47a243b2d5 Mon Sep 17 00:00:00 2001 From: Joshua Fogg Date: Thu, 25 Apr 2024 13:15:49 +0100 Subject: [PATCH] Bypassed Gurobi's issue with `np.inner` --- Example2/n2-test.ipynb | 74 ++++++------- Scratch/numerical-error.ipynb | 16 ++- Scratch/risk-minimization.ipynb | 191 +++++++++++++++++--------------- robust-genetics.ipynb | 176 +++++++++++++++-------------- 4 files changed, 247 insertions(+), 210 deletions(-) diff --git a/Example2/n2-test.ipynb b/Example2/n2-test.ipynb index b93e15f..345b6cd 100644 --- a/Example2/n2-test.ipynb +++ b/Example2/n2-test.ipynb @@ -142,7 +142,7 @@ " 6 1.74999618e+00 1.75000030e+00 0.00e+00 0.00e+00 2.06e-06 0s\n", " 7 1.75000000e+00 1.75000000e+00 0.00e+00 0.00e+00 2.06e-09 0s\n", "\n", - "Barrier solved model in 7 iterations and 0.03 seconds (0.00 work units)\n", + "Barrier solved model in 7 iterations and 0.02 seconds (0.00 work units)\n", "Optimal objective 1.75000000e+00\n", "\n", "Standard: [0.00000001 0.99999999]\n" @@ -197,9 +197,9 @@ "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n", "\n", "Optimize a model with 1 rows, 3 columns and 2 nonzeros\n", - "Model fingerprint: 0xd120428b\n", + "Model fingerprint: 0xea4b60b6\n", "Model has 2 quadratic objective terms\n", - "Model has 2 quadratic constraints\n", + "Model has 1 quadratic constraint\n", "Coefficient statistics:\n", " Matrix range [1e+00, 1e+00]\n", " QMatrix range [1e-01, 4e+00]\n", @@ -208,30 +208,31 @@ " Bounds range [0e+00, 0e+00]\n", " RHS range [1e+00, 1e+00]\n", "Presolve time: 0.01s\n", - "Presolved: 7 rows, 8 columns, 12 nonzeros\n", - "Presolved model has 3 second-order cone constraints\n", + "Presolved: 6 rows, 7 columns, 10 nonzeros\n", + "Presolved model has 2 second-order cone constraints\n", "Ordering time: 0.00s\n", "\n", "Barrier statistics:\n", - " AA' NZ : 1.700e+01\n", - " Factor NZ : 2.800e+01\n", - " Factor Ops : 1.400e+02 (less than 1 second per iteration)\n", + " AA' NZ : 1.500e+01\n", + " Factor NZ : 2.100e+01\n", + " Factor Ops : 9.100e+01 (less than 1 second per iteration)\n", " Threads : 1\n", "\n", " Objective Residual\n", "Iter Primal Dual Primal Dual Compl Time\n", - " 0 -2.28881626e-02 1.22656250e+00 9.39e-01 1.43e-01 3.26e-01 0s\n", - " 1 6.91883948e-01 1.06699842e+00 1.98e-01 1.57e-07 5.84e-02 0s\n", - " 2 6.44324078e-01 7.21385315e-01 2.17e-07 1.73e-13 7.71e-03 0s\n", - " 3 6.67762956e-01 6.69433547e-01 6.21e-09 6.66e-16 1.67e-04 0s\n", - " 4 6.68355745e-01 6.68371340e-01 1.53e-13 2.00e-15 1.56e-06 0s\n", - " 5 6.68366724e-01 6.68367382e-01 3.57e-11 6.05e-13 6.58e-08 0s\n", + " 0 -5.73852130e-02 1.22656250e+00 8.97e-01 1.41e-01 3.51e-01 0s\n", + " 1 5.91715790e-01 9.52821847e-01 2.13e-01 1.51e-07 6.38e-02 0s\n", + " 2 5.15372316e-01 5.96285470e-01 2.34e-07 6.30e-10 8.99e-03 0s\n", + " 3 5.49225634e-01 5.53202935e-01 4.81e-10 1.05e-15 4.42e-04 0s\n", + " 4 5.52127256e-01 5.52382277e-01 1.12e-13 5.55e-16 2.83e-05 0s\n", + " 5 5.52206847e-01 5.52232334e-01 3.52e-12 3.66e-15 2.83e-06 0s\n", + " 6 5.52216328e-01 5.52218236e-01 1.38e-11 2.96e-14 2.12e-07 0s\n", "\n", - "Barrier solved model in 5 iterations and 0.03 seconds (0.00 work units)\n", - "Optimal objective 6.68366724e-01\n", + "Barrier solved model in 6 iterations and 0.03 seconds (0.00 work units)\n", + "Optimal objective 5.52216328e-01\n", "\n", "Standard: [0.00000001 0.99999999]\n", - "Robust: [0.85714285 0.14285715]\n" + "Robust: [0.83155134 0.16844866]\n" ] } ], @@ -253,7 +254,7 @@ "# add sum-to-half constraints to both models\n", "model_rbs.addConstr(np.ones([1,2]) @ w_std == 1, name=\"sum-to-one\")\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**2 >= w_rbs.transpose()@omega@w_rbs, name=\"uncertainty\")\n", "\n", "# solve problem with Gurobi, print alongside standard solution for comparison\n", "model_rbs.optimize()\n", @@ -265,7 +266,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Note that if $\\kappa = 0$, the standard solution not _exactly_ recovered, only accurate to eight decimal places." + "Note that if $\\kappa = 0$, the standard solution not _exactly_ recovered, only accurate to seven decimal places." ] }, { @@ -283,9 +284,9 @@ "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n", "\n", "Optimize a model with 1 rows, 3 columns and 2 nonzeros\n", - "Model fingerprint: 0x6c4430f5\n", + "Model fingerprint: 0xcd1e4542\n", "Model has 2 quadratic objective terms\n", - "Model has 2 quadratic constraints\n", + "Model has 1 quadratic constraint\n", "Coefficient statistics:\n", " Matrix range [1e+00, 1e+00]\n", " QMatrix range [1e-01, 4e+00]\n", @@ -294,31 +295,30 @@ " Bounds range [0e+00, 0e+00]\n", " RHS range [1e+00, 1e+00]\n", "Presolve time: 0.01s\n", - "Presolved: 7 rows, 8 columns, 12 nonzeros\n", - "Presolved model has 3 second-order cone constraints\n", + "Presolved: 6 rows, 7 columns, 10 nonzeros\n", + "Presolved model has 2 second-order cone constraints\n", "Ordering time: 0.00s\n", "\n", "Barrier statistics:\n", - " AA' NZ : 1.700e+01\n", - " Factor NZ : 2.800e+01\n", - " Factor Ops : 1.400e+02 (less than 1 second per iteration)\n", + " AA' NZ : 1.500e+01\n", + " Factor NZ : 2.100e+01\n", + " Factor Ops : 9.100e+01 (less than 1 second per iteration)\n", " Threads : 1\n", "\n", " Objective Residual\n", "Iter Primal Dual Primal Dual Compl Time\n", - " 0 1.22656250e+00 1.22656250e+00 1.04e+00 6.81e-01 5.08e-01 0s\n", - " 1 1.54200502e+00 1.52822126e+00 4.89e-02 1.42e-01 7.24e-02 0s\n", - " 2 1.64909037e+00 1.63852570e+00 5.38e-08 3.21e-02 1.66e-02 0s\n", - " 3 1.74359514e+00 1.75057927e+00 3.23e-10 8.27e-05 7.48e-04 0s\n", - " 4 1.74999141e+00 1.75002072e+00 1.01e-12 9.47e-08 2.99e-06 0s\n", - " 5 1.74999468e+00 1.75000778e+00 3.11e-12 3.00e-08 1.33e-06 0s\n", - " 6 1.74999882e+00 1.75000092e+00 1.06e-11 1.95e-09 2.11e-07 0s\n", + " 0 1.22656250e+00 1.22656250e+00 1.08e+00 6.94e-01 4.95e-01 0s\n", + " 1 1.61470687e+00 1.61496484e+00 5.67e-02 1.55e-01 6.97e-02 0s\n", + " 2 1.67968985e+00 1.70168242e+00 6.24e-08 2.02e-02 1.10e-02 0s\n", + " 3 1.74763426e+00 1.75020504e+00 7.06e-14 2.52e-05 2.97e-04 0s\n", + " 4 1.74999548e+00 1.75001176e+00 1.87e-14 2.84e-08 1.82e-06 0s\n", + " 5 1.74999978e+00 1.75000067e+00 4.92e-13 9.33e-10 9.90e-08 0s\n", "\n", - "Barrier solved model in 6 iterations and 0.03 seconds (0.00 work units)\n", - "Optimal objective 1.74999882e+00\n", + "Barrier solved model in 5 iterations and 0.03 seconds (0.00 work units)\n", + "Optimal objective 1.74999978e+00\n", "\n", "Standard: [0.0000000076 0.9999999924]\n", - "Robust: [0.0000000230 0.9999999770]\n" + "Robust: [0.0000001531 0.9999998469]\n" ] } ], @@ -340,7 +340,7 @@ "# add sum-to-half constraints to both models\n", "model_rbs.addConstr(np.ones([1,2]) @ w_std == 1, name=\"sum-to-one\")\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**2 >= w_rbs.transpose()@omega@w_rbs, name=\"uncertainty\")\n", "\n", "# solve problem with Gurobi, print alongside standard solution for comparison\n", "model_rbs.optimize()\n", diff --git a/Scratch/numerical-error.ipynb b/Scratch/numerical-error.ipynb index 791a356..56a658e 100644 --- a/Scratch/numerical-error.ipynb +++ b/Scratch/numerical-error.ipynb @@ -6,7 +6,16 @@ "source": [ "# Numerical Error\n", "\n", - "The purpose of this notebook is to illustrate the numerical error which comes from setting $\\kappa = 0$." + "The purpose of this notebook is to illustrate the numerical error which comes from setting $\\kappa = 0$ and using \n", + "\n", + "```python3\n", + "model.addConstr(z**2 <= np.inner(w, omega@w), name=\"uncertainty\")\n", + "```\n", + "as the cone constraint rather than\n", + "```python3\n", + "model.addConstr(z**2 <= w.transpose()@omega@w, name=\"uncertainty\")\n", + "```\n", + "This has been fixed in the latest version of all other notebooks." ] }, { @@ -150,7 +159,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.04 seconds (0.00 work units)\n", + "Barrier solved model in 8 iterations and 0.03 seconds (0.00 work units)\n", "Optimal objective -8.12500000e-01\n", "\n", "w = [0.50000 0.37500 0.12500]\n" @@ -287,6 +296,7 @@ "\n", "# add quadratic uncertainty constraint\n", "model.addConstr(z**2 <= np.inner(w, omega@w), name=\"uncertainty\")\n", + "# model.addConstr(z**2 <= w.transpose()@omega@w, name=\"uncertainty\")\n", "# add sub-to-half constraints\n", "model.addConstr(M @ w == m, name=\"sum-to-half\")\n", "# add weight-bound constraints~\n", @@ -302,7 +312,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Note that even though $\\kappa = 0$, we have $z\\neq0$. This suggests that a small amount of numerical error has been introduced as a result of asking Gurobi to optimize $z$ as well as $w$. " + "Note that even though $\\kappa = 0$, we have $z\\neq0$. This suggests that a small amount of numerical error has been introduced." ] } ], diff --git a/Scratch/risk-minimization.ipynb b/Scratch/risk-minimization.ipynb index 337a592..7a6d8c0 100644 --- a/Scratch/risk-minimization.ipynb +++ b/Scratch/risk-minimization.ipynb @@ -302,7 +302,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", "w = [0.50000 0.37500 0.12500]\n" @@ -402,9 +402,9 @@ "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n", "\n", "Optimize a model with 22 rows, 4 columns and 24 nonzeros\n", - "Model fingerprint: 0x5b10d4aa\n", + "Model fingerprint: 0xf7819a35\n", "Model has 3 quadratic objective terms\n", - "Model has 3 quadratic constraints\n", + "Model has 1 quadratic constraint\n", "Coefficient statistics:\n", " Matrix range [1e+00, 1e+00]\n", " QMatrix range [1e-01, 4e+00]\n", @@ -413,32 +413,46 @@ " Bounds range [0e+00, 0e+00]\n", " RHS range [5e-01, 1e+00]\n", "Presolve removed 21 rows and 1 columns\n", - "Presolve time: 0.01s\n", - "Presolved: 9 rows, 8 columns, 14 nonzeros\n", - "Presolved model has 3 second-order cone constraints\n", - "Ordering time: 0.00s\n", "\n", - "Barrier statistics:\n", - " AA' NZ : 2.200e+01\n", - " Factor NZ : 4.500e+01\n", - " Factor Ops : 2.850e+02 (less than 1 second per iteration)\n", - " Threads : 1\n", + "Continuous model is non-convex -- solving as a MIP\n", "\n", - " Objective Residual\n", - "Iter Primal Dual Primal Dual Compl Time\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.24036608e-01\n", - "\n", - "w = [0.50000 0.36395 0.13605],\n", - "z = 0.04810269528692041.\n" + "Presolve removed 21 rows and 1 columns\n", + "Presolve time: 0.00s\n", + "Presolved: 4 rows, 6 columns, 9 nonzeros\n", + "Presolved model has 2 quadratic objective terms\n", + "Presolved model has 1 quadratic constraint(s)\n", + "Presolved model has 2 bilinear constraint(s)\n", + "Variable types: 6 continuous, 0 integer (0 binary)\n", + "Found heuristic solution: objective -1.0491801\n", + "\n", + "Root relaxation: objective -1.095481e+00, 4 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.09548 0 - -1.04918 -1.09548 4.41% - 0s\n", + " 0 0 -1.09548 0 1 -1.04918 -1.09548 4.41% - 0s\n", + " 0 0 -1.08932 0 2 -1.04918 -1.08932 3.83% - 0s\n", + " 0 0 -1.08772 0 2 -1.04918 -1.08772 3.67% - 0s\n", + " 0 0 -1.08664 0 2 -1.04918 -1.08664 3.57% - 0s\n", + " 0 0 -1.08587 0 2 -1.04918 -1.08587 3.50% - 0s\n", + " 0 0 -1.08527 0 1 -1.04918 -1.08527 3.44% - 0s\n", + " 0 0 -1.08481 0 1 -1.04918 -1.08481 3.40% - 0s\n", + " 0 0 -1.08443 0 1 -1.04918 -1.08443 3.36% - 0s\n", + " 0 0 -1.08415 0 2 -1.04918 -1.08415 3.33% - 0s\n", + " 0 0 -1.04960 0 2 -1.04918 -1.04960 0.04% - 0s\n", + " 0 0 -1.04955 0 2 -1.04918 -1.04955 0.04% - 0s\n", + " 0 2 -1.04955 0 2 -1.04918 -1.04955 0.04% - 0s\n", + "\n", + "Explored 3 nodes (22 simplex iterations) in 0.04 seconds (0.00 work units)\n", + "Thread count was 8 (of 8 available processors)\n", + "\n", + "Solution count 1: -1.04918 \n", + "\n", + "Optimal solution found (tolerance 1.00e-04)\n", + "Best objective -1.049180143116e+00, best bound -1.049222787703e+00, gap 0.0041%\n", + "w = [0.50000 0.42869 0.07131],\n", + "z = 0.9928459846964108.\n" ] } ], @@ -456,7 +470,7 @@ "GRB.MINIMIZE)\n", "\n", "# add quadratic uncertainty constraint\n", - "model.addConstr(z**2 <= np.inner(w, omega@w), name=\"uncertainty\")\n", + "model.addConstr(z**2 <= w.transpose()@omega@w, name=\"uncertainty\")\n", "# add sub-to-half constraints\n", "model.addConstr(M @ w == m, name=\"sum-to-half\")\n", "# add weight-bound constraints\n", @@ -474,20 +488,15 @@ "source": [ "We can repeat this experiment with varying $\\kappa$ to see how how different tolerances for uncertainty impact the robust turning point returned.\n", "\n", - "| $\\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.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", + "| $\\kappa$ | $w$ | $z$ | $f(w,z)$ | Gap |\n", + "| ---: | ------------------------: | ------: | -------: | --: |\n", + "| 0.0 | [0.50000 0.37500 0.12500] | 0.0 | -8.12500000e-01 | 0% |\n", + "| 0.5 | [0.50000 0.42869 0.07131] | 0.9928459846964108 | -1.049180143116e+00 | 0.0041% |\n", + "| 1.0 | [0.50000 0.48606 0.01394] | 1.0931753188733573 | -1.309752491846e+00 | 0.0023% |\n", + "| 2.0 | [0.50000 0.50000 0.00000] | 1.1180339796124110 | -1.868033983797e+00 | 0.0005% |\n", + "| 4.0 | [0.50000 0.50000 0.00000] | 1.1180339856336740 | -2.986067972548e+00 | 0.0021%\n", "\n", - "We also note that there are not further changes once $\\kappa$ falls above a certain threshold." + "We can see that in the case $\\kappa = 0$ we return the standard optimization solution, as expected. We also note that there are not further changes once $\\kappa$ falls above a certain threshold, indicating that we have successfully converged to the \"riskiest\" portfolio." ] }, { @@ -619,7 +628,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.02 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", @@ -628,12 +637,12 @@ "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n", "\n", "Optimize a model with 4 rows, 51 columns and 100 nonzeros\n", - "Model fingerprint: 0x74debc94\n", + "Model fingerprint: 0xd47c065e\n", "Model has 1275 quadratic objective terms\n", - "Model has 50 quadratic constraints\n", + "Model has 1 quadratic constraint\n", "Coefficient statistics:\n", " Matrix range [1e+00, 1e+00]\n", - " QMatrix range [7e-06, 1e+00]\n", + " QMatrix range [1e-05, 1e+00]\n", " Objective range [7e-01, 1e+00]\n", " QObjective range [5e-02, 1e+00]\n", " Bounds range [0e+00, 0e+00]\n", @@ -643,66 +652,70 @@ "Continuous model is non-convex -- solving as a MIP\n", "\n", "Presolve removed 2 rows and 0 columns\n", - "Presolve time: 0.01s\n", - "Presolved: 2552 rows, 1327 columns, 7600 nonzeros\n", + "Presolve time: 0.00s\n", + "Presolved: 2503 rows, 1327 columns, 6326 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", + "Found heuristic solution: objective -1.1617295\n", "\n", - "Root relaxation: objective -1.285097e+00, 31 iterations, 0.00 seconds (0.00 work units)\n", + "Root relaxation: objective -3.490954e+00, 43 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 - 0 -0.97636 -1.02501 4.98% - 0s\n", + " 0 0 -3.49095 0 38 -1.16173 -3.49095 200% - 0s\n", + " 0 0 -3.49095 0 42 -1.16173 -3.49095 200% - 0s\n", + " 0 0 -1.58157 0 231 -1.16173 -1.58157 36.1% - 0s\n", + " 0 0 -1.34540 0 512 -1.16173 -1.34540 15.8% - 0s\n", + " 0 0 -1.32278 0 619 -1.16173 -1.32278 13.9% - 0s\n", + " 0 0 -1.32278 0 59 -1.16173 -1.32278 13.9% - 1s\n", + " 0 0 -1.32278 0 59 -1.16173 -1.32278 13.9% - 1s\n", + " 0 2 -1.32278 0 59 -1.16173 -1.32278 13.9% - 1s\n", + " 110 118 -1.23808 14 191 -1.16173 -1.27940 10.1% 343 5s\n", + " 340 353 -1.19438 40 319 -1.16173 -1.27940 10.1% 377 10s\n", + "H 551 573 -1.1617295 -1.27240 9.53% 438 12s\n", "\n", - "Cutting planes:\n", - " RLT: 6\n", - "\n", - "Explored 1 nodes (1828 simplex iterations) in 0.22 seconds (0.10 work units)\n", + "Explored 600 nodes (256270 simplex iterations) in 13.46 seconds (12.04 work units)\n", "Thread count was 8 (of 8 available processors)\n", "\n", - "Solution count 1: -0.976358 \n", + "Solution count 1: -1.16173 \n", "\n", "Optimal solution found (tolerance 5.00e-02)\n", - "Best objective -9.763583540500e-01, best bound -1.025012483088e+00, gap 4.9832%\n", + "Best objective -1.161729466120e+00, best bound -1.217554010902e+00, gap 4.8053%\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", - "\n", - "Maximum change: 0.00209\n", - "Average change: 0.00029\n", + "00 0.00000 0.00000 01 0.08319 0.14976\n", + "02 0.11649 0.49337 03 0.00000 0.00000\n", + "04 0.00000 0.00000 05 0.00000 0.00000\n", + "06 0.00000 0.00000 07 0.06817 0.06251\n", + "08 0.00000 0.00000 09 0.03762 0.00000\n", + "10 0.00000 0.00000 11 0.00000 0.00000\n", + "12 0.00000 0.00000 13 0.00355 0.00000\n", + "14 0.11542 0.00518 15 0.06966 0.00000\n", + "16 0.00000 0.00000 17 0.00000 0.00000\n", + "18 0.00000 0.00000 19 0.00000 0.00000\n", + "20 0.00000 0.00000 21 0.00001 0.00000\n", + "22 0.00000 0.00000 23 0.00000 0.00000\n", + "24 0.03537 0.00000 25 0.00000 0.00000\n", + "26 0.00000 0.00000 27 0.00000 0.00000\n", + "28 0.00000 0.00000 29 0.00000 0.00000\n", + "30 0.00000 0.00000 31 0.03580 0.00000\n", + "32 0.08447 0.00000 33 0.00000 0.00000\n", + "34 0.00000 0.00000 35 0.02849 0.09487\n", + "36 0.13489 0.00144 37 0.07446 0.10716\n", + "38 0.00000 0.00000 39 0.00000 0.00000\n", + "40 0.00000 0.00000 41 0.07312 0.08570\n", + "42 0.00000 0.00000 43 0.00000 0.00000\n", + "44 0.01336 0.00000 45 0.00000 0.00000\n", + "46 0.00000 0.00000 47 0.02593 0.00000\n", + "48 0.00000 0.00000 49 0.00000 0.00000\n", + "\n", + "Maximum change: 0.37689\n", + "Average change: 0.02220\n", "Minimum change: 0.00000\n" ] } @@ -750,7 +763,7 @@ "model_rbs.addConstr(M @ w_rbs == 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**2 <= w_rbs.transpose()@omega@w_rbs, name=\"uncertainty\")\n", "\n", "# since working with non-trivial size, set a time limit\n", "time_limit = 60*5 # 5 minutes\n", @@ -783,7 +796,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "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." + "We can see that with parameters $\\lambda = 0.5, \\kappa = 2$, there's a more significant shift in the portfolios produced and we see in the return maximization version of the problem. It also takes Gurobi significantly longer to compute. More testing would be needed to confirm why this was the case." ] } ], diff --git a/robust-genetics.ipynb b/robust-genetics.ipynb index d22db74..1507c41 100644 --- a/robust-genetics.ipynb +++ b/robust-genetics.ipynb @@ -339,7 +339,7 @@ " Bounds range [0e+00, 0e+00]\n", " RHS range [5e-01, 1e+00]\n", "Presolve removed 21 rows and 1 columns\n", - "Presolve time: 0.01s\n", + "Presolve time: 0.00s\n", "Presolved: 1 rows, 2 columns, 2 nonzeros\n", "Presolved model has 2 quadratic objective terms\n", "Ordering time: 0.00s\n", @@ -362,7 +362,7 @@ " 7 1.62500000e+00 1.62501541e+00 0.00e+00 5.55e-17 4.81e-07 0s\n", " 8 1.62500000e+00 1.62500002e+00 0.00e+00 0.00e+00 4.81e-10 0s\n", "\n", - "Barrier solved model in 8 iterations and 0.02 seconds (0.00 work units)\n", + "Barrier solved model in 8 iterations and 0.01 seconds (0.00 work units)\n", "Optimal objective 1.62500000e+00\n", "\n", "w = [0.50000 0.37500 0.12500]\n" @@ -466,9 +466,9 @@ "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n", "\n", "Optimize a model with 22 rows, 4 columns and 24 nonzeros\n", - "Model fingerprint: 0x5341cd5b\n", + "Model fingerprint: 0xe9e68b83\n", "Model has 3 quadratic objective terms\n", - "Model has 3 quadratic constraints\n", + "Model has 1 quadratic constraint\n", "Coefficient statistics:\n", " Matrix range [1e+00, 1e+00]\n", " QMatrix range [1e-01, 4e+00]\n", @@ -478,32 +478,30 @@ " RHS range [5e-01, 1e+00]\n", "Presolve removed 21 rows and 1 columns\n", "Presolve time: 0.01s\n", - "Presolved: 6 rows, 6 columns, 9 nonzeros\n", + "Presolved: 7 rows, 8 columns, 11 nonzeros\n", "Presolved model has 2 second-order cone constraints\n", "Ordering time: 0.00s\n", "\n", "Barrier statistics:\n", - " AA' NZ : 1.100e+01\n", - " Factor NZ : 2.100e+01\n", - " Factor Ops : 9.100e+01 (less than 1 second per iteration)\n", + " AA' NZ : 1.700e+01\n", + " Factor NZ : 2.800e+01\n", + " Factor Ops : 1.400e+02 (less than 1 second per iteration)\n", " Threads : 1\n", "\n", " Objective Residual\n", "Iter Primal Dual Primal Dual Compl Time\n", - " 0 6.75206688e-01 1.66071429e+00 1.95e+00 7.60e-01 1.40e+00 0s\n", - " 1 9.18871990e-01 2.22754663e+00 2.08e-01 8.36e-07 1.99e-01 0s\n", - " 2 1.26521618e+00 1.51422516e+00 1.85e-02 9.20e-13 3.06e-02 0s\n", - " 3 1.27457950e+00 1.28458784e+00 2.03e-08 4.44e-16 1.11e-03 0s\n", - " 4 1.28096283e+00 1.28138761e+00 8.88e-13 6.66e-16 4.72e-05 0s\n", - " 5 1.28122931e+00 1.28125496e+00 3.34e-12 5.11e-15 2.85e-06 0s\n", - " 6 1.28124988e+00 1.28125699e+00 1.16e-10 2.63e-13 7.90e-07 0s\n", - " 7 1.28125000e+00 1.28125001e+00 1.41e-09 4.38e-12 9.12e-10 0s\n", + " 0 6.73751896e-01 1.79000000e+00 1.90e+00 6.81e-01 1.00e+00 0s\n", + " 1 1.20057128e+00 2.19855901e+00 1.99e-01 7.50e-07 1.45e-01 0s\n", + " 2 1.11075775e+00 1.33966849e+00 2.19e-07 8.25e-13 2.29e-02 0s\n", + " 3 1.17894103e+00 1.20876338e+00 3.10e-13 2.22e-16 2.98e-03 0s\n", + " 4 1.19359138e+00 1.19386663e+00 1.07e-14 2.22e-16 2.75e-05 0s\n", + " 5 1.19381906e+00 1.19382014e+00 3.32e-12 6.66e-16 1.08e-07 0s\n", "\n", - "Barrier solved model in 7 iterations and 0.02 seconds (0.00 work units)\n", - "Optimal objective 1.28125000e+00\n", + "Barrier solved model in 5 iterations and 0.02 seconds (0.00 work units)\n", + "Optimal objective 1.19381906e+00\n", "\n", - "w = [0.50000 0.31249 0.18751],\n", - "z = 0.62499.\n" + "w = [0.50000 0.32623 0.17377],\n", + "z = 0.82431.\n" ] } ], @@ -520,11 +518,11 @@ " w.transpose()@expected_breeding_values - (lam/2)*w@(relationship_matrix@w) - kappa*z,\n", "GRB.MAXIMIZE)\n", "\n", - "# add quadratic uncertainty constraint\n", - "model.addConstr(z**2 >= np.inner(w, omega@w), name=\"uncertainty\")\n", + "# add quadratic uncertainty constraintāµ\n", + "model.addConstr(z**2 >= w.transpose()@omega@w, name=\"uncertainty\")\n", "# add sub-to-half constraints\n", "model.addConstr(M @ w == m, name=\"sum-to-half\")\n", - "# add weight-bound constraints~\n", + "# add weight-bound constraints\n", "model.addConstr(w >= lower_bound, name=\"lower-bound\")\n", "model.addConstr(w <= upper_bound, name=\"upper-bound\")\n", "\n", @@ -541,17 +539,21 @@ "\n", "| $\\kappa$ | $w$ | $z$ | $f(w,z)$ |\n", "| -------: | ------------------------: | ------: | -------: |\n", - "| 0.00 | [0.50000 0.37500 0.12500] | 3.47856 | 1.62500 |\n", - "| 0.25 | [0.50000 0.34374 0.15626] | 0.68748 | 1.44531 |\n", - "| 0.50 | [0.50000 0.31249 0.18751] | 0.62499 | 1.28125 |\n", - "| 0.75 | [0.50000 0.28140 0.21860] | 0.56280 | 1.13281 |\n", - "| 1.00 | [0.50000 0.25032 0.24968] | 0.50063 | 1.00000 |\n", - "| 2.00 | [0.50000 0.25000 0.25000] | 0.50000 | 0.50000 |\n", - "| 4.00 | [0.50000 0.25000 0.25000] | 0.50000 | -0.50000 |\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$. We also note that there are not further changes once $\\kappa$ falls above a certain threshold.\n", - "\n", - "TODO _do a comparison of the value $z$ against $\\sqrt{w\\Omega w}$, analyzing the gap between the two._" + "| 0.00 | [0.50000 0.37500 0.12500] | 3.52875 | 1.62499880e+00\n", + "| 0.25 | [0.50000 0.34998 0.15002] | 0.86184 | 1.40453330e+00\n", + "| 0.50 | [0.50000 0.32623 0.17377] | 0.82431 | 1.19381906e+00\n", + "| 0.75 | [0.50000 0.30443 0.19557] | 0.79089 | 9.91998750e-01\n", + "| 1.00 | [0.50000 0.28391 0.21609] | 0.75044 | 7.98188968e-01\n", + "| 2.00 | [0.50000 0.21877 0.28123] | 0.67181 | 8.61190139e-02\n", + "| 4.00 | [0.50000 0.14644 0.35356] | 0.59279 | -1.16409198e+00\n", + "| 8.00 | [0.50000 0.09102 0.40898] | 0.55140 | -3.43139120e+00\n", + "| 16.00 | [0.50000 0.05658 0.44342] | 0.53608 | -7.76342696e+00\n", + "| 32.00 | [0.50000 0.03695 0.46305] | 0.53128 | -1.62903237e+01\n", + "| 64.00 | [0.50000 0.02635 0.47365] | 0.52992 | -3.32626440e+01\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$. We also not that as $\\kappa\\rightarrow\\infty$, we tend towards the portfolio of maximum risk.\n", + "\n", + "TODO _do a comparison of the value $z$ against $\\sqrt{w\\Omega w}$, analyzing the gap between the two. I suspect this will explain away what I originally thought was \"numerical error\"_" ] }, { @@ -694,7 +696,7 @@ " 10 2.03643494e+00 2.03643500e+00 6.20e-13 2.22e-15 1.31e-09 0s\n", " 11 2.03643499e+00 2.03643499e+00 1.10e-12 1.33e-15 8.34e-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.02 seconds (0.00 work units)\n", "Optimal objective 2.03643499e+00\n", "\n", "Gurobi Optimizer version 11.0.0 build v11.0.0rc2 (linux64 - \"Ubuntu 22.04.4 LTS\")\n", @@ -703,76 +705,75 @@ "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n", "\n", "Optimize a model with 4 rows, 51 columns and 100 nonzeros\n", - "Model fingerprint: 0xacf3991a\n", + "Model fingerprint: 0x806a4c54\n", "Model has 1275 quadratic objective terms\n", - "Model has 50 quadratic constraints\n", + "Model has 1 quadratic constraint\n", "Coefficient statistics:\n", " Matrix range [1e+00, 1e+00]\n", - " QMatrix range [7e-06, 1e+00]\n", + " QMatrix range [1e-05, 1e+00]\n", " Objective range [1e+00, 1e+01]\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 removed 2 rows and 1 columns\n", "\n", "Continuous model is non-convex -- solving as a MIP\n", "\n", - "Presolve removed 2 rows and 0 columns\n", - "Presolve time: 0.01s\n", - "Presolved: 2552 rows, 1327 columns, 7600 nonzeros\n", + "Presolve removed 2 rows and 1 columns\n", + "Presolve time: 0.00s\n", + "Presolved: 2503 rows, 1325 columns, 6325 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 2.0353486\n", + "Variable types: 1325 continuous, 0 integer (0 binary)\n", + "Found heuristic solution: objective 2.0364350\n", "\n", - "Root relaxation: interrupted, 25 iterations, 0.00 seconds (0.00 work units)\n", + "Root relaxation: interrupted, 37 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 - 0 2.03535 2.03643 0.05% - 0s\n", + " 0 0 - 0 2.03643 2.03643 0.00% - 0s\n", "\n", - "Explored 1 nodes (25 simplex iterations) in 0.07 seconds (0.03 work units)\n", + "Explored 1 nodes (37 simplex iterations) in 0.03 seconds (0.01 work units)\n", "Thread count was 8 (of 8 available processors)\n", "\n", - "Solution count 1: 2.03535 \n", + "Solution count 1: 2.03643 \n", "\n", "Optimal solution found (tolerance 5.00e-02)\n", - "Best objective 2.035348603347e+00, best bound 2.036434994087e+00, gap 0.0534%\n", + "Best objective 2.036434950287e+00, best bound 2.036434994087e+00, gap 0.0000%\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.00010 01 0.10276 0.10148\n", - "02 0.14509 0.14457 03 0.00000 0.00012\n", - "04 0.00000 0.00005 05 0.00000 0.00020\n", - "06 0.00000 0.00012 07 0.10172 0.10016\n", - "08 0.00000 0.00017 09 0.02222 0.02115\n", - "10 0.00000 0.00011 11 0.00000 0.00009\n", - "12 0.00000 0.00010 13 0.00000 0.00059\n", - "14 0.11196 0.11075 15 0.07168 0.07042\n", - "16 0.00000 0.00004 17 0.00000 0.00006\n", - "18 0.00000 0.00008 19 0.00000 0.00010\n", - "20 0.00000 0.00009 21 0.00000 0.00038\n", - "22 0.00000 0.00011 23 0.00000 0.00019\n", - "24 0.00000 0.00093 25 0.00000 0.00023\n", - "26 0.00000 0.00011 27 0.00000 0.00012\n", - "28 0.00000 0.00006 29 0.00000 0.00014\n", - "30 0.00000 0.00011 31 0.02845 0.02737\n", - "32 0.07177 0.07175 33 0.00000 0.00039\n", - "34 0.00000 0.00012 35 0.00000 0.00335\n", - "36 0.17118 0.16970 37 0.08277 0.07995\n", - "38 0.00000 0.00007 39 0.00000 0.00016\n", - "40 0.00000 0.00010 41 0.09038 0.08832\n", - "42 0.00000 0.00009 43 0.00000 0.00018\n", - "44 0.00000 0.00035 45 0.00000 0.00008\n", - "46 0.00000 0.00020 47 0.00001 0.00472\n", - "48 0.00000 0.00011 49 0.00000 0.00005\n", + "00 0.00000 0.00000 01 0.10276 0.10276\n", + "02 0.14509 0.14509 03 0.00000 0.00000\n", + "04 0.00000 0.00000 05 0.00000 0.00000\n", + "06 0.00000 0.00000 07 0.10172 0.10172\n", + "08 0.00000 0.00000 09 0.02222 0.02222\n", + "10 0.00000 0.00000 11 0.00000 0.00000\n", + "12 0.00000 0.00000 13 0.00000 0.00000\n", + "14 0.11196 0.11196 15 0.07168 0.07168\n", + "16 0.00000 0.00000 17 0.00000 0.00000\n", + "18 0.00000 0.00000 19 0.00000 0.00000\n", + "20 0.00000 0.00000 21 0.00000 0.00000\n", + "22 0.00000 0.00000 23 0.00000 0.00000\n", + "24 0.00000 0.00000 25 0.00000 0.00000\n", + "26 0.00000 0.00000 27 0.00000 0.00000\n", + "28 0.00000 0.00000 29 0.00000 0.00000\n", + "30 0.00000 0.00000 31 0.02845 0.02845\n", + "32 0.07177 0.07177 33 0.00000 0.00000\n", + "34 0.00000 0.00000 35 0.00000 0.00000\n", + "36 0.17118 0.17118 37 0.08277 0.08277\n", + "38 0.00000 0.00000 39 0.00000 0.00000\n", + "40 0.00000 0.00000 41 0.09038 0.09038\n", + "42 0.00000 0.00000 43 0.00000 0.00000\n", + "44 0.00000 0.00000 45 0.00000 0.00000\n", + "46 0.00000 0.00000 47 0.00001 0.00002\n", + "48 0.00000 0.00000 49 0.00000 0.00000\n", "\n", - "Maximum change: 0.00471\n", - "Average change: 0.00058\n", - "Minimum change: 0.00002\n" + "Maximum change: 0.00001\n", + "Average change: 0.00000\n", + "Minimum change: 0.00000\n" ] } ], @@ -820,7 +821,7 @@ "model_rbs.addConstr(M @ w_rbs == 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**2 <= w_rbs.transpose()@omega@w_rbs, name=\"uncertainty\")\n", "\n", "# since working with non-trivial size, set a time limit\n", "time_limit = 60*5 # 5 minutes\n", @@ -857,7 +858,9 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We can see that with parameters $\\lambda = 1, \\kappa = 10$, we go from having most candidates having zero contribution in the standard solution to most candidates having a contribution of order 0.0001 in the robust solution." + "We can see that with parameters $\\lambda = 1, \\kappa = 10$, we go from having most candidates having zero contribution in the standard solution to most candidates having a contribution of order 0.0001 in the robust solution.\n", + "\n", + "BUG: _with issue regarding Gurobi and `np.inner` fixed, there's now little-if-any difference between the standard and robust solutions. Need to look at what's going on here._" ] }, { @@ -893,7 +896,18 @@ " &= {\\|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}$." + " 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", + "5.\n", + " Due to issues with Gurobi, it's important that our uncertainty constraint is modelled as\n", + " ```python3\n", + " model.addConstr(z**2 <= w.transpose()@omega@w, name=\"uncertainty\")\n", + " ```\n", + " rather than something like\n", + " ```python3\n", + " model.addConstr(z**2 <= np.inner(w, omega@w), name=\"uncertainty\")\n", + " ```\n", + " For reasons unknown, Gurobi treats these two constraints differently." ] } ],