4.scipy/8.scipy-sparse.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import pylab as pl\n",
    "from scipy import sparse\n",
    "from scipy.sparse import csgraph"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 稀疏矩阵-sparse"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 稀疏矩阵的储存形式"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dict_keys([(2, 3), (3, 3), (4, 3)])\n",
      "dict_values([1.0, 2.0, 3.0])\n"
     ]
    }
   ],
   "source": [
    "from scipy import sparse\n",
    "a = sparse.dok_matrix((10, 5))\n",
    "\n",
    "a[2, 3] = 1.0\n",
    "a[3, 3] = 2.0\n",
    "a[4, 3] = 3.0\n",
    "print(a.keys())\n",
    "print(a.values())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[list([]) list([]) list([1.0]) list([3.0, 2.0]) list([]) list([]) list([])\n",
      " list([]) list([]) list([])]\n",
      "[list([]) list([]) list([3]) list([2, 4]) list([]) list([]) list([])\n",
      " list([]) list([]) list([])]\n"
     ]
    }
   ],
   "source": [
    "b = sparse.lil_matrix((10, 5))\n",
    "b[2, 3] = 1.0\n",
    "b[3, 4] = 2.0\n",
    "b[3, 2] = 3.0\n",
    "print(b.data)\n",
    "print(b.rows)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[3 4 2 3] [2 3 3 2] [ 1  2  3 10]\n",
      "[[ 0  0  0  0  0  0]\n",
      " [ 0  0  0  0  0  0]\n",
      " [ 0  0  0 11  0  0]\n",
      " [ 0  0  3  0  2  0]\n",
      " [ 0  0  0  0  0  0]]\n"
     ]
    }
   ],
   "source": [
    "row = [2, 3, 3, 2]\n",
    "col = [3, 4, 2, 3]\n",
    "data = [1, 2, 3, 10]\n",
    "c = sparse.coo_matrix((data, (row, col)), shape=(5, 6))\n",
    "print (c.col, c.row, c.data)\n",
    "print (c.toarray())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 矩阵向量相乘"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1, -3, -1], dtype=int32)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    " import numpy as np\n",
    "from scipy.sparse import csr_matrix\n",
    "A = csr_matrix([[1, 2, 0], [0, 0, 3], [4, 0, 5]])\n",
    "v = np.array([1, 0, -1])\n",
    "A.dot(v)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 示例1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "构造一个1000x1000 lil_matrix并添加值："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.sparse import lil_matrix\n",
    "from scipy.sparse.linalg import spsolve\n",
    "from numpy.linalg import solve, norm\n",
    "from numpy.random import rand"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = lil_matrix((1000, 1000))\n",
    "A[0, :100] = rand(100)\n",
    "A[1, 100:200] = A[0, :100]\n",
    "A.setdiag(rand(1000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "现在将其转换为CSR格式，并求解$A x = b$的$x$："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = A.tocsr()\n",
    "b = rand(1000)\n",
    "x = spsolve(A, b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "将其转换为密集矩阵并求解，并检查结果是否相同："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_ = solve(A.toarray(), b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "现在我们可以使用以下公式计算错误的范数："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "err = norm(x-x_)\n",
    "err < 1e-10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 示例2\n",
    "构造COO格式的矩阵："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import sparse\n",
    "from numpy import array\n",
    "I = array([0,3,1,0])\n",
    "J = array([0,3,1,2])\n",
    "V = array([4,5,7,9])\n",
    "A = sparse.coo_matrix((V,(I,J)),shape=(4,4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "注意，索引不需要排序。\n",
    "\n",
    "转换为CSR或CSC时，将对重复的（i，j）条目进行求和。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "I = array([0,0,1,3,1,0,0])\n",
    "J = array([0,2,1,3,1,0,0])\n",
    "V = array([1,1,1,1,1,1,1])\n",
    "B = sparse.coo_matrix((V,(I,J)),shape=(4,4)).tocsr()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这对于构造有限元刚度矩阵和质量矩阵是有用的。"
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}