{ "cells": [ { "cell_type": "markdown", "metadata": { "id": "48HONqnFpIvI" }, "source": [ "# Problem 7\n", "Let $A$ be the $20000 \\times 20000$ matrix whose entries are 0 everywhere except for the primes 2,3,5,7,...,224737 along the main diagonal and the number 1 in all the positions $ a_{ij} $ with $ |i-j| = 1,2,4,8,\\cdots,16384$. What is the (1,1) entry of $A^{-1}$" ] }, { "cell_type": "markdown", "metadata": { "id": "JiZSIfm1qR-F" }, "source": [ "### Some introductory code " ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 34 }, "executionInfo": { "elapsed": 1558, "status": "ok", "timestamp": 1603893832461, "user": { "displayName": "Husnain Raza", "photoUrl": "", "userId": "16146658556003895350" }, "user_tz": 300 }, "id": "wvpmKjZMqRiB", "outputId": "9ad6fbfd-f2b3-4ad4-eca6-97830b9d01b0" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "554466\n" ] } ], "source": [ "def is_prime(N):\n", " if N < 2: return False\n", " if N == 2: return True\n", " for i in range(2, int(N**0.5)+1):\n", " if N%i == 0: return False\n", " return True\n", "def is_prime(N):\n", " return N == 2 or any([N%i == 0 for i in range()])\n", "entries = []\n", "primes = [i for i in range(2,250000) if is_prime(i)][:20000]\n", "entries.extend([(i,i,primes[i]) for i in range(len(primes))])\n", "for diff in [2**j for j in range(15)]:\n", " for i in range(20000):\n", " j1, j2 = i+diff, i-diff\n", " if 0<=j1<20000: entries.append((i,j1,1))\n", " if 0<=j2<20000: entries.append((i,j2,1))\n", "print(len(entries))" ] }, { "cell_type": "markdown", "metadata": { "id": "vsdtHzxT3WQc" }, "source": [ "### Attempt 1: Using standard libraries\n", "\n", "We avoid calculating $A^{-1}$ explicitly. Instead, note that if we have $\\hat{v}$ be the first column of $ A^{-1} $, then we get that $ A \\hat{v} = (1,0,\\cdots, 0)^{T}$. Let's solve this system. " ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 229 }, "executionInfo": { "elapsed": 789, "status": "error", "timestamp": 1603893823785, "user": { "displayName": "Husnain Raza", "photoUrl": "", "userId": "16146658556003895350" }, "user_tz": 300 }, "id": "3QZE1MQao2Zz", "outputId": "91fb2579-b63b-4656-d445-18ff6f771c42" }, "outputs": [ { "ename": "NameError", "evalue": "ignored", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mscipy\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msparse\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinalg\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mspsolve\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 4\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 5\u001b[0;31m \u001b[0mrow\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mentries\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 6\u001b[0m \u001b[0mcol\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mentries\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 7\u001b[0m \u001b[0mdata\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mentries\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mNameError\u001b[0m: name 'entries' is not defined" ] } ], "source": [ "import numpy as np\n", "from scipy.sparse import coo_matrix\n", "from scipy.sparse.linalg import spsolve\n", "\n", "row = np.array([i[0] for i in entries])\n", "col = np.array([i[1] for i in entries])\n", "data = np.array([i[2] for i in entries])\n", "\n", "A = coo_matrix((data, (row, col))).tocsc()\n", "b = np.zeros(A.shape[0])\n", "b[0] = 1\n", "v = spsolve(A,b)\n", "# Timeout" ] }, { "cell_type": "markdown", "metadata": { "id": "NsvKncT7ibZS" }, "source": [ "The above code times out - so instead, let's try to use an iterative method to solve this linear system." ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 1000 }, "executionInfo": { "elapsed": 13545, "status": "ok", "timestamp": 1603894098393, "user": { "displayName": "Husnain Raza", "photoUrl": "", "userId": "16146658556003895350" }, "user_tz": 300 }, "id": "NjRjx9eWr4C7", "outputId": "df01ffbf-12b1-42b3-8826-607105d1f96a" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Iteration 10: current val = 0.5160474805069852\n", "Iteration 20: current val = 0.5424188205227113\n", "Iteration 30: current val = 0.5649984704641765\n", "Iteration 40: current val = 0.5809367071267324\n", "Iteration 50: current val = 0.6096730237142656\n", "Iteration 60: current val = 0.6190393486130085\n", "Iteration 70: current val = 0.6263245512238234\n", "Iteration 80: current val = 0.6508621776392641\n", "Iteration 90: current val = 0.6577928020000656\n", "Iteration 100: current val = 0.6608853920165397\n", "Iteration 110: current val = 0.6646370853741161\n", "Iteration 120: current val = 0.6769767188475972\n", "Iteration 130: current val = 0.686598661153931\n", "Iteration 140: current val = 0.6950283972972651\n", "Iteration 150: current val = 0.6985776576929462\n", "Iteration 160: current val = 0.700549226696957\n", "Iteration 170: current val = 0.7016993660417654\n", "Iteration 180: current val = 0.7037530670049756\n", "Iteration 190: current val = 0.70578202448334\n", "Iteration 200: current val = 0.7080325880243894\n", "Iteration 210: current val = 0.7107390862052434\n", "Iteration 220: current val = 0.7131266080608784\n", "Iteration 230: current val = 0.7154864307143483\n", "Iteration 240: current val = 0.7169760206590395\n", "Iteration 250: current val = 0.7177222840589774\n", "Iteration 260: current val = 0.7182998666211887\n", "Iteration 270: current val = 0.7189625336279692\n", "Iteration 280: current val = 0.7195270338881113\n", "Iteration 290: current val = 0.7198970070328818\n", "Iteration 300: current val = 0.7202715787882576\n", "Iteration 310: current val = 0.720548747555294\n", "Iteration 320: current val = 0.7208623428688858\n", "Iteration 330: current val = 0.7211885297932319\n", "Iteration 340: current val = 0.7216667311611941\n", "Iteration 350: current val = 0.7219940973832849\n", "Iteration 360: current val = 0.7224389483334973\n", "Iteration 370: current val = 0.7228169933348944\n", "Iteration 380: current val = 0.7231983252626363\n", "Iteration 390: current val = 0.7234404924750876\n", "Iteration 400: current val = 0.7236478043297795\n", "Iteration 410: current val = 0.7238808716336398\n", "Iteration 420: current val = 0.7241166500666772\n", "Iteration 430: current val = 0.7243659045725263\n", "Iteration 440: current val = 0.7244940472057783\n", "Iteration 450: current val = 0.724633608525676\n", "Iteration 460: current val = 0.7247020541105452\n", "Iteration 470: current val = 0.7247680447799154\n", "Iteration 480: current val = 0.7248237539591655\n", "Iteration 490: current val = 0.7248677969693484\n", "Iteration 500: current val = 0.7248916848307266\n", "Iteration 510: current val = 0.7249208479337073\n", "Iteration 520: current val = 0.7249533593512273\n", "Iteration 530: current val = 0.7249736254056114\n", "Iteration 540: current val = 0.7249887123438771\n", "Iteration 550: current val = 0.7250067816490801\n", "Iteration 560: current val = 0.7250207131481929\n", "Iteration 570: current val = 0.7250304954206683\n", "Iteration 580: current val = 0.7250427286641773\n", "Iteration 590: current val = 0.72505125564296\n", "Iteration 600: current val = 0.7250561238493789\n", "Iteration 610: current val = 0.7250605939852699\n", "Iteration 620: current val = 0.7250640020707945\n", "Iteration 630: current val = 0.7250660938618454\n", "Iteration 640: current val = 0.7250677953689088\n", "Iteration 650: current val = 0.7250693490905185\n", "Iteration 660: current val = 0.725070446117053\n", "Iteration 670: current val = 0.7250713225273709\n", "Iteration 680: current val = 0.7250726985038631\n", "Iteration 690: current val = 0.7250739362881828\n", "Iteration 700: current val = 0.7250748779230066\n", "Iteration 710: current val = 0.7250755450783181\n", "Iteration 720: current val = 0.7250760283885833\n", "Iteration 730: current val = 0.7250765020947028\n", "Iteration 740: current val = 0.7250770006680382\n", "Iteration 750: current val = 0.7250773496606184\n", "Iteration 760: current val = 0.7250776006228311\n", "Iteration 770: current val = 0.7250777636744249\n", "Iteration 780: current val = 0.7250778756720382\n", "Iteration 790: current val = 0.7250779688442739\n", "Iteration 800: current val = 0.7250780439687791\n", "Iteration 810: current val = 0.7250781067493376\n", "Iteration 820: current val = 0.7250781516265811\n", "Iteration 830: current val = 0.7250781883759363\n", "Iteration 840: current val = 0.7250782241898078\n", "Iteration 850: current val = 0.7250782614326434\n", "Iteration 860: current val = 0.7250782837700643\n", "Iteration 870: current val = 0.7250782994297919\n", "Iteration 880: current val = 0.7250783112476437\n", "Iteration 890: current val = 0.7250783221719436\n", "Iteration 900: current val = 0.725078329800566\n", "Iteration 910: current val = 0.7250783351385289\n", "Iteration 920: current val = 0.7250783385028293\n", "Iteration 930: current val = 0.7250783410043907\n", "Iteration 940: current val = 0.725078342638123\n", "Iteration 950: current val = 0.7250783438331366\n", "Iteration 960: current val = 0.7250783446050801\n", "Iteration 970: current val = 0.7250783450804342\n", "Iteration 980: current val = 0.7250783453953678\n", "Iteration 990: current val = 0.725078345614237\n", "Iteration 1000: current val = 0.7250783457924285\n", "Iteration 1010: current val = 0.7250783459042496\n", "Iteration 1020: current val = 0.7250783459839163\n", "Iteration 1030: current val = 0.7250783460699098\n", "Iteration 1040: current val = 0.7250783461385248\n", "Iteration 1050: current val = 0.7250783461804372\n", "Iteration 1060: current val = 0.7250783462047834\n", "Iteration 1070: current val = 0.7250783462184939\n", "Iteration 1080: current val = 0.7250783462318636\n", "Iteration 1090: current val = 0.7250783462415041\n", "Iteration 1100: current val = 0.7250783462487344\n", "Iteration 1110: current val = 0.7250783462541678\n", "Iteration 1120: current val = 0.7250783462581615\n", "Iteration 1130: current val = 0.72507834626106\n", "Iteration 1140: current val = 0.725078346263089\n", "Iteration 1150: current val = 0.72507834626454\n", "Iteration 1160: current val = 0.7250783462657099\n", "Iteration 1170: current val = 0.7250783462665752\n", "Iteration 1180: current val = 0.7250783462671782\n", "Iteration 1190: current val = 0.7250783462676181\n", "Iteration 1200: current val = 0.7250783462679155\n", "Iteration 1210: current val = 0.7250783462681067\n", "Iteration 1220: current val = 0.7250783462682162\n", "Iteration 1230: current val = 0.7250783462682798\n", "Iteration 1240: current val = 0.7250783462683179\n", "Iteration 1250: current val = 0.725078346268345\n", "Iteration 1260: current val = 0.7250783462683673\n", "Iteration 1270: current val = 0.7250783462683824\n", "Iteration 1280: current val = 0.7250783462683897\n", "Iteration 1290: current val = 0.7250783462683944\n", "Iteration 1300: current val = 0.7250783462683967\n", "Iteration 1310: current val = 0.7250783462683981\n", "Iteration 1320: current val = 0.7250783462683992\n", "Iteration 1330: current val = 0.7250783462684001\n", "Iteration 1340: current val = 0.7250783462684001\n", "Iteration 1350: current val = 0.7250783462684001\n", "Iteration 1360: current val = 0.7250783462684001\n", "Iteration 1370: current val = 0.7250783462684001\n", "Iteration 1380: current val = 0.7250783462684001\n", "Iteration 1390: current val = 0.7250783462684001\n", "Iteration 1400: current val = 0.7250783462684001\n", "Iteration 1410: current val = 0.7250783462684001\n", "Iteration 1420: current val = 0.7250783462684001\n", "Iteration 1430: current val = 0.7250783462684001\n", "Iteration 1440: current val = 0.7250783462684001\n", "Iteration 1450: current val = 0.7250783462684001\n", "Iteration 1460: current val = 0.7250783462684001\n", "Iteration 1470: current val = 0.7250783462684001\n", "Iteration 1480: current val = 0.7250783462684001\n", "Iteration 1490: current val = 0.7250783462684001\n", "Iteration 1500: current val = 0.7250783462684001\n", "Iteration 1510: current val = 0.7250783462684001\n", "Iteration 1520: current val = 0.7250783462684001\n", "Iteration 1530: current val = 0.7250783462684001\n", "Iteration 1540: current val = 0.7250783462684001\n", "Iteration 1550: current val = 0.7250783462684001\n", "Iteration 1560: current val = 0.7250783462684001\n", "Iteration 1570: current val = 0.7250783462684001\n", "Iteration 1580: current val = 0.7250783462684001\n", "Iteration 1590: current val = 0.7250783462684001\n", "Iteration 1600: current val = 0.7250783462684001\n", "Iteration 1610: current val = 0.7250783462684001\n", "Iteration 1620: current val = 0.7250783462684001\n", "Iteration 1630: current val = 0.7250783462684001\n", "Iteration 1640: current val = 0.7250783462684001\n", "Iteration 1650: current val = 0.7250783462684001\n", "Iteration 1660: current val = 0.7250783462684001\n", "Iteration 1670: current val = 0.7250783462684001\n", "Iteration 1680: current val = 0.7250783462684001\n", "Iteration 1690: current val = 0.7250783462684001\n", "Iteration 1700: current val = 0.7250783462684001\n", "Iteration 1710: current val = 0.7250783462684001\n", "Iteration 1720: current val = 0.7250783462684001\n", "Iteration 1730: current val = 0.7250783462684001\n", "Iteration 1740: current val = 0.7250783462684001\n", "Iteration 1750: current val = 0.7250783462684001\n", "Iteration 1760: current val = 0.7250783462684001\n", "Iteration 1770: current val = 0.7250783462684001\n", "Iteration 1780: current val = 0.7250783462684001\n", "Iteration 1790: current val = 0.7250783462684001\n", "Iteration 1800: current val = 0.7250783462684001\n", "Iteration 1810: current val = 0.7250783462684001\n", "Iteration 1820: current val = 0.7250783462684001\n", "Iteration 1830: current val = 0.7250783462684001\n", "Iteration 1840: current val = 0.7250783462684001\n", "Iteration 1850: current val = 0.7250783462684001\n", "Iteration 1860: current val = 0.7250783462684001\n", "Iteration 1870: current val = 0.7250783462684001\n", "Iteration 1880: current val = 0.7250783462684001\n", "Iteration 1890: current val = 0.7250783462684001\n", "Iteration 1900: current val = 0.7250783462684001\n", "Iteration 1910: current val = 0.7250783462684001\n", "Iteration 1920: current val = 0.7250783462684001\n", "Iteration 1930: current val = 0.7250783462684001\n", "Iteration 1940: current val = 0.7250783462684001\n", "Iteration 1950: current val = 0.7250783462684001\n", "Iteration 1960: current val = 0.7250783462684001\n", "Iteration 1970: current val = 0.7250783462684001\n", "Iteration 1980: current val = 0.7250783462684001\n", "Iteration 1990: current val = 0.7250783462684001\n", "Iteration 2000: current val = 0.7250783462684001\n", "Iteration 2010: current val = 0.7250783462684001\n", "Iteration 2020: current val = 0.7250783462684001\n", "Iteration 2030: current val = 0.7250783462684001\n", "Iteration 2040: current val = 0.7250783462684001\n", "Iteration 2050: current val = 0.7250783462684001\n", "Iteration 2060: current val = 0.7250783462684001\n", "Iteration 2070: current val = 0.7250783462684001\n", "Iteration 2080: current val = 0.7250783462684001\n", "Iteration 2090: current val = 0.7250783462684001\n", "Iteration 2100: current val = 0.7250783462684001\n", "Iteration 2110: current val = 0.7250783462684001\n", "Iteration 2120: current val = 0.7250783462684001\n", "Iteration 2130: current val = 0.7250783462684001\n", "Iteration 2140: current val = 0.7250783462684001\n", "Iteration 2150: current val = 0.7250783462684001\n", "Iteration 2160: current val = 0.7250783462684001\n", "Iteration 2170: current val = 0.7250783462684001\n", "Iteration 2180: current val = 0.7250783462684001\n", "Iteration 2190: current val = 0.7250783462684001\n", "Iteration 2200: current val = 0.7250783462684001\n", "Iteration 2210: current val = 0.7250783462684001\n", "Iteration 2220: current val = 0.7250783462684001\n", "Iteration 2230: current val = 0.7250783462684001\n", "Iteration 2240: current val = 0.7250783462684001\n", "Iteration 2250: current val = 0.7250783462684001\n", "Iteration 2260: current val = 0.7250783462684001\n", "Iteration 2270: current val = 0.7250783462684001\n", "Iteration 2280: current val = 0.7250783462684001\n", "Iteration 2290: current val = 0.7250783462684001\n", "0.7250783462684001\n" ] } ], "source": [ "import numpy as np\n", "from scipy.sparse import coo_matrix\n", "from scipy.sparse.linalg import bicg\n", "row = np.array([i[0] for i in entries])\n", "col = np.array([i[1] for i in entries])\n", "data = np.array([i[2] for i in entries])\n", "\n", "A = coo_matrix((data, (row, col))).tocsc()\n", "b = np.zeros(A.shape[0])\n", "b[0] = 1\n", "\n", "iterations = 0\n", "def callback(xk):\n", " global iterations\n", " iterations += 1\n", " if iterations % 10 == 0: print(\"Iteration {}: current val = {}\".format(iterations, xk[0]))\n", "x, exit_code = bicg(A,b,tol=1e-50,callback=callback)\n", "print(x[0])" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "QWXAQN7K6gH7" }, "outputs": [], "source": [] } ], "metadata": { "colab": { "authorship_tag": "ABX9TyPEr5R+vBWEUgLKt8e+FhsB", "name": "Problem 7", "provenance": [] }, "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.10.6" } }, "nbformat": 4, "nbformat_minor": 1 }