Skip to content

Instantly share code, notes, and snippets.

@andyfaff

andyfaff/cobyqa.ipynb

Created May 16, 2024 22:33
Show Gist options
  • Save andyfaff/4cfdc554f9f520fb1a6a6f983025996f to your computer and use it in GitHub Desktop.
Save andyfaff/4cfdc554f9f520fb1a6a6f983025996f to your computer and use it in GitHub Desktop.
Download ZIP
Display the source blob
Display the rendered blob
Raw
cobyqa.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"id": "ef22a576-f902-4da1-a785-1e6b522c634d",
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"%load_ext line_profiler\n",
"import numpy as np\n",
"from scipy.optimize import rosen, differential_evolution, minimize\n",
"from scipy.optimize._cobyqa_py import _minimize_cobyqa\n",
"from scipy._lib.cobyqa import minimize as cby\n",
"from scipy._lib.cobyqa.problem import Problem\n",
"from scipy._lib.cobyqa.models import Models, Quadratic\n",
"from scipy._lib.cobyqa.framework import TrustRegion\n",
"from scipy._lib.cobyqa.subsolvers.geometry import spider_geometry\n",
"import time"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "85d20172-0bec-48f5-a146-6a0e152b6e1e",
"metadata": {},
"outputs": [],
"source": [
"bounds = [(0, 10)] * 10\n",
"rng = np.random.default_rng(1)\n",
"x0 = rng.random(size=10)*10"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "32a783da-1e87-4744-a581-f4a4414515b6",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(5,)"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a = np.arange(50).reshape(10, 5)\n",
"np.linalg.norm(a, axis=0).shape"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "e929aeea-3519-4a78-94d3-c367f9967f36",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.0028153699710471362 1093 1.8765268325805664\n"
]
}
],
"source": [
"start = time.time()\n",
"res = minimize(rosen, x0, method='cobyqa')\n",
"finish = time.time()\n",
"sftime = time.time()\n",
"_ = [rosen(res.x) for i in range(res.nfev)]\n",
"fftime = time.time()\n",
"\n",
"print((fftime - sftime)/(finish - start), res.nfev, finish-start)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "34a2a933-0895-49e0-9a31-aa394b0933d0",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.3194187963354983 1898 0.026467084884643555\n"
]
}
],
"source": [
"start = time.time()\n",
"res = minimize(rosen, x0, method='nelder-mead')\n",
"finish = time.time()\n",
"sftime = time.time()\n",
"_ = [rosen(res.x) for i in range(res.nfev)]\n",
"fftime = time.time()\n",
"\n",
"print((fftime - sftime)/(finish - start), res.nfev, finish-start)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "140544e4-a7c1-4058-83fe-77afbe92594c",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Timer unit: 1e-09 s\n",
"\n",
"Total time: 0.882836 s\n",
"File: /Users/andrew/Documents/Andy/programming/scipy/build-install/lib/python3.12/site-packages/scipy/_lib/cobyqa/framework.py\n",
"Function: get_geometry_step at line 634\n",
"\n",
"Line # Hits Time Per Hit % Time Line Contents\n",
"==============================================================\n",
" 634 def get_geometry_step(self, k_new, options):\n",
" 635 \"\"\"\n",
" 636 Get the geometry-improving step.\n",
" 637 \n",
" 638 Three different geometry-improving steps are computed and the best one\n",
" 639 is returned. For more details, see Section 5.2.7 of [1]_.\n",
" 640 \n",
" 641 Parameters\n",
" 642 ----------\n",
" 643 k_new : int\n",
" 644 Index of the interpolation point to be modified.\n",
" 645 options : dict\n",
" 646 Options of the solver.\n",
" 647 \n",
" 648 Returns\n",
" 649 -------\n",
" 650 `numpy.ndarray`, shape (n,)\n",
" 651 Geometry-improving step.\n",
" 652 \n",
" 653 Raises\n",
" 654 ------\n",
" 655 `numpy.linalg.LinAlgError`\n",
" 656 If the computation of a determinant fails.\n",
" 657 \n",
" 658 References\n",
" 659 ----------\n",
" 660 .. [1] T. M. Ragonneau. *Model-Based Derivative-Free Optimization\n",
" 661 Methods and Software*. PhD thesis, Department of Applied\n",
" 662 Mathematics, The Hong Kong Polytechnic University, Hong Kong, China,\n",
" 663 2022. URL: https://theses.lib.polyu.edu.hk/handle/200/12294.\n",
" 664 \"\"\"\n",
" 665 327 95000.0 290.5 0.0 if options[Options.DEBUG]:\n",
" 666 assert k_new != self.best_index, \\\n",
" 667 \"The index `k_new` must be different from the best index.\"\n",
" 668 \n",
" 669 # Build the k_new-th Lagrange polynomial.\n",
" 670 327 1744000.0 5333.3 0.2 coord_vec = np.squeeze(np.eye(1, self.models.npt, k_new))\n",
" 671 654 117963000.0 180371.6 13.4 lag = Quadratic(\n",
" 672 327 144000.0 440.4 0.0 self.models.interpolation,\n",
" 673 327 43000.0 131.5 0.0 coord_vec,\n",
" 674 327 57000.0 174.3 0.0 options[Options.DEBUG],\n",
" 675 )\n",
" 676 327 2632000.0 8048.9 0.3 g_lag = lag.grad(self.x_best, self.models.interpolation)\n",
" 677 \n",
" 678 # Compute a simple constrained Cauchy step.\n",
" 679 327 934000.0 2856.3 0.1 xl = self._pb.bounds.xl - self.x_best\n",
" 680 327 902000.0 2758.4 0.1 xu = self._pb.bounds.xu - self.x_best\n",
" 681 654 47588000.0 72764.5 5.4 step = cauchy_geometry(\n",
" 682 327 41000.0 125.4 0.0 0.0,\n",
" 683 327 42000.0 128.4 0.0 g_lag,\n",
" 684 327 68000.0 208.0 0.0 lambda v: lag.curv(v, self.models.interpolation),\n",
" 685 327 26000.0 79.5 0.0 xl,\n",
" 686 327 36000.0 110.1 0.0 xu,\n",
" 687 327 109000.0 333.3 0.0 self.radius,\n",
" 688 327 67000.0 204.9 0.0 options[Options.DEBUG],\n",
" 689 )\n",
" 690 327 222675000.0 680963.3 25.2 sigma = self.models.determinants(self.x_best + step, k_new)\n",
" 691 \n",
" 692 # Compute the solution on the straight lines joining the interpolation\n",
" 693 # points to the k-th one, and choose it if it provides a larger value\n",
" 694 # of the determinant of the interpolation system in absolute value.\n",
" 695 327 42000.0 128.4 0.0 xpt = (\n",
" 696 654 743000.0 1136.1 0.1 self.models.interpolation.xpt\n",
" 697 327 367000.0 1122.3 0.0 - self.models.interpolation.xpt[:, self.best_index, np.newaxis]\n",
" 698 )\n",
" 699 327 1882000.0 5755.4 0.2 xpt[:, [0, self.best_index]] = xpt[:, [self.best_index, 0]]\n",
" 700 654 258084000.0 394623.9 29.2 step_alt = spider_geometry(\n",
" 701 327 31000.0 94.8 0.0 0.0,\n",
" 702 327 43000.0 131.5 0.0 g_lag,\n",
" 703 327 76000.0 232.4 0.0 lambda v: lag.curv(v, self.models.interpolation),\n",
" 704 327 111000.0 339.4 0.0 xpt[:, 1:],\n",
" 705 327 30000.0 91.7 0.0 xl,\n",
" 706 327 39000.0 119.3 0.0 xu,\n",
" 707 327 96000.0 293.6 0.0 self.radius,\n",
" 708 327 70000.0 214.1 0.0 options[Options.DEBUG],\n",
" 709 )\n",
" 710 327 222076000.0 679131.5 25.2 sigma_alt = self.models.determinants(self.x_best + step_alt, k_new)\n",
" 711 327 120000.0 367.0 0.0 if abs(sigma_alt) > abs(sigma):\n",
" 712 15 2000.0 133.3 0.0 step = step_alt\n",
" 713 15 2000.0 133.3 0.0 sigma = sigma_alt\n",
" 714 \n",
" 715 # Compute a Cauchy step on the tangent space of the active constraints.\n",
" 716 327 2799000.0 8559.6 0.3 if self._pb.type in [\n",
" 717 \"linearly constrained\",\n",
" 718 \"nonlinearly constrained\",\n",
" 719 ]:\n",
" 720 aub, bub, aeq, beq = self.get_constraint_linearizations(\n",
" 721 self.x_best\n",
" 722 )\n",
" 723 tol_bd = get_arrays_tol(xl, xu)\n",
" 724 tol_ub = get_arrays_tol(bub)\n",
" 725 free_xl = xl <= -tol_bd\n",
" 726 free_xu = xu >= tol_bd\n",
" 727 free_ub = bub >= tol_ub\n",
" 728 \n",
" 729 # Compute the Cauchy step.\n",
" 730 n_act, q = qr_tangential_byrd_omojokun(\n",
" 731 aub,\n",
" 732 aeq,\n",
" 733 free_xl,\n",
" 734 free_xu,\n",
" 735 free_ub,\n",
" 736 )\n",
" 737 g_lag_proj = q[:, n_act:] @ (q[:, n_act:].T @ g_lag)\n",
" 738 norm_g_lag_proj = np.linalg.norm(g_lag_proj)\n",
" 739 if (\n",
" 740 0 < n_act < self._pb.n\n",
" 741 and norm_g_lag_proj > TINY * self.radius\n",
" 742 ):\n",
" 743 step_alt = (self.radius / norm_g_lag_proj) * g_lag_proj\n",
" 744 if lag.curv(step_alt, self.models.interpolation) < 0.0:\n",
" 745 step_alt = -step_alt\n",
" 746 \n",
" 747 # Evaluate the constraint violation at the Cauchy step.\n",
" 748 cbd = np.block([xl - step_alt, step_alt - xu])\n",
" 749 cub = aub @ step_alt - bub\n",
" 750 ceq = aeq @ step_alt - beq\n",
" 751 maxcv_val = max(\n",
" 752 np.max(array, initial=0.0)\n",
" 753 for array in [cbd, cub, np.abs(ceq)]\n",
" 754 )\n",
" 755 \n",
" 756 # Accept the new step if it is nearly feasible and do not\n",
" 757 # drastically worsen the determinant of the interpolation\n",
" 758 # system in absolute value.\n",
" 759 tol = np.max(np.abs(step_alt[~free_xl]), initial=0.0)\n",
" 760 tol = np.max(np.abs(step_alt[~free_xu]), initial=tol)\n",
" 761 tol = np.max(np.abs(aub[~free_ub, :] @ step_alt), initial=tol)\n",
" 762 tol = min(10.0 * tol, 1e-2 * np.linalg.norm(step_alt))\n",
" 763 if maxcv_val <= tol:\n",
" 764 sigma_alt = self.models.determinants(\n",
" 765 self.x_best + step_alt,\n",
" 766 k_new\n",
" 767 )\n",
" 768 if abs(sigma_alt) >= 0.1 * abs(sigma):\n",
" 769 step = np.clip(step_alt, xl, xu)\n",
" 770 \n",
" 771 327 71000.0 217.1 0.0 if options[Options.DEBUG]:\n",
" 772 tol = get_arrays_tol(xl, xu)\n",
" 773 if np.any(step + tol < xl) or np.any(xu < step - tol):\n",
" 774 warnings.warn(\n",
" 775 \"The geometry step does not respect the bound \"\n",
" 776 \"constraints.\",\n",
" 777 RuntimeWarning,\n",
" 778 2,\n",
" 779 )\n",
" 780 if np.linalg.norm(step) > 1.1 * self.radius:\n",
" 781 warnings.warn(\n",
" 782 \"The geometry step does not respect the \"\n",
" 783 \"trust-region constraint.\",\n",
" 784 RuntimeWarning,\n",
" 785 2,\n",
" 786 )\n",
" 787 327 986000.0 3015.3 0.1 return step\n",
"\n",
"Total time: 0.03948 s\n",
"File: /Users/andrew/Documents/Andy/programming/scipy/build-install/lib/python3.12/site-packages/scipy/_lib/cobyqa/models.py\n",
"Function: build_system at line 450\n",
"\n",
"Line # Hits Time Per Hit % Time Line Contents\n",
"==============================================================\n",
" 450 @staticmethod\n",
" 451 def build_system(xpt):\n",
" 452 \"\"\"\n",
" 453 Build the left-hand side matrix of the interpolation system.\n",
" 454 \n",
" 455 Parameters\n",
" 456 ----------\n",
" 457 xpt : `numpy.ndarray`, shape (n, npt)\n",
" 458 Interpolation points.\n",
" 459 \n",
" 460 Returns\n",
" 461 -------\n",
" 462 `numpy.ndarray`, shape (npt + n + 1, npt + n + 1)\n",
" 463 Left-hand side matrix of the interpolation system.\n",
" 464 \"\"\"\n",
" 465 4375 818000.0 187.0 2.1 n, npt = xpt.shape\n",
" 466 4375 5117000.0 1169.6 13.0 a = np.zeros((npt + n + 1, npt + n + 1))\n",
" 467 4375 19135000.0 4373.7 48.5 a[:npt, :npt] = 0.5 * (xpt.T @ xpt) ** 2.0\n",
" 468 4375 2721000.0 621.9 6.9 a[:npt, npt] = 1.0\n",
" 469 4375 3687000.0 842.7 9.3 a[:npt, npt + 1:] = xpt.T\n",
" 470 4375 2389000.0 546.1 6.1 a[npt, :npt] = 1.0\n",
" 471 4375 3122000.0 713.6 7.9 a[npt + 1:, :npt] = xpt\n",
" 472 4375 2491000.0 569.4 6.3 return a\n",
"\n",
"Total time: 1.42702 s\n",
"File: /Users/andrew/Documents/Andy/programming/scipy/build-install/lib/python3.12/site-packages/scipy/_lib/cobyqa/models.py\n",
"Function: solve_systems at line 474\n",
"\n",
"Line # Hits Time Per Hit % Time Line Contents\n",
"==============================================================\n",
" 474 @staticmethod\n",
" 475 def solve_systems(interpolation, rhs):\n",
" 476 \"\"\"\n",
" 477 Solve the interpolation systems.\n",
" 478 \n",
" 479 Parameters\n",
" 480 ----------\n",
" 481 interpolation : `cobyqa.models.Interpolation`\n",
" 482 Interpolation set.\n",
" 483 rhs : `numpy.ndarray`, shape (npt + n + 1, m)\n",
" 484 Right-hand side vectors of the ``m`` interpolation systems.\n",
" 485 \n",
" 486 Returns\n",
" 487 -------\n",
" 488 `numpy.ndarray`, shape (npt + n + 1, m)\n",
" 489 Solutions of the interpolation systems.\n",
" 490 `numpy.ndarray`, shape (m, )\n",
" 491 Whether the interpolation systems are ill-conditioned.\n",
" 492 \n",
" 493 Raises\n",
" 494 ------\n",
" 495 `numpy.linalg.LinAlgError`\n",
" 496 If the interpolation systems are ill-defined.\n",
" 497 \"\"\"\n",
" 498 4375 1546000.0 353.4 0.1 n, npt = interpolation.xpt.shape\n",
" 499 4375 1136000.0 259.7 0.1 assert rhs.ndim == 2 and rhs.shape[0] == npt + n + 1, \\\n",
" 500 \"The shape of `rhs` is not valid.\"\n",
" 501 \n",
" 502 # Compute the scaled directions from the base point to the\n",
" 503 # interpolation points. We scale the directions to avoid numerical\n",
" 504 # difficulties.\n",
" 505 \n",
" 506 4375 107293000.0 24524.1 7.5 scale = np.max(np.linalg.norm(interpolation.xpt, axis=0), initial=EPS)\n",
" 507 4375 4616000.0 1055.1 0.3 xpt_scale = interpolation.xpt / scale\n",
" 508 \n",
" 509 # Build the left-hand side matrix of the interpolation system. The\n",
" 510 # matrix below stores diag(left_scaling) * W * diag(right_scaling),\n",
" 511 # where W is the theoretical matrix of the interpolation system. The\n",
" 512 # left and right scaling matrices are chosen to keep the elements in\n",
" 513 # the matrix well-balanced.\n",
" 514 4375 47907000.0 10950.2 3.4 a = Quadratic.build_system(xpt_scale)\n",
" 515 \n",
" 516 # Build the left and right scaling diagonal matrices.\n",
" 517 4375 1485000.0 339.4 0.1 left_scaling = np.empty(npt + n + 1)\n",
" 518 4375 1223000.0 279.5 0.1 right_scaling = np.empty(npt + n + 1)\n",
" 519 4375 2995000.0 684.6 0.2 left_scaling[:npt] = 1.0 / scale**2.0\n",
" 520 4375 1172000.0 267.9 0.1 left_scaling[npt] = scale**2.0\n",
" 521 4375 2013000.0 460.1 0.1 left_scaling[npt + 1:] = scale\n",
" 522 4375 2583000.0 590.4 0.2 right_scaling[:npt] = 1.0 / scale**2.0\n",
" 523 4375 1040000.0 237.7 0.1 right_scaling[npt] = scale**2.0\n",
" 524 4375 1981000.0 452.8 0.1 right_scaling[npt + 1:] = scale\n",
" 525 \n",
" 526 # Build the solution. After a discussion with Mike Saunders and Alexis\n",
" 527 # Montoison during their visit to the Hong Kong Polytechnic University\n",
" 528 # in 2024, we decided to use the eigendecomposition of the symmetric\n",
" 529 # matrix a. This is more stable than the previously employed LBL\n",
" 530 # decomposition, and allows us to directly detect ill-conditioning of\n",
" 531 # the system and to build the least-squares solution if necessary.\n",
" 532 # Numerical experiments have shown that this strategy improves the\n",
" 533 # performance of the solver.\n",
" 534 4375 4837000.0 1105.6 0.3 rhs_scaled = rhs * right_scaling[:, np.newaxis]\n",
" 535 4375 40424000.0 9239.8 2.8 if not (np.all(np.isfinite(a)) and np.all(np.isfinite(rhs_scaled))):\n",
" 536 raise np.linalg.LinAlgError(\n",
" 537 \"The interpolation system is ill-defined.\"\n",
" 538 )\n",
" 539 4375 1135764000.0 259603.2 79.6 eig_values, eig_vectors = eigh(a, check_finite=False)\n",
" 540 4375 6668000.0 1524.1 0.5 large_eig_values = np.abs(eig_values) > EPS\n",
" 541 4375 12130000.0 2772.6 0.9 eig_vectors = eig_vectors[:, large_eig_values]\n",
" 542 4375 4908000.0 1121.8 0.3 inv_eig_values = 1.0 / eig_values[large_eig_values]\n",
" 543 4375 18073000.0 4131.0 1.3 ill_conditioned = ~np.all(large_eig_values, 0)\n",
" 544 8750 5683000.0 649.5 0.4 left_scaled_solutions = eig_vectors @ (\n",
" 545 4375 11554000.0 2640.9 0.8 (eig_vectors.T @ rhs_scaled) * inv_eig_values[:, np.newaxis]\n",
" 546 )\n",
" 547 4375 5719000.0 1307.2 0.4 return (\n",
" 548 4375 3843000.0 878.4 0.3 left_scaled_solutions * right_scaling[:, np.newaxis],\n",
" 549 4375 426000.0 97.4 0.0 ill_conditioned,\n",
" 550 )\n",
"\n",
"Total time: 0.557575 s\n",
"File: /Users/andrew/Documents/Andy/programming/scipy/build-install/lib/python3.12/site-packages/scipy/_lib/cobyqa/models.py\n",
"Function: _get_model at line 552\n",
"\n",
"Line # Hits Time Per Hit % Time Line Contents\n",
"==============================================================\n",
" 552 @staticmethod\n",
" 553 def _get_model(interpolation, values):\n",
" 554 \"\"\"\n",
" 555 Solve the interpolation system.\n",
" 556 \n",
" 557 Parameters\n",
" 558 ----------\n",
" 559 interpolation : `cobyqa.models.Interpolation`\n",
" 560 Interpolation set.\n",
" 561 values : `numpy.ndarray`, shape (npt,)\n",
" 562 Values of the interpolated function at the interpolation points.\n",
" 563 \n",
" 564 Returns\n",
" 565 -------\n",
" 566 float\n",
" 567 Constant term of the quadratic model.\n",
" 568 `numpy.ndarray`, shape (n,)\n",
" 569 Gradient of the quadratic model at ``interpolation.x_base``.\n",
" 570 `numpy.ndarray`, shape (npt,)\n",
" 571 Implicit Hessian matrix of the quadratic model.\n",
" 572 \n",
" 573 Raises\n",
" 574 ------\n",
" 575 `numpy.linalg.LinAlgError`\n",
" 576 If the interpolation system is ill-defined.\n",
" 577 \"\"\"\n",
" 578 1577 1130000.0 716.6 0.2 assert values.shape == (interpolation.npt,), \\\n",
" 579 \"The shape of `values` is not valid.\"\n",
" 580 1577 501000.0 317.7 0.1 n, npt = interpolation.xpt.shape\n",
" 581 3154 519816000.0 164811.7 93.2 x, ill_conditioned = Quadratic.solve_systems(\n",
" 582 1577 129000.0 81.8 0.0 interpolation,\n",
" 583 3154 33272000.0 10549.1 6.0 np.block([[\n",
" 584 1577 126000.0 79.9 0.0 values,\n",
" 585 1577 512000.0 324.7 0.1 np.zeros(n + 1),\n",
" 586 1577 399000.0 253.0 0.1 ]]).T,\n",
" 587 )\n",
" 588 1577 1690000.0 1071.7 0.3 return x[npt, 0], x[npt + 1:, 0], x[:npt, 0], ill_conditioned\n",
"\n",
"Total time: 0.965335 s\n",
"File: /Users/andrew/Documents/Andy/programming/scipy/build-install/lib/python3.12/site-packages/scipy/_lib/cobyqa/models.py\n",
"Function: determinants at line 1301\n",
"\n",
"Line # Hits Time Per Hit % Time Line Contents\n",
"==============================================================\n",
" 1301 def determinants(self, x_new, k_new=None):\n",
" 1302 \"\"\"\n",
" 1303 Compute the normalized determinants of the new interpolation systems.\n",
" 1304 \n",
" 1305 Parameters\n",
" 1306 ----------\n",
" 1307 x_new : `numpy.ndarray`, shape (n,)\n",
" 1308 New interpolation point. Its value is interpreted as relative to\n",
" 1309 the origin, not the base point.\n",
" 1310 k_new : int, optional\n",
" 1311 Index of the updated interpolation point. If `k_new` is not\n",
" 1312 specified, all the possible determinants are computed.\n",
" 1313 \n",
" 1314 Returns\n",
" 1315 -------\n",
" 1316 {float, `numpy.ndarray`, shape (npt,)}\n",
" 1317 Determinant(s) of the new interpolation system.\n",
" 1318 \n",
" 1319 Raises\n",
" 1320 ------\n",
" 1321 `numpy.linalg.LinAlgError`\n",
" 1322 If the interpolation system is ill-defined.\n",
" 1323 \n",
" 1324 Notes\n",
" 1325 -----\n",
" 1326 The determinants are normalized by the determinant of the current\n",
" 1327 interpolation system. For stability reasons, the calculations are done\n",
" 1328 using the formula (2.12) in [1]_.\n",
" 1329 \n",
" 1330 References\n",
" 1331 ----------\n",
" 1332 .. [1] M. J. D. Powell. On updating the inverse of a KKT matrix.\n",
" 1333 Technical Report DAMTP 2004/NA01, Department of Applied Mathematics\n",
" 1334 and Theoretical Physics, University of Cambridge, Cambridge, UK,\n",
" 1335 2004.\n",
" 1336 \"\"\"\n",
" 1337 1399 232000.0 165.8 0.0 if self._debug:\n",
" 1338 assert x_new.shape == (self.n,), \\\n",
" 1339 \"The shape of `x_new` is not valid.\"\n",
" 1340 assert k_new is None or 0 <= k_new < self.npt, \\\n",
" 1341 \"The index `k_new` is not valid.\"\n",
" 1342 \n",
" 1343 # Compute the values independent of k_new.\n",
" 1344 1399 1194000.0 853.5 0.1 shift = x_new - self.interpolation.x_base\n",
" 1345 1399 2998000.0 2143.0 0.3 new_col = np.empty((self.npt + self.n + 1, 1))\n",
" 1346 1399 1971000.0 1408.9 0.2 new_col[: self.npt, 0] = (\n",
" 1347 1399 3534000.0 2526.1 0.4 0.5 * (self.interpolation.xpt.T @ shift) ** 2.0\n",
" 1348 )\n",
" 1349 1399 1461000.0 1044.3 0.2 new_col[self.npt, 0] = 1.0\n",
" 1350 1399 1666000.0 1190.9 0.2 new_col[self.npt + 1:, 0] = shift\n",
" 1351 1399 465930000.0 333045.0 48.3 inv_new_col = Quadratic.solve_systems(self.interpolation, new_col)[0]\n",
" 1352 1399 3092000.0 2210.2 0.3 beta = 0.5 * (shift @ shift) ** 2.0 - new_col[:, 0] @ inv_new_col[:, 0]\n",
" 1353 \n",
" 1354 # Compute the values that depend on k.\n",
" 1355 1399 226000.0 161.5 0.0 if k_new is None:\n",
" 1356 745 4750000.0 6375.8 0.5 coord_vec = np.eye(self.npt + self.n + 1, self.npt)\n",
" 1357 2235 256186000.0 114624.6 26.5 alpha = np.diag(Quadratic.solve_systems(\n",
" 1358 745 231000.0 310.1 0.0 self.interpolation,\n",
" 1359 745 73000.0 98.0 0.0 coord_vec,\n",
" 1360 745 104000.0 139.6 0.0 )[0])\n",
" 1361 745 1016000.0 1363.8 0.1 tau = inv_new_col[:self.npt, 0]\n",
" 1362 else:\n",
" 1363 654 3228000.0 4935.8 0.3 coord_vec = np.eye(self.npt + self.n + 1, 1, -k_new)\n",
" 1364 2616 213179000.0 81490.4 22.1 alpha = Quadratic.solve_systems(\n",
" 1365 654 191000.0 292.0 0.0 self.interpolation,\n",
" 1366 654 59000.0 90.2 0.0 coord_vec,\n",
" 1367 1308 162000.0 123.9 0.0 )[0][k_new, 0]\n",
" 1368 654 163000.0 249.2 0.0 tau = inv_new_col[k_new, 0]\n",
" 1369 1399 3689000.0 2636.9 0.4 return alpha * beta + tau**2.0\n",
"\n",
"Total time: 0.224351 s\n",
"File: /Users/andrew/Documents/Andy/programming/scipy/build-install/lib/python3.12/site-packages/scipy/_lib/cobyqa/subsolvers/geometry.py\n",
"Function: spider_geometry at line 106\n",
"\n",
"Line # Hits Time Per Hit % Time Line Contents\n",
"==============================================================\n",
" 106 def spider_geometry(const, grad, curv, xpt, xl, xu, delta, debug):\n",
" 107 r\"\"\"\n",
" 108 Maximize approximately the absolute value of a quadratic function subject\n",
" 109 to bound constraints in a trust region.\n",
" 110 \n",
" 111 This function solves approximately\n",
" 112 \n",
" 113 .. math::\n",
" 114 \n",
" 115 \\max_{s \\in \\mathbb{R}^n} \\quad \\bigg\\lvert c + g^{\\mathsf{T}} s +\n",
" 116 \\frac{1}{2} s^{\\mathsf{T}} H s \\bigg\\rvert \\quad \\text{s.t.} \\quad\n",
" 117 \\left\\{ \\begin{array}{l}\n",
" 118 l \\le s \\le u,\\\\\n",
" 119 \\lVert s \\rVert \\le \\Delta,\n",
" 120 \\end{array} \\right.\n",
" 121 \n",
" 122 by maximizing the objective function along given straight lines.\n",
" 123 \n",
" 124 Parameters\n",
" 125 ----------\n",
" 126 const : float\n",
" 127 Constant :math:`c` as shown above.\n",
" 128 grad : `numpy.ndarray`, shape (n,)\n",
" 129 Gradient :math:`g` as shown above.\n",
" 130 curv : callable\n",
" 131 Curvature of :math:`H` along any vector.\n",
" 132 \n",
" 133 ``curv(s) -> float``\n",
" 134 \n",
" 135 returns :math:`s^{\\mathsf{T}} H s`.\n",
" 136 xpt : `numpy.ndarray`, shape (n, npt)\n",
" 137 Points defining the straight lines. The straight lines considered are\n",
" 138 the ones passing through the origin and the points in `xpt`.\n",
" 139 xl : `numpy.ndarray`, shape (n,)\n",
" 140 Lower bounds :math:`l` as shown above.\n",
" 141 xu : `numpy.ndarray`, shape (n,)\n",
" 142 Upper bounds :math:`u` as shown above.\n",
" 143 delta : float\n",
" 144 Trust-region radius :math:`\\Delta` as shown above.\n",
" 145 debug : bool\n",
" 146 Whether to make debugging tests during the execution.\n",
" 147 \n",
" 148 Returns\n",
" 149 -------\n",
" 150 `numpy.ndarray`, shape (n,)\n",
" 151 Approximate solution :math:`s`.\n",
" 152 \n",
" 153 Notes\n",
" 154 -----\n",
" 155 This function is described as the second alternative in Section 6.5 of\n",
" 156 [1]_. It is assumed that the origin is feasible with respect to the bound\n",
" 157 constraints and that `delta` is finite and positive.\n",
" 158 \n",
" 159 References\n",
" 160 ----------\n",
" 161 .. [1] T. M. Ragonneau. *Model-Based Derivative-Free Optimization Methods\n",
" 162 and Software*. PhD thesis, Department of Applied Mathematics, The Hong\n",
" 163 Kong Polytechnic University, Hong Kong, China, 2022. URL:\n",
" 164 https://theses.lib.polyu.edu.hk/handle/200/12294.\n",
" 165 \"\"\"\n",
" 166 327 51000.0 156.0 0.0 if debug:\n",
" 167 assert isinstance(const, float)\n",
" 168 assert isinstance(grad, np.ndarray) and grad.ndim == 1\n",
" 169 assert inspect.signature(curv).bind(grad)\n",
" 170 assert (\n",
" 171 isinstance(xpt, np.ndarray)\n",
" 172 and xpt.ndim == 2\n",
" 173 and xpt.shape[0] == grad.size\n",
" 174 )\n",
" 175 assert isinstance(xl, np.ndarray) and xl.shape == grad.shape\n",
" 176 assert isinstance(xu, np.ndarray) and xu.shape == grad.shape\n",
" 177 assert isinstance(delta, float)\n",
" 178 assert isinstance(debug, bool)\n",
" 179 tol = get_arrays_tol(xl, xu)\n",
" 180 assert np.all(xl <= tol)\n",
" 181 assert np.all(xu >= -tol)\n",
" 182 assert np.isfinite(delta) and delta > 0.0\n",
" 183 327 311000.0 951.1 0.1 xl = np.minimum(xl, 0.0)\n",
" 184 327 219000.0 669.7 0.1 xu = np.maximum(xu, 0.0)\n",
" 185 \n",
" 186 # Iterate through the straight lines.\n",
" 187 327 1011000.0 3091.7 0.5 step = np.zeros_like(grad)\n",
" 188 327 41000.0 125.4 0.0 q_val = const\n",
" 189 327 7335000.0 22431.2 3.3 s_norm = np.linalg.norm(xpt, axis=0)\n",
" 190 \n",
" 191 327 1393000.0 4259.9 0.6 i_xl_pos = ((xl > -np.inf) & (xpt.T > -TINY * xl))\n",
" 192 327 1152000.0 3522.9 0.5 i_xl_neg = ((xl > -np.inf) & (xpt.T < TINY * xl))\n",
" 193 327 1053000.0 3220.2 0.5 i_xu_pos = ((xu < np.inf) & (xpt.T > TINY * xu))\n",
" 194 327 1058000.0 3235.5 0.5 i_xu_neg = ((xu < np.inf) & (xpt.T < -TINY * xu))\n",
" 195 6867 1198000.0 174.5 0.5 for k in range(xpt.shape[1]):\n",
" 196 # Set alpha_tr to the step size for the trust-region constraint.\n",
" 197 # s_norm = np.linalg.norm(xpt[:, k])\n",
" 198 6540 1745000.0 266.8 0.8 if s_norm[k] > TINY * delta:\n",
" 199 6540 2186000.0 334.3 1.0 alpha_tr = max(delta / s_norm[k], 0.0)\n",
" 200 else:\n",
" 201 # The current straight line is basically zero.\n",
" 202 continue\n",
" 203 \n",
" 204 # Set alpha_xl to the step size for the lower-bound constraint and\n",
" 205 # alpha_xu to the step size for the upper-bound constraint.\n",
" 206 # i_xl_pos = (xl > -np.inf) & (xpt[:, k] > -TINY * xl)\n",
" 207 # i_xl_neg = (xl > -np.inf) & (xpt[:, k] < TINY * xl)\n",
" 208 # i_xu_pos = (xu < np.inf) & (xpt[:, k] > TINY * xu)\n",
" 209 # i_xu_neg = (xu < np.inf) & (xpt[:, k] < -TINY * xu)\n",
" 210 6540 33981000.0 5195.9 15.1 alpha_xl_pos = np.max(xl[i_xl_pos[k]] / xpt[i_xl_pos[k], k], initial=-np.inf)\n",
" 211 6540 32645000.0 4991.6 14.6 alpha_xl_neg = np.min(xl[i_xl_neg[k]] / xpt[i_xl_neg[k], k], initial=np.inf)\n",
" 212 6540 32466000.0 4964.2 14.5 alpha_xu_neg = np.max(xu[i_xu_neg[k]] / xpt[i_xu_neg[k], k], initial=-np.inf)\n",
" 213 6540 32154000.0 4916.5 14.3 alpha_xu_pos = np.min(xu[i_xu_pos[k]] / xpt[i_xu_pos[k], k], initial=np.inf)\n",
" 214 6540 2480000.0 379.2 1.1 alpha_bd_pos = max(min(alpha_xu_pos, alpha_xl_neg), 0.0)\n",
" 215 6540 2123000.0 324.6 0.9 alpha_bd_neg = min(max(alpha_xl_pos, alpha_xu_neg), 0.0)\n",
" 216 \n",
" 217 # Set alpha_quad_pos and alpha_quad_neg to the step size to the extrema\n",
" 218 # of the quadratic function along the positive and negative directions.\n",
" 219 6540 6937000.0 1060.7 3.1 grad_step = grad @ xpt[:, k]\n",
" 220 6540 30971000.0 4735.6 13.8 curv_step = curv(xpt[:, k])\n",
" 221 if (\n",
" 222 6540 1112000.0 170.0 0.5 grad_step >= 0.0\n",
" 223 2682 647000.0 241.2 0.3 and curv_step < -TINY * grad_step\n",
" 224 3943 583000.0 147.9 0.3 or grad_step <= 0.0\n",
" 225 3858 847000.0 219.5 0.4 and curv_step > -TINY * grad_step\n",
" 226 ):\n",
" 227 6451 2133000.0 330.6 1.0 alpha_quad_pos = max(-grad_step / curv_step, 0.0)\n",
" 228 else:\n",
" 229 89 14000.0 157.3 0.0 alpha_quad_pos = np.inf\n",
" 230 if (\n",
" 231 6540 959000.0 146.6 0.4 grad_step >= 0.0\n",
" 232 2682 527000.0 196.5 0.2 and curv_step > TINY * grad_step\n",
" 233 6455 943000.0 146.1 0.4 or grad_step <= 0.0\n",
" 234 3858 747000.0 193.6 0.3 and curv_step < TINY * grad_step\n",
" 235 ):\n",
" 236 89 30000.0 337.1 0.0 alpha_quad_neg = min(-grad_step / curv_step, 0.0)\n",
" 237 else:\n",
" 238 6451 985000.0 152.7 0.4 alpha_quad_neg = -np.inf\n",
" 239 \n",
" 240 # Select the step that provides the largest value of the objective\n",
" 241 # function if it improves the current best. The best positive step is\n",
" 242 # either the one that reaches the constraints or the one that reaches\n",
" 243 # the extremum of the objective function along the current direction\n",
" 244 # (only possible if the resulting step is feasible). We test both, and\n",
" 245 # we perform similar calculations along the negative step.\n",
" 246 # N.B.: we select the largest possible step among all the ones that\n",
" 247 # maximize the objective function. This is to avoid returning the zero\n",
" 248 # step in some extreme cases.\n",
" 249 6540 1468000.0 224.5 0.7 alpha_pos = min(alpha_tr, alpha_bd_pos)\n",
" 250 6540 1523000.0 232.9 0.7 alpha_neg = max(-alpha_tr, alpha_bd_neg)\n",
" 251 6540 553000.0 84.6 0.2 q_val_pos = (\n",
" 252 6540 2271000.0 347.2 1.0 const + alpha_pos * grad_step + 0.5 * alpha_pos**2.0 * curv_step\n",
" 253 )\n",
" 254 6540 638000.0 97.6 0.3 q_val_neg = (\n",
" 255 6540 2011000.0 307.5 0.9 const + alpha_neg * grad_step + 0.5 * alpha_neg**2.0 * curv_step\n",
" 256 )\n",
" 257 6540 973000.0 148.8 0.4 if alpha_quad_pos < alpha_pos:\n",
" 258 2214 220000.0 99.4 0.1 q_val_quad_pos = (\n",
" 259 6642 920000.0 138.5 0.4 const\n",
" 260 2214 270000.0 122.0 0.1 + alpha_quad_pos * grad_step\n",
" 261 2214 535000.0 241.6 0.2 + 0.5 * alpha_quad_pos**2.0 * curv_step\n",
" 262 )\n",
" 263 2214 503000.0 227.2 0.2 if abs(q_val_quad_pos) > abs(q_val_pos):\n",
" 264 1919 184000.0 95.9 0.1 alpha_pos = alpha_quad_pos\n",
" 265 1919 204000.0 106.3 0.1 q_val_pos = q_val_quad_pos\n",
" 266 6540 1008000.0 154.1 0.4 if alpha_quad_neg > alpha_neg:\n",
" 267 61 2000.0 32.8 0.0 q_val_quad_neg = (\n",
" 268 183 20000.0 109.3 0.0 const\n",
" 269 61 6000.0 98.4 0.0 + alpha_quad_neg * grad_step\n",
" 270 61 14000.0 229.5 0.0 + 0.5 * alpha_quad_neg**2.0 * curv_step\n",
" 271 )\n",
" 272 61 14000.0 229.5 0.0 if abs(q_val_quad_neg) > abs(q_val_neg):\n",
" 273 17 2000.0 117.6 0.0 alpha_neg = alpha_quad_neg\n",
" 274 17 0.0 0.0 0.0 q_val_neg = q_val_quad_neg\n",
" 275 6540 1523000.0 232.9 0.7 if abs(q_val_pos) >= abs(q_val_neg) and abs(q_val_pos) > abs(q_val):\n",
" 276 54 194000.0 3592.6 0.1 step = np.clip(alpha_pos * xpt[:, k], xl, xu)\n",
" 277 54 8000.0 148.1 0.0 q_val = q_val_pos\n",
" 278 6486 2149000.0 331.3 1.0 elif abs(q_val_neg) > abs(q_val_pos) and abs(q_val_neg) > abs(q_val):\n",
" 279 1035 3860000.0 3729.5 1.7 step = np.clip(alpha_neg * xpt[:, k], xl, xu)\n",
" 280 1035 138000.0 133.3 0.1 q_val = q_val_neg\n",
" 281 \n",
" 282 327 42000.0 128.4 0.0 if debug:\n",
" 283 assert np.all(xl <= step)\n",
" 284 assert np.all(step <= xu)\n",
" 285 assert np.linalg.norm(step) < 1.1 * delta\n",
" 286 327 1061000.0 3244.6 0.5 return step"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%lprun -f Quadratic.build_system -f Quadratic._get_model -f Quadratic.solve_systems -f TrustRegion.get_geometry_step -f spider_geometry -f Models.determinants minimize(rosen, x0, method='cobyqa')"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "3cab8725-8ed8-4326-a0fb-a747fe15fba8",
"metadata": {},
"outputs": [],
"source": [
"a = np.arange(200.).reshape(10, 20)\n",
"b = np.arange(10.)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "cc217e20-d47f-4567-b79f-277e7624a9d1",
"metadata": {},
"outputs": [
{
"ename": "ValueError",
"evalue": "operands could not be broadcast together with shapes (10,20) (10,) ",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)",
"Cell \u001b[0;32mIn[8], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m ((b \u001b[38;5;241m>\u001b[39m \u001b[38;5;241m-\u001b[39mnp\u001b[38;5;241m.\u001b[39minf) \u001b[38;5;241m&\u001b[39m (\u001b[43ma\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m<\u001b[39;49m\u001b[43m \u001b[49m\u001b[43mb\u001b[49m))\u001b[38;5;241m.\u001b[39mshape\n",
"\u001b[0;31mValueError\u001b[0m: operands could not be broadcast together with shapes (10,20) (10,) "
]
}
],
"source": [
"((b > -np.inf) & (a < b)).shape"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "420e356f-2ead-49ec-8edf-bd3b6720baf4",
"metadata": {},
"outputs": [],
"source": [
"np.max("
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "36585e02-6c16-4cbd-a09b-097b4095e741",
"metadata": {},
"outputs": [],
"source": [
"np.broadcast_to("
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"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.12.3"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment