Created
May 16, 2024 22:33
-
-
Save andyfaff/4cfdc554f9f520fb1a6a6f983025996f to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| { | |
| "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