{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# [NTDS'18] tutorial 5: sparse matrices in scipy\n", "[ntds'18]: https://github.com/mdeff/ntds_2018\n", "\n", "[Eda Bayram](http://lts4.epfl.ch/bayram), [EPFL LTS4](http://lts4.epfl.ch)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Ojective\n", "\n", "This is a short tutorial on the `scipy.sparse` module. We will talk about:\n", "\n", "1. What is sparsity?\n", "2. Sparse matrix storage schemes\n", "3. Linear operations on sparse matrices" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy import sparse\n", "import scipy.sparse.linalg\n", "from scipy import linalg\n", "import pandas as pd" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. Sparsity" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Why do we need data structures for sparse matrices?\n", "\n", "* Less memory usage\n", "* More efficiency computations\n", "\n", "Most real-world graphs / networks are sparse!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us create a random sparse matrix and analyze the sparsity." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Number of non-zeros: 625, density: 0.01\n" ] } ], "source": [ "N = 250 \n", "dummy = sparse.random(N, N, density=0.01)\n", "density = dummy.nnz / N**2\n", "print('Number of non-zeros: {}, density: {}'.format(dummy.nnz, density))" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQsAAAD/CAYAAADmIfPpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJztnX2MX9V5579PnDHyEJt4xoa6zBBParMSpOVlLLuyUcSKDSloFTaSiUilLjGt3D9AcmIi4Wz+AO1qtezKGLlLxdYRpkmbhhInFQglqqmVNsKuzXpYQyBobIexPGMsbGZcD/FM8GDO/vG7d3J9576ce16fc+d8pNHM/H73nvuc557znOc8542EEIhEIpE6PuFbgEgkEgbRWEQiESmisYhEIlJEYxGJRKSIxiISiUgRjUUkEpHCu7Egoj8iomEiOk5E23zLUwQRnSCiXxDRESI6nHzWQ0QvE9Gx5PdST7LtJqIzRPRm5rNC2ajDXyS6foOIbmUi72NEdCrR7xEiujvz3bcSeYeJ6IuOZe0nop8R0dtE9BYRbUk+Z6ffClnN6VYI4e0HwAIAvwLwWQALAbwO4AafMpXIeQLAstxn/wvAtuTvbQD+pyfZPg/gVgBv1skG4G4APwVAAP4QwCEm8j4G4JsF196QlIkrAAwkZWWBQ1lXALg1+XsxgKOJTOz0WyGrMd369izWAjguhHhHCHERwHMA7vEskyz3APhu8vd3AfwnH0IIIX4OYCL3cZls9wD4nuhwEMCniWiFG0k7lMhbxj0AnhNCfCiEGAFwHJ0y4wQhxGkhxGvJ3x8AeBvAtWCo3wpZy2isW9/G4loAo5n/x1CdQV8IAHuJaIiINiefXSOEOA10XhSAq71JN5cy2Tjr+6HEdd+d6dKxkZeIVgK4BcAhMNdvTlbAkG59Gwsq+Izj/PMNQohbAdwF4EEi+rxvgRThqu+nAfwegJsBnAbwRPI5C3mJ6FMAfgTg60KIyapLCz5zKm+BrMZ069tYjAHoz/zfB+BdT7KUIoR4N/l9BsA/oOOuvZe6mMnvM/4knEOZbCz1LYR4TwhxSQjxMYDv4LfusHd5iagLncr3fSHEj5OPWeq3SFaTuvVtLP4vgNVENEBECwHcB+BFzzJdBhFdSUSL078B3AngTXTkvD+57H4AL/iRsJAy2V4E8J+TqP0fAjifutM+yfXrv4yOfoGOvPcR0RVENABgNYBXHcpFAJ4B8LYQYkfmK3b6LZPVqG5dRWsrorh3oxO5/RWAb/uWp0C+z6ITNX4dwFupjAB6AewDcCz53eNJvh+g417OoNNa/GmZbOi4nn+Z6PoXANYwkfdvEnneSArxisz1307kHQZwl2NZb0PHNX8DwJHk526O+q2Q1ZhuKbkpEolEKvHdDYlEIoEQjUUkEpEiGotIJCJFNBaRSESKaCwikYgU1oxF09WkmWnU7AlJViAseUOSFQhLXl1ZrRgLIlqAznjzXeisbvsqEd1Qc1swSkdYsgJhyRuSrEBY8vIzFgh7NWkkEinAyqQsItoI4I+EEH+W/P8nANYJIR7KXLMZiaW78sorB3t7e7F8+XLjspjivcnf4MwHH+LqxVfgEx9+wFrWLBc/+hin3v83XLvs01j4Sf4hqrNnz3rR7cWPPsa5qYtY2r1QWk+pbq9eugS//vCjRvf6YGho6NdCiMWq93/SpDAZale0CSF2AdgFAGvWrBGHDx+2JIoZRiemsGdoDBsH+9Df0+1bHGmefPkodu47hk13rMY3vnB97fW288lVj031lL3nDwZ6cGhkotG9PiCiYZ37bRkL76sFTdPf011ZELhWgo2DfZf9rmPP0Bh27jsGAMoFv0oXJtK3QVM9Za/dsKoX+4+PN7o3RGwZi9nVpABOobOa9I8tPYsFXCtBnZHLo1Jp8lTpQiZ9H4a3qZ7y96wd6LUhFiusdLCEEB8BeAjAP6KzvdfzQoi3bDyriNGJKTz58lGMTky5eiQ2DvZhyx2rtVsXH7JnyVaAKjmq5KzSRZp+lRFIjc2eoTHFXERsYMuzgBDiJwB+Yiv9MkYnpvDNH76OQyOdbR5dtfIqLVMRXDyUMjnSVn9yegbPHjgx53tAXxcmvBuTuPJ0uHZlU6wZC1/sGRrDoZEJrBvoYVPYmuC6opQV0DI5UiOyaf1KI55UEaYMrylcGXAuDUUZrTMW2ULO0Tpzo6yA5itsalQ2rOr0zTnp13aL7MqAc/Oo8rTOWHBplVQLsOvWRbaA+mj1ZHVoWzZXZcrWc1I9YkHXQp10WmUsOPUtVQuw69ZFtoD6aPVkdci9RfZNqscF3Vct00mnVcaCU98yOwb/5MtHLxuLrzJkploX04bTh8cmawS4eJMmMfn+Uv1t3X7+fZ10+M5NVcDU8KWJ56QFeP/xcezcdwxP7D1qbTiwaBhTd/jR9xCuK/L55JJvk8PHs8b00sxFnXRa5Vlw7Fu6mOVX5OnouuZpmpPTM1iyqMtLQNOWp5httfPP4DIiUfX+yrwO293woI2FqnJcjme7mOWXL1jZ/AGdyVVN85reOzk9463yNDV4Re+16LOsQcg/Y+NgHyanZzA5PYPRiSlvIz5VDVKZQbNt6II2FqrKqWo1uU+MKSJfsLJ6AaCkozTN0YmpWR25pr+ne7b1l3kfReUh/ezgO+N4+M5OtzA//JsfMl6yqAs79x3DkkVdTgxk0zJXZkRtB3qDNhaqytmwqhcH3xnHB7+ZOwuRixuqQ5FeVAuQz+Bh1WzcogpWlu+D74zj0MgEnth7VGpmr+lKV2cMmpa5sndi+10FbSxUlbP/eKfw3LBi7izENgzD5fViu7LreGN1K1TLZuPmK1hZOv093dh+702zE8qK4kb5e01XujpjEEqZC9pYyJIvDGm/9IPfzMy5to3DcFWY6HbpeGOyK1TzsuUrWFU6dXEj295knTEIpczNC2OR7bduv/em2X5p2gVx1TfliImKotMyVt1bVYny3+nKkAY1Xx0Zl5oP04RQjEEd88JYZPute4bG8I0vXD9bQNLvZZFtiUMJlJpwgXUqg6mKpCtDGtT85elJ6RhJGaG8+6bMC2OR7bemlaK/pxuPfunGxmnJtsShBErb0urpUjcfpsn7DOXdNyVIY6FiuU1VCtmWOJSglUlCblHr4hpN3mdb372V3b2b0nTD3nSj1C3MN0idb6Tv5ea+q3DLdUvxwG0DpUYjZMNiE5t6IaIhIcQa1ftZrg3Jzs8vmqvvag1IpBkbB/uwbqAHR8bO49kDJyrXNcS1K8Vw3lKQZTekbgZidnZhkxWdpoitYjFpbGj3KyMAqt1wU2tXgPDiAlXlh3MXhqWxqJqBWLQIKB3pAOYWHBuLblwscArZCMksPNONIaWzcNOp2zaw9T5k54Rwg6WxqJqBWLQIqGpFp41FN7asf8itZYqrPKSzcPcfH7e2QE92lmhTOHsPlQghvP8MDg6KlJPjF8SOvcPi5PgFUUTd97LXF33eNG3T+H6+CVTyYPIeU+Wj6Lsde4fFZx55SezYOywtpylMlA0Ah4VGPfVuKETOWPh8IbrPbkNll8F0Pk2+86ZpNbne5/s1oSNdY8GuG+LTRVN9tsxZGj4ZnZiaDTpWDWfKYrqrYfKdN02ryfU+4wkcui5BzrPgxuyhuutXettVqoz8Mm8Tc1PaEog1Tfa4BJejc7LozrNg51mESNXqSN+ky7zTiVImWibOEXufyIzOhUzwxoJDK8e58pg2ZBz0rYutPMiMzjWFk75ZzuBsAucZbxyQOYi4CVX6DmVWpa0yk+p67UCvMZ1zKt/BexYcAj/ziSp9hzJPJKQyw0nWGOCMGIOTyxyZSysXkqnA0QXmKJMpivJmussTMm18960xFpz6dikcZTJFm/NmAlX9cDYywccsUnz27crcb079TdNdBE5584Xq6tGqxY1lRx9wIGjPImuFfbrAZa0IJ7fctCfAKW++qNJplX7K7qs6+qAOFx6JlmdBRCcAfADgEoCPhBBriKgHwN8DWAngBICvCCHO6YlZDJfoewitbFZGboFI09PRXaH63svu05kT46IumOiG/HshRPYo920A9gkhHieibcn/j9QlolKAuVRSW5OyTO4knpUxnZ4O8HB19wyNBXksg+n3rpOei7pgI2ZxD4Dbk7+/C+CfIWEsZCyj7ZOjuGFrJ3EuRjZF9ViGULHhBVTVhbTeYEHXQp1n6BoLAWAvEQkAfyWE2AXgGiHEaQAQQpwmoqtlEpIpwDZdLW6uOSBfqZvuGuXbEyqSR+VYhlBxbazTerOg+6plOunoBjg3CCFuBXAXgAeJ6POyNxLRZiI6TESHz549KxUwa7JRb1nAp+xzjkOBskHE7K5RPuGoQ06kZQ9A4+CwTgAzrTeXps6/X391OVrGQgjxbvL7DIB/ALAWwHtEtAIAkt9nSu7dJYRYI4RYs3z5cqnnNYnAV0Wciz4PbcfwtPC8OjKOyekZbFq/0rvsoelQlSYVN3utjjHVuXfWk7w0c7HxzRmUuyFEdCWATwghPkj+vhPAfwXwIoD7ATye/H5BR0BVZCLOWUKLf+SXQ2+5Y7X37lOdDjl29VRQPZ1Mp/vBIc6kvDaEiD6LjjcBdIzO3wkh/jsR9QJ4HsB1AE4CuFcIMVGVVhvWhriuCNw3WimiLYdDhXruqe7akLiQzBAyFYFTwfHBfM+/b+JOWUzwPZoTAvOlm+ITmzoMerq3b5pON9+wqhfrBnqsHowTMvNpNMXW9GybOmTjWbhoVUw/o6mn4OJgnJDhEMRzhS0vs0iHXCZlGcOFi56PTOsaDpvbzvvGR5fAxYiUSr5s3GOqLMjMajY1KYuNsXBRkbLPaGKcyl5808Id0vBsW+MrKvmycY+psiAjW1rut27Xm5TFxli4qEjZZzQxTm2tOGWMTkyxmehlGpVGydU9Ksg8Jy33WzUnZcWhUwnmW5Te1HyINusthLzlZYx7cDqgaKSD8/ZnspTlwdS07TaPbvjIW9MyZ1pGNt2Q0Ai9a1K1hZupLqFJV5xbS+4jWF1X5vI6Mi1ja42F7cLFdWRDNt86W7jJYjIOxc04+whW15W5vI5My9haY2G7cHEd2ZDNd5Mt3Di06jLGWVVODvmToa7M2W7AWmssuLb8tpHNdxNjx6FVl5F39ysjePbACfzL8Bn87z++Vbric8ifCWw3YK01Flxbftuo5lt1W3uOHBk7j92vjGDJoi4pb8FE/kLxTnQIZjTExOhDG0YwbKG6rT0nHrhtAJvWr8Sm9SsBQHokwET+bIyOcCuvwXgWuq4i9wNcsvhopULzHorI7uU5OjE161m4wIb+shscbb/3Ju/GOhhjofoy0oo3OT1jPfpvCh99aJXuC2fX23U31MbzNg72ze6EtmdozHsDx7YbknfBVF3FbMXbcsfqUgvNyeWzuZdlNp9NNzXOY9L1dq1/Tu+7jP6ebmy/9yalsmAjf2w9C1Ota7pN/l2//ztzloVnW0ZOEXGbrWI2nwAK86wy/GpSLhf65/S+q1AtCzbyx9ZYmCqIVXtImNpMNSSK8im7qXEe2YIs012peqaN7k7b37eV/AkhvP8MDg4KFU6OXxA79g6Lk+MXpK/J/i9zf0SfHXuHxWceeUns2Dvs5f4UlffdpjIC4LDQqKdsYhYqfSyZPnM+1pG9RyUOEkJfVxZXedGNwTS93+RBUlwWw3Eod2y6ISp9LB/7DITS15XBZF6qugoy53BWdTGa9tvL8uWyvGTzlcqk041SeVemu29sjIXKS5EpRKYPU25TX5dDgNKG8S3Ll8q7NxFgBIoDyU1QeVemdcvGWNgaATCtsDZNI0+XMedbH5UWSdXw2DC+HN6RTCC5CSp5Mq3b1u+UFerpUa4o2hWrLSeHRS6nFTtlXfzoY6ngjUqQx8Rhyi7IHnTsMpBVFDy0OSlMFVsBPg6Bw1Bg0Q05N3VRqqtgO7joMx6RP+gYcBNALXJvObjxeWy9+zYFrG3Dwlgs7V6ITRItWZPKrNKl8FlJ0jxlDzqO/BZbhrxNAWvbsI5Z6MQQTPe752M8I8RT4iLltCJmUYZODMF0v7upLFV94VD6yS5iOHHSk11M5otFN6QMHRfRdJciK4tMa1jVFzaxN4eL1rip/l0OuZqmrbGL/HwPLXTmipv6UV0b4guZtQpVawp01xuYWivRhKo1Nj7lMkXdOwl1jUhWbmiuDWHtWXBFpjWs8mxczyI14YnkW96iXZy4eAkq1L2TUD0Pkx52NBYK+B5aNLVWognZ0ZonXz46u09Idhcn33qxiQ8DrcvoxBR2vzICoLM/qS7RWCTUvdyy7zkUirws2eFXUy1+agjSUSYA2H7vTZctlmozPgy0LnuGxvDsgRMAgCWLurTTqzUWRLQbwH8EcEYI8bnksx4Afw9gJYATAL4ihDhHRARgJ4C7AUwB+JoQ4jVtKR2Q3wgnbwDKXj6HQpGXJT+xy2SLnzU8bfYkdOHQJds42IfJ6ZnZv7dqpifjWfw1gKcAfC/z2TYA+4QQjxPRtuT/RwDcBWB18rMOwNPJb/ZkX26RAUgVPzk9g9GJqVkj4rpQyJzvYXNiVzQQcnDQU3a3cyPIREHR8SDezPw/DGBF8vcKAMPJ338F4KtF11X9cBsNKYt8m4z2q0bXQx5xiPgFnkZDrhFCnE6MzWkiujr5/FoAo5nrxpLPTucTIKLNADYDwHXXXacohh3KWgUO+z9wcG9DhlOMSQcf+TA9g5MKPiucTy6E2CWEWCOEWLN8+XLDYtjB5DZ8qjNMQzkdjCu2Z4y6mgnqY+arqmfxHhGtSLyKFQDOJJ+PAejPXNcH4F0dAUOnyIMIpXULRc4m2PbMXAW8fXiYqsbiRQD3A3g8+f1C5vOHiOg5dAKb59Puynyl6KVyGkGpgpucJoyX7cCjq0qcz4cLwy4zdPoDALcDWEZEYwAeRcdIPE9EfwrgJIB7k8t/gs6w6XF0hk43yQrSxlYMKC6ctgtUfjKOqj7zcppKVxVuxquIJsbIZJlPdTM5PSN9enxTao2FEOKrJV/dUXCtAPCgiiAhFARTqLZuMoUrfwD0kkVdyvrMy5mf5OP6Pbkwsib2I5XF5I7dqU4mp2es1SM2Mzi5R/mrWlVXXpFM4dozNIZDIxO4ue8q3HLdUqP6zE/ycY3tLkSRfm02YkVlvq4slcmT6sbm6fFsjAWHSSxVVLWqnIJa+RmWJjExycd1692Eqh25TVa+bH6bbl9QJ4/NesTGWHCnqlX1FdRSvcYn+cqQ7zZxW6BnQ59VBsGnMagjeGPhqlWqalW5V1BO5CtD2m1aN9DDtgtqmiqDwLkssd6DU4Z0FeS6gZ7ZfRUi4cClC6JKSPK3eg9OGTYO9mHdQM/svgoR/zSZxRj6jFQue4gWYXo2afDGor+nGw/feT3WDfRgw6peq8/iuqmrCblM5i2tQN/84evsdFWETt5NbwxtEtOGLHhjAQD7j3f2b9h/fHzOdyYrwe5XRrBz37HZIVQumCgUJgtWaN6ebN6LypJpz8hkeTVtyIIPcALVASPZYU3ufU+ZfSyK8i+bL5MjOv093UHtoiWbdxfrfEwOwxsPluqsbzf1Y3M/C9l9I3R37LaN7D4WeRnj/hfm0N3RvKz8ZD+3WcYQd/euRta66u7YbRvV1s/2pCKOXpgtdNf5yG7N6GNtiRQ6lsbUz+/fdEuQZzJwxIX3M9+8FVM6lfEsmtzX9D2gDZ6F7CnqrnBlsW08x4b3k5eT+zoe0+x+ZQTPHjiByekZPPqlG5XfW9m7UT2zxPV7YDEasrR7IavhJ1dj5zrPcTmMm5cz9LkRurieW5Ef1UjfPYDa99C6s04XfvITLDyKFFcWW+c5Lpf0c/YkbHuBaSXbtH7l7EE9rvVRtFWA7LuPZ506xucoiGo/d75gO37CMT7T5N3Hs04N0KRF8rUxT9WKTM4LjlwxOjGFyekZbFq/0lorr+NF2PJ6mrx7k+WERczCB036nboz4VT7jfNxRWYT0j1GlizqkqqMKu+hKj5Tlx7ndSMqzFvPokmLoWudTZwR0uZgomoL3LTVN+0h6m5UExpsjIXrCSYu3XjVQjNfuhqqlbipfkxXXpmNajYOzj031zVp3cKCroU66bDphrTNZcsy34ca65Dt5jXtRuSvN/0eZNKrK9e2hsCz6aYyLOi+aplOmmw8i7a5bE0JZQq1z4lkTT0QG4Fp2fyn16XbJpSVa1vB82y66bO3bj//vk6abIyFze3xQyCUoxB8ytm0QWlyvWw5ks2/zHU2R3PS9Das6p3N19ZLMxd10mRjLFQxXXh9GR/bnpWpfPn0AJs2KE2uly1HsvmXuS4dzdlyx2orO7F/4wvXz247aYKgjYUNy+yr5fRxJkYVZcalrUFXWSMgm3+Z61wY3uwztmqmFbSxyFvmbD9x//FxpVa0rbET38OM3DFlBJt4cC4Mr8lnBG0siraV37nvGA6+M658DkVWuSHGQ0x5BBtW9eLgO+NS+5qGqCdbtNnIshk6VSE/dJUOwT185/WFQ3FNh6lCHM6VkVlGD1X7mqo8sw4fmyHbeCanDXxN5y9ozyJPtvVcOzC3RWxq9UPsksgG1qr00DQWZEJPPlpkG8/00Z0pw3T+gjUWKspsWqhDDOaZCKw1jdKb0JMPw8y5MTBR0U3nL9gTydIhoS13rNYqqLG/PReTOon6VcOG3nRPJAvWszBlNdsckFLFpEcV9asGR682WGORKjMN4mQt8OjE1OxBQA/cNuDsvIy2YLJV46Tf6OXowWo0pCp6W/ZdUSQ+7XM/e+BEbYQ+LvKai8m9Qev063Mv0aZwPb7SFbXGgoh2E9EZInoz89ljRHSKiI4kP3dnvvsWER0nomEi+mITYbIvM/9iyl500VDVxsE+bFq/cjaaP99fclN0hv/K3lMTY28L3WHNJrKqlDnu5VSmG/LXAJ4C8L3c508KIbZnPyCiGwDcB+BGAL8L4J+I6HohxCUZYbIuq+xhOUV9u/6ebjz6pRt/K2hmfjy3fiBHdPrLZe+Jw3b2unGAJrKqxGq4x3dqjYUQ4udEtFIyvXsAPCeE+BDACBEdB7AWwL/K3Jx9mfkXY6MAR8xT9p6aGHuuNJFVpcxxL6dSQ6eJsXhJCPG55P/HAHwNwCSAwwAeFkKcI6KnABwUQvxtct0zAH4qhNhTlb7K0KkP5lOAbD7lVRcfulJ5pu7QqWqA82kAvwfgZgCnATyRylNwbaE1IqLNRHSYiA6fPXtWUQy3hDj9WxUueeXejwf86MrHM5WGToUQ76V/E9F3ALyU/DsGoD9zaR+Ad0vS2AVgF9DxLFTkqMKGtS9yE8tWutY9v+x7Li06F5eYez8emD+zT5WMBRGtEEKcTv79MoB0pORFAH9HRDvQCXCuBvCqtpQKuJr3X7bSte75Zd9zqRxcYglcjFYVPnTl45m1xoKIfgDgdgDLiGgMwKMAbieim9HpYpwA8OcAIIR4i4ieB/BLAB8BeFB2JMQ0rgpZdvuy1LOQeX7Z9yFUDpPUeVJcjFYk4LUhQNinnZvAhVy2n2FqjU+knnm5NiQtwJPTM3j2wAkAdl12Ll2DPC7ksv2M+eZJNYFbIxWksUgL8Kb1K51sNMK1QLvew9EGIXUzXFdebo1UkMYiW4A5nV4W0qlqsrJWPYNby2cb15WXWyMVpLHg2hrlD3bhWpGyp7MffGcc2++9SUnGosrTZgNSV3lN551bOWe16jR0sguVuExqKmLPUOd09msWX4FDIxOzy/mbUrQwqyjfXCZW6cpRt4KW8zs3QZCeBVeq1rbYpGmLlsp06twU9rx2Svm5RS1fUb5V3XfTLXVbg7WuvLlWGAuOrq9LF7JpJchuHHTt0u5GhVtlXoRqJTJdudsarHUVS2mFseAWNXaNaiVQKdwqulatRNwCfLZRbfRc6akVxqKJsjh6IbJwOFIwr2tb+rSRLvdGRVU+V++fTYBTNviUv65poWoahOISnAPmyu5DtnyQz9ShRnnq0lVJc8OqXqwb6JE6Zc0Hujt52YaNZyFrVfPX2T44iFNrlJc9u4hNdfgT0GvFZfSposO6dFXSzJ6ytnagl52XyW2odA5CCO8/g4OD4uT4BbFj77A4OX7hsr/z5L+rutYEptKXTafJ806OXxBf+T8HxGceeUns2DusLNuOvcPaaVRh4x2ppJm/x3a+uQHgsNCop94NhUiMRZY2vkTZPDXNu2pFlDHOtg2xb5o0Sm1A11iw6YZkaWMUXDZPTfOu6rrm3fiiGZiuu2Ccpstz6n5ygaWxqDtAiFM/U5Ymayx8rjvIT1kvukaHqvdnooJWpd+k7LSxwdKFpbFIKSo8bWztbFeSImR24bZhuKryaqKCVqXfRM/sg40eYG0sigqPa4vvwjiZrCST0zNYsqhL2bjZriRVeW3y7DLjWJV+iPNxuMgBBL5Tlgs4vawq8hsC2dx5ioNObO+wxWUHL5NyzMudslzCyR2tqqTZOE/qWdiCQ/DPtoepk76rg6VdG+1oLAJCppLaMm7ZgllUgDmNZPhO36Qx5TRi0wpjwcEtdoHPCH3ZUGvZ923DxUiK6lYDrspDK4xF2wtqis8uUV3BbPtQo6mRFJNDx67LQyuMRdsLKgeq5r5kv28rpsqY7aFjm7BZdapD3XZnJuC0+tQFZfltuhpUR2+cdG6qjFWtLHVRjnVohbFwAff9FU1XrLL81i2jzt+nozfuOleBu0GoohXdEBc0dRFDP2OiLL913Y38fTquNXe3fL4RJ2UpIGMIXE/q4TAixEEGG6jki6Mu4qQsD8i04q5bRQ4BxraOSqnkq426iMZCARlDwKHy2qCqxWxjt2F0YgqT0zPYtH5lo3xxmQFqkhjgVKC/p3t2vweZgGJR8NFVpN9m4DOftk7wzoU+VPcCffbACSxZ1NUoXzq64BrYjZ6FIk3cTJ9L7W0GPk2m7UIfNvYCtQFXDy0aC0WavFCfS+1NP8fWqWsu9KHyDB/dSa5d2DgaMs/h2j9WhWt+OMilOxoSYxbzHK79Y1W4ngtjS88uZ7nGbohFOLQmdaQu+YZVvXjy5aPYsKoX+4+Ps5a5iqZdDVexI91dusrKkssh2lrPgoj6iehnRPQ2Eb1FRFuSz3uI6GWzZRDHAAAJZklEQVQiOpb8Xpp8TkT0F0R0nIjeIKJbbQhucg2CLXa/MoKd+45h9ysjvkUpJe0f7z8+jp37juGJvUeD9jSajkJsHOzDpvUrMTk9Y7XsNJGryAtRnX5vEhnP4iMADwshXiOixQCGiOhlAF8DsE8I8TgRbQOwDcAjAO4CsDr5WQfg6eS3UVRPJjPR2ofgMTQl62GknoUJuOuqv6cbSxZ1Yee+Y1iyqGt2Ze2eoTFlL0s3z00C4i6DobXGQghxGsDp5O8PiOhtANcCuAfA7cll3wXwz+gYi3sAfC851OQgEX2aiFYk6RhDdQ2CCbctn0ZZ4XjgtgHrW9wBZipkttCtHbj8LNCy9GWeG8JMxnzZyR4LeWhkAkAz2XXzXGQAOIyQNIpZENFKALcAOATgmtQACCFOE9HVyWXXAhjN3DaWfHaZsSCizQA2A8B1113XWPC88mSVaWKIrqxwAVCSSRfbFbIo/dGJKXzzh6/XVqYNq3px8J3x0sOIOXge2Ul22W0DVb0slTLGQQ91SBsLIvoUgB8B+LoQYpKISi8t+GzO+KwQYheAXUBn6FRWjhRVV9FEBc6n4XsSje3nF6W/Z2gMh0YmsG6gp/K5+cOI8/g+Byb9/9S5Kex57RQmp2fwwG0DAIAVVy0y5hnUUacHGWNi2+BIGQsi6kLHUHxfCPHj5OP30u4FEa0AcCb5fAxAf+b2PgDvmhI4RddVBMwpN1845sPmtVkDUpVHbtvxlcW6bu67qvSaKuretWxZqNODjEzZazYO9mH3KyP44MMZLL6ia9YA6lBrLKjjQjwD4G0hxI7MVy8CuB/A48nvFzKfP0REz6ET2DxvOl4BmAnI2WrVbKVr0wiVtbhlz5I1UHXXue6Ll8W4isqRTJmqe9eyZaHpPiF116RrWlKWLOoqz4QsdScnA7gNnW7EGwCOJD93A+gFsA/AseR3T3I9AfhLAL8C8AsAa+qekT9F3RW2Tsq2la7N0+XzabfxJHsb1L3rpmXBVNk5OX5BPPbCm+Lh5/+feOyFN8XJ8Qvap6jH6d4BwcmziNjB5qZJutO9o7GIBEebDZnNvMW1IRH2uNpM2AYuZwZzN4LRWESMUlS5TFfu7BRn25XZpWHivqgvLiSLGKUo+p+P5Ou2oNmRg7SPn32eSUwP7epuS+jT+4jGYp5juvAVFfj8sKDJoeWi55nMU5OhXd3p7zLP8jl9np2x4Npv4yqXLmWFTzW/MgXeZGtdNCFOZhq6LE304GLX942DfZicnpldJeuyLLKLWXDtt5mQi+My+rIlzjbfg81TuWSnoTdJT1YPMsvFm+S9qLykq2SfPXDCeR1h51n4XmdRhgm5bLqQpj0Bru+hDtlp6Crp1WF6NmpZefH1buI8C4fY7Mq4PgEtYh/T5SVOyooAaG9MZT5h+x3GSVkRAHP7whzjI6ES+qa+pojGoqVwL3ghUaVLk4ZEJkDaBNNGjl2AM2KGUAOUHKnSpcm9L1QDpK52/o7GwgAc4wVVBY+jvJyp0mUTo+x6/xTTDca8MxY2KorpQmC7MoewiS5Hit5LE2/AlrdXlq7podx5ZyxsVBTThcCEjLprEELBlmEtStfGrt0mcLXbWNDGQqWg1O02rYLpl2V7AhiHbeVN4dK1b5ORVSHo0RCViH92t+k6fA0/1k0JlpFLJ7JuI986J8hVXWt6BKEqXZvT1EMgaM9CxdJzCEjpIiOXjvdgI99lu2rLPMOHl9Qm7yv1wLGga6FWQjobeJr6yW7Ya2uzWxU4yZIlK5cNGZukKXtt/rr0/0PvvF97v833YCLtojTK8uujLKWbLy9YvOxdoVFP2XVDfE0mKlvhx9HtzMplQ1+6h/jKpJk/kLnqfturVHX1V5RG/jOfk+TSLtWlqfPv66TDrhviK4hU5eqOTkzNnoSeHtbCZZ6C76Cbif0ZdO7XxcTzi9IoO5/ERz5TY7v10sxFnXTiQrKEqpGV7NZtW+5YDQBxhSdT2jThjNuqU3aehS/qZulNTs/M/p39PMILrkFpFbjlJRoLCfp7uvHol2687DMOLy8yF9/dGhlMnX/qGnYBThlczX/wMc/C9DPn21J1rkHpLEXBzqJ5KNy6U0EaC93IsmwF8hHBNv3M0Jaq+zJuLp9bNOGryeiJLx0F2Q3Rdc9k+4I+3EDTz+Tmytbhq5/u8rlF8bEmoye+dDQvR0M4ungpnGVzga/8h6R3VVnjtnoKcO7XhtZtMI2Jd5N109sYs/G1hWKQ3ZA2E1q3IUva4m1Y1Yv9x8e1Wmmdlj7rpgOQctm5DVM2wZXs0Vgwg+MCJtmKmxbafxk+gyNj5zE5PTNnyFkWnQpQNaOyyT2h4Er2aCwQVn/VB00DwqfOTeHI2PnLvmuqY50KkDe4MsaGo5GWxZXs8zJmkcfmUGwb5mrI7hmRFtot/+F6bLlj9ew6GkBOx1m5bcSVfMcvfD9fl1pjQUT9RPQzInqbiN4ioi3J548R0SkiOpL83J2551tEdJyIhonoizYzYALdDVSqKkIb5mo0rbhF18vo2LaufAePfT9fF5luyEcAHhZCvEZEiwEMEdHLyXdPCiG2Zy8mohsA3AfgRgC/C+CfiOh6IcQlk4Ln0elK6LpxVS5z2Xc2uz4c+982Tld32bUxge/n69J4ngURvQDgKQAbAPy6wFh8CwCEEP8j+f8fATwmhPjXsjRNzLMI7azP0OTlSNRhM5yedUpEKwH8HMDnAGwF8DUAkwAOo+N9nCOipwAcFEL8bXLPMwB+KoTYU5HuWQAXAKhvzrGga+GC7quWXZo6/z401+1LsAw6sgLhyesOeVnd6rCMkHT774QQi1Vvlh4NIaJPAfgRgK8LISaJ6GkA/w2ASH4/AeABAFRw+xyLRESbAWxO/v02gM06Vs8lRHQ4FFmBsOQNSVYgLHmJSMt9lxoNIaIudAzF94UQPwYAIcR7QohLQoiPAXwHwNrk8jEA/Znb+wC8m09TCLFLCLEm+dmlk4lIJGIfmdEQAvAMgLeFEDsyn6/IXPZlAG8mf78I4D4iuoKIBgCsBvCqOZEjkYgPZLohGwD8CYBfENGR5LP/AuCrRHQzOl2MEwD+HACEEG8R0fMAfonOSMqDkiMhIXkXIckKhCVvSLICYcmrJSuLVaeRSIQ/cQZnJBKRIhqLSCQiRTQWkUhEimgsIpGIFNFYRCIRKaKxiEQiUkRjEYlEpIjGIhKJSPH/AVCpMvbyj1c/AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.spy(dummy, markersize=1);" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " (137, 129)\t0.05409563715292309\n", " (219, 194)\t0.12676575598664064\n", " (236, 77)\t0.1297133685024946\n", " (170, 78)\t0.9913235653850567\n", " (129, 211)\t0.8319650607418313\n", " (179, 121)\t0.4852727096924244\n", " (50, 38)\t0.7331880317797461\n", " (129, 249)\t0.07947388309988046\n", " (173, 169)\t0.8207981993058538\n", " (46, 241)\t0.675461135234256\n", " (64, 136)\t0.21879256797804525\n", " (84, 110)\t0.03396910126513719\n", " (118, 214)\t0.6867001717374005\n", " (101, 236)\t0.6995301472373171\n", " (191, 101)\t0.5915974314523315\n", " (75, 181)\t0.8773385273151388\n", " (25, 14)\t0.09964382302244934\n", " (198, 137)\t0.4576743324358349\n", " (59, 2)\t0.6872116290619612\n", " (162, 138)\t0.005728883655892636\n", " (79, 17)\t0.5707762190206077\n", " (109, 142)\t0.1411587614633989\n", " (224, 31)\t0.41961256130139424\n", " (226, 17)\t0.897923863750546\n", " (89, 118)\t0.6354325627216016\n", " :\t:\n", " (229, 4)\t0.8557838376021782\n", " (223, 56)\t0.42981863607576776\n", " (167, 13)\t0.1363390090049006\n", " (176, 229)\t0.6047778926232239\n", " (44, 223)\t0.8098997267890242\n", " (16, 33)\t0.6915070379679459\n", " (106, 193)\t0.949574956031497\n", " (48, 92)\t0.5688923600606588\n", " (196, 63)\t0.034631305534336354\n", " (133, 219)\t0.6168090033350935\n", " (189, 192)\t0.17081354791601855\n", " (122, 72)\t0.5787819602985719\n", " (194, 213)\t0.06466836426599487\n", " (58, 67)\t0.16853843827033066\n", " (33, 58)\t0.6559115059826436\n", " (53, 198)\t0.5550195847459326\n", " (106, 85)\t0.20177755561196198\n", " (75, 115)\t0.46980875381404874\n", " (171, 107)\t0.45934086279992803\n", " (86, 75)\t0.29816299296978654\n", " (86, 77)\t0.7790217312936337\n", " (115, 49)\t0.497554416365638\n", " (70, 168)\t0.39036906182339637\n", " (122, 120)\t0.8379818823021193\n", " (150, 116)\t0.7450357965166439\n" ] } ], "source": [ "print(dummy)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us convert the sparse array to some dense formats." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "numpy.ndarray" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "type(dummy.A)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "numpy.ndarray" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "type(dummy.toarray())" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "numpy.matrixlib.defmatrix.matrix" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "type(dummy.todense())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. Sparse matrix storage schemes\n", "\n", "The `scipy.sparse` module provides several formats to store sparse matrices.\n", "Each format has pros and cons, and some are better for some tasks, such as matrix construction, indexing, or linear operations." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.1 List of lists format (LIL)\n", "\n", "* Supports indexing, which cannot be done with other sparse matrix formats.\n", "* Changing sparsity structure is efficient, e.g., reading a sparse matrix from a text file." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "# Create an empty lil matrix.\n", "mtx = sparse.lil_matrix((4, 5))" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "# Assign some of the indices, i.e., changing the sparsity.\n", "mtx[:2, [1, 3]] = np.array([[1, 2], [3, 4]])" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[0., 1., 0., 2., 0.],\n", " [0., 3., 0., 4., 0.],\n", " [0., 0., 0., 0., 0.],\n", " [0., 0., 0., 0., 0.]])" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mtx.toarray()" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[0., 1., 0., 2., 0.],\n", " [0., 3., 0., 4., 0.]])" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Read some of the indices.\n", "mtx[:2].toarray()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.2 Coordinate format (COO)\n", "\n", "A COO matrix is constructed from three lists:\n", "* a list of column indices,\n", "* a list of row indices,\n", "* a list of values,\n", "where each element of those lists represents a non-zero element in the resulting sparse matrix.\n", "\n", "This format is well-adapted to build a sparse adjacency matrix from an edge list." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "row = np.array([0, 3, 1, 0]) # row coordinates\n", "col = np.array([0, 3, 1, 2]) # column coordinates\n", "data = np.array([4, 5, 7, 9]) # values\n", "\n", "mtx = sparse.coo_matrix((data, (row, col)), shape=(4, 4))" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[4, 0, 9, 0],\n", " [0, 7, 0, 0],\n", " [0, 0, 0, 0],\n", " [0, 0, 0, 5]])" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mtx.toarray()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Advantages:\n", "* Fast element-wise operations.\n", "* Fast conversion to other sparse formats." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[2. , 0. , 3. , 0. ],\n", " [0. , 2.64575131, 0. , 0. ],\n", " [0. , 0. , 0. , 0. ],\n", " [0. , 0. , 0. , 2.23606798]])" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Element-wise power.\n", "mtx.power(0.5).toarray()" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "mtx_csr = mtx.tocsr()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Disadvantages:\n", "* Indexing is not possible. (Use LIL instead!)\n", "* Slow at arithmetic operations. (Use CSR, CSC instead!)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Exercise:** Can you construct the sparse adjacency matrix in `COO` and `LIL` formats for a network given by the following edge list ?" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
node_1node_2weights
0130.6
1140.5
2150.7
3260.1
4340.6
5350.1
6360.9
\n", "
" ], "text/plain": [ " node_1 node_2 weights\n", "0 1 3 0.6\n", "1 1 4 0.5\n", "2 1 5 0.7\n", "3 2 6 0.1\n", "4 3 4 0.6\n", "5 3 5 0.1\n", "6 3 6 0.9" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "edges = pd.DataFrame(\n", " {\"node_1\": [1,1,1,2,3,3,3],\n", " \"node_2\": [3,4,5,6,4,5,6],\n", " \"weights\": [0.6,0.5,0.7,0.1,0.6,0.1,0.9]\n", " })\n", "edges" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.3 Compressed sparse row & column formats (CSR & CSC)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([4, 9, 7, 5], dtype=int64)" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Get the data array\n", "mtx_csr.data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`CSR` is row oriented:\n", "* efficient row slicing\n", "* fast matrix vector products, the right multiplication `CSR * v`" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([0, 2, 1, 3], dtype=int32)" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Get array of column indices for CSR.\n", "mtx_csr.indices" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([13, 7, 0, 5], dtype=int64)" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Matrix-vector product from the right.\n", "v = np.array([1, 1, 1, 1])\n", "mtx_csr.dot(v)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`CSC` is column oriented:\n", "* efficient column slicing\n", "* fast matrix vector products, the left multiplication `v * CSC`" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([0, 1, 0, 3], dtype=int32)" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mtx_csc = mtx.tocsc()\n", "# Get array of row indices for CSC\n", "mtx_csc.indices" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([4, 7, 9, 5], dtype=int64)" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# vectro-matrix product\n", "v * mtx_csc" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Efficient arithmetic operations `CSC + CSC`, `CSR * CSR`, etc." ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[16, 0, 36, 0],\n", " [ 0, 49, 0, 0],\n", " [ 0, 0, 0, 0],\n", " [ 0, 0, 0, 25]], dtype=int64)" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Matrix-Matrix product (* is elementwise product on Numpy!)\n", "prod = mtx_csc * mtx_csc\n", "prod.toarray()" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[16, 0, 36, 0],\n", " [ 0, 49, 0, 0],\n", " [ 0, 0, 0, 0],\n", " [ 0, 0, 0, 25]], dtype=int64)" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "prod = mtx_csr @ mtx_csr # @ is matrix product both on numpy and scipy!\n", "prod.toarray()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can read more about sparse matrix storage schemes [on Wikipedia](https://en.wikipedia.org/wiki/Sparse_matrix#Storing_a_sparse_matrix)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3. Linear agebra on sparse matrices" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3.1 Some basic operations" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[1, 0, 3, 0],\n", " [1, 2, 0, 4],\n", " [0, 2, 3, 0],\n", " [0, 0, 3, 4]])" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# sparse matrix from diagonals\n", "A = sparse.spdiags(np.array([[1,2,3,4], [1,2,3,4], [1,2,3,4]]), [-1,0,2], 4, 4)\n", "A.toarray()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Inversion of a sparse matrix**" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 0.66666667, 0.33333333, -0.33333333, -0.33333333],\n", " [-0.16666667, 0.16666667, 0.33333333, -0.16666667],\n", " [ 0.11111111, -0.11111111, 0.11111111, 0.11111111],\n", " [-0.08333333, 0.08333333, -0.08333333, 0.16666667]])" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A = A.tocsc() # Convert it to CSC matrix for efficiency.\n", "Ainv = sparse.linalg.inv(A)\n", "Ainv.toarray()" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "8.306623862918075" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sparse.linalg.norm(A) # Default to Frobenius norm." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Solve $A x = b$**" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([0.33333333, 0.16666667, 0.22222222, 0.08333333])" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "b = np.array([1, 1, 1, 1])\n", "x = sparse.linalg.spsolve(A, b)\n", "x" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3.2 Eigenvalue decomposition" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For the full eigendecomposition of an array, you can use the functions provided by Numpy:\n", "* `numpy.linalg.eig`\n", "* `numpy.linalg.eigvals`\n", "* `numpy.linalg.eigh`\n", "* `numpy.linalg.eighvals`\n", "\n", "Scipy presents more functionality (read [here](https://www.scipy.org/scipylib/faq.html#why-both-numpy-linalg-and-scipy-linalg-what-s-the-difference)) such as solving generalized eigenvalue problem, you can use the functions from Scipy:\n", "* `scipy.linalg.eig`\n", "* `scipy.linalg.eigvals`\n", "* `scipy.linalg.eigh`\n", "* `scipy.linalg.eighvals`" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([1.16822694+2.48096949j, 1.16822694-2.48096949j,\n", " 1.57169108+0.j , 6.09185505+0.j ])" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "linalg.eigvals(A.toarray())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Decomposition of an Hermitian matrix:" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([0.17157288, 5.82842712])" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A = np.array([[1, -2j], [2j, 5]])\n", "linalg.eigvalsh(A)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "However, for quickly finding a few eigenvalues of a large sparse matrix, you should use the corresponding functions from the [sparse module](https://docs.scipy.org/doc/scipy/reference/tutorial/arpack.html):\n", "\n", "* `scipy.sparse.eigs`\n", "* `scipy.sparse.eigsh`" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 2.77793090e-18+0.00000000e+00j, -4.94872981e-17+0.00000000e+00j,\n", " 7.58308618e-19+0.00000000e+00j, 1.73472373e-18+1.74175441e-17j,\n", " 1.73472373e-18-1.74175441e-17j])" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dummy = sparse.random(30, 30, density=0.01)\n", "evals, evecs = sparse.linalg.eigs(dummy, k=5, which='SM')\n", "evals" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.6" } }, "nbformat": 4, "nbformat_minor": 2 }