04-scipy/04.09-linear-algbra.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 线性代数"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`numpy` 和 `scipy` 中，负责进行线性代数部分计算的模块叫做 `linalg`。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.linalg\n",
    "import scipy as sp\n",
    "import scipy.linalg\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import linalg\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## numpy.linalg VS scipy.linalg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "一方面`scipy.linalg` 包含 `numpy.linalg` 中的所有函数，同时还包含了很多 `numpy.linalg` 中没有的函数。\n",
    "\n",
    "另一方面，`scipy.linalg` 能够保证这些函数使用 BLAS/LAPACK 加速，而 `numpy.linalg` 中这些加速是可选的。\n",
    "\n",
    "因此，在使用时，我们一般使用 `scipy.linalg` 而不是 `numpy.linalg`。\n",
    "\n",
    "我们可以简单看看两个模块的差异："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number of items in numpy.linalg: 36\n",
      "number of items in scipy.linalg: 115\n"
     ]
    }
   ],
   "source": [
    "print \"number of items in numpy.linalg:\", len(dir(numpy.linalg))\n",
    "print \"number of items in scipy.linalg:\", len(dir(scipy.linalg))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## numpy.matrix VS 2D numpy.ndarray"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "线性代数的基本操作对象是矩阵，而矩阵的表示方法主要有两种：`numpy.matrix` 和 2D `numpy.ndarray`。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### numpy.matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`numpy.matrix` 是一个矩阵类，提供了一些方便的矩阵操作：\n",
    "- 支持类似 `MATLAB` 创建矩阵的语法\n",
    "- 矩阵乘法默认用 `*` 号\n",
    "- `.I` 表示逆，`.T` 表示转置\n",
    "\n",
    "可以用 `mat` 或者 `matrix` 来产生矩阵："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "matrix([[1, 2],\n",
      "        [3, 4]])\n",
      "matrix([[1, 2],\n",
      "        [3, 4]])\n"
     ]
    }
   ],
   "source": [
    "A = np.mat(\"[1, 2; 3, 4]\")\n",
    "print repr(A)\n",
    "\n",
    "A = np.matrix(\"[1, 2; 3, 4]\")\n",
    "print repr(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "转置和逆："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "matrix([[-2. ,  1. ],\n",
      "        [ 1.5, -0.5]])\n",
      "matrix([[1, 3],\n",
      "        [2, 4]])\n"
     ]
    }
   ],
   "source": [
    "print repr(A.I)\n",
    "print repr(A.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "矩阵乘法："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "matrix([[17],\n",
      "        [39]])\n"
     ]
    }
   ],
   "source": [
    "b = np.mat('[5; 6]')\n",
    "print repr(A * b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2 维 numpy.ndarray"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "虽然 `numpy.matrix` 有着上面的好处，但是一般不建议使用，而是用 2 维 `numpy.ndarray` 对象替代，这样可以避免一些不必要的困惑。\n",
    "\n",
    "我们可以使用 `array` 复现上面的操作："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "array([[1, 2],\n",
      "       [3, 4]])\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[1,2], [3,4]])\n",
    "print repr(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "逆和转置："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "array([[-2. ,  1. ],\n",
      "       [ 1.5, -0.5]])\n",
      "array([[1, 3],\n",
      "       [2, 4]])\n"
     ]
    }
   ],
   "source": [
    "print repr(linalg.inv(A))\n",
    "print repr(A.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "矩阵乘法："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "array([17, 39])\n"
     ]
    }
   ],
   "source": [
    "b = np.array([5, 6])\n",
    "\n",
    "print repr(A.dot(b))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "普通乘法："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "array([[ 5, 12],\n",
      "       [15, 24]])\n"
     ]
    }
   ],
   "source": [
    "print repr(A * b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`scipy.linalg` 的操作可以作用到两种类型的对象上，没有区别。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 基本操作"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 求逆"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "矩阵 $\\mathbf{A}$ 的逆 $\\mathbf{B}$ 满足：$\\mathbf{BA}=\\mathbf{AB}=I$，记作 $\\mathbf{B} = \\mathbf{A}^{-1}$。\n",
    "\n",
    "事实上，我们已经见过求逆的操作，`linalg.inv` 可以求一个可逆矩阵的逆："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-2.   1. ]\n",
      " [ 1.5 -0.5]]\n",
      "[[  1.00000000e+00   0.00000000e+00]\n",
      " [  8.88178420e-16   1.00000000e+00]]\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[1,2],[3,4]])\n",
    "\n",
    "print linalg.inv(A)\n",
    "\n",
    "print A.dot(scipy.linalg.inv(A))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 求解线性方程组"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "例如，下列方程组\n",
    "$$\n",
    "\\begin{eqnarray*} \n",
    "x + 3y + 5z & = & 10 \\\\\n",
    "2x + 5y + z & = & 8  \\\\\n",
    "2x + 3y + 8z & = & 3\n",
    "\\end{eqnarray*}\n",
    "$$\n",
    "的解为：\n",
    "$$\n",
    "\\begin{split}\\left[\\begin{array}{c} x\\\\ y\\\\ z\\end{array}\\right]=\\left[\\begin{array}{ccc} 1 & 3 & 5\\\\ 2 & 5 & 1\\\\ 2 & 3 & 8\\end{array}\\right]^{-1}\\left[\\begin{array}{c} 10\\\\ 8\\\\ 3\\end{array}\\right]=\\frac{1}{25}\\left[\\begin{array}{c} -232\\\\ 129\\\\ 19\\end{array}\\right]=\\left[\\begin{array}{c} -9.28\\\\ 5.16\\\\ 0.76\\end{array}\\right].\\end{split}\n",
    "$$\n",
    "\n",
    "我们可以使用 `linalg.solve` 求解方程组，也可以先求逆再相乘，两者中 `solve` 比较快。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-9.28  5.16  0.76]\n",
      "[  0.00000000e+00  -1.77635684e-15  -8.88178420e-16]\n",
      "inv and dot: 0.0353579521179 s\n",
      "[-9.28  5.16  0.76]\n",
      "[  0.00000000e+00  -1.77635684e-15  -1.77635684e-15]\n",
      "solve: 0.0284671783447 s\n"
     ]
    }
   ],
   "source": [
    "import time\n",
    "\n",
    "A = np.array([[1, 3, 5],\n",
    "              [2, 5, 1],\n",
    "              [2, 3, 8]])\n",
    "b = np.array([10, 8, 3])\n",
    "\n",
    "tic = time.time()\n",
    "\n",
    "for i in xrange(1000):\n",
    "    x = linalg.inv(A).dot(b)\n",
    "\n",
    "print x\n",
    "print A.dot(x)-b\n",
    "print \"inv and dot: {} s\".format(time.time() - tic)\n",
    "\n",
    "tic = time.time()\n",
    "\n",
    "for i in xrange(1000):\n",
    "    x = linalg.solve(A, b)\n",
    "\n",
    "print x\n",
    "print A.dot(x)-b\n",
    "print \"solve: {} s\".format(time.time() - tic)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 计算行列式"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "方阵的行列式为\n",
    "$$\n",
    "\\left|\\mathbf{A}\\right|=\\sum_{j}\\left(-1\\right)^{i+j}a_{ij}M_{ij}.\n",
    "$$\n",
    "\n",
    "其中 $a_{ij}$ 表示 $\\mathbf{A}$ 的第 $i$ 行 第 $j$ 列的元素，$M_{ij}$ 表示矩阵 $\\mathbf{A}$ 去掉第 $i$ 行 第 $j$ 列的新矩阵的行列式。\n",
    "\n",
    "例如，矩阵\n",
    "$$\n",
    "\\begin{split}\\mathbf{A=}\\left[\\begin{array}{ccc} 1 & 3 & 5\\\\ 2 & 5 & 1\\\\ 2 & 3 & 8\\end{array}\\right]\\end{split}\n",
    "$$\n",
    "的行列式是：\n",
    "$$\n",
    "\\begin{eqnarray*} \\left|\\mathbf{A}\\right| & = & 1\\left|\\begin{array}{cc} 5 & 1\\\\ 3 & 8\\end{array}\\right|-3\\left|\\begin{array}{cc} 2 & 1\\\\ 2 & 8\\end{array}\\right|+5\\left|\\begin{array}{cc} 2 & 5\\\\ 2 & 3\\end{array}\\right|\\\\  & = & 1\\left(5\\cdot8-3\\cdot1\\right)-3\\left(2\\cdot8-2\\cdot1\\right)+5\\left(2\\cdot3-2\\cdot5\\right)=-25.\\end{eqnarray*}\n",
    "$$\n",
    "\n",
    "可以用 `linalg.det` 计算行列式："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-25.0\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[1, 3, 5],\n",
    "              [2, 5, 1],\n",
    "              [2, 3, 8]])\n",
    "\n",
    "print linalg.det(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 计算矩阵或向量的模"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "矩阵的模定义如下：\n",
    "$$\n",
    "\\begin{split}\\left\\Vert \\mathbf{A}\\right\\Vert =\\left\\{ \\begin{array}{cc} \\max_{i}\\sum_{j}\\left|a_{ij}\\right| & \\textrm{ord}=\\textrm{inf}\\\\ \\min_{i}\\sum_{j}\\left|a_{ij}\\right| & \\textrm{ord}=-\\textrm{inf}\\\\ \\max_{j}\\sum_{i}\\left|a_{ij}\\right| & \\textrm{ord}=1\\\\ \\min_{j}\\sum_{i}\\left|a_{ij}\\right| & \\textrm{ord}=-1\\\\ \\max\\sigma_{i} & \\textrm{ord}=2\\\\ \\min\\sigma_{i} & \\textrm{ord}=-2\\\\ \\sqrt{\\textrm{trace}\\left(\\mathbf{A}^{H}\\mathbf{A}\\right)} & \\textrm{ord}=\\textrm{'fro'}\\end{array}\\right.\\end{split}\n",
    "$$\n",
    "其中，$\\sigma_i$ 是矩阵的奇异值。\n",
    "\n",
    "向量的模定义如下：\n",
    "$$\n",
    "\\begin{split}\\left\\Vert \\mathbf{x}\\right\\Vert =\\left\\{ \\begin{array}{cc} \\max\\left|x_{i}\\right| & \\textrm{ord}=\\textrm{inf}\\\\ \\min\\left|x_{i}\\right| & \\textrm{ord}=-\\textrm{inf}\\\\ \\left(\\sum_{i}\\left|x_{i}\\right|^{\\textrm{ord}}\\right)^{1/\\textrm{ord}} & \\left|\\textrm{ord}\\right|<\\infty.\\end{array}\\right.\\end{split}\n",
    "$$\n",
    "\n",
    "`linalg.norm` 可以计算向量或者矩阵的模："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5.47722557505\n",
      "5.47722557505\n",
      "6\n",
      "4\n",
      "7\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[1, 2],\n",
    "              [3, 4]])\n",
    "\n",
    "print linalg.norm(A)\n",
    "\n",
    "print linalg.norm(A,'fro') # frobenius norm 默认值\n",
    "\n",
    "print linalg.norm(A,1) # L1 norm 最大列和\n",
    "\n",
    "print linalg.norm(A,-1) # L -1 norm 最小列和\n",
    "\n",
    "print linalg.norm(A,np.inf) # L inf norm 最大行和"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 最小二乘解和伪逆"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 问题描述"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "所谓最小二乘问题的定义如下：\n",
    "\n",
    "假设 $y_i$ 与 $\\mathbf{x_i}$ 的关系可以用一组系数 $c_j$ 和对应的模型函数 $f_j(\\mathbf{x_i})$ 的模型表示：\n",
    "\n",
    "$$\n",
    "y_{i}=\\sum_{j}c_{j}f_{j}\\left(\\mathbf{x}_{i}\\right)+\\epsilon_{i}\n",
    "$$\n",
    "\n",
    "其中 $\\epsilon_i$ 表示数据的不确定性。最小二乘就是要优化这样一个关于 $c_j$ 的问题：\n",
    "$$\n",
    "J\\left(\\mathbf{c}\\right)=\\sum_{i}\\left|y_{i}-\\sum_{j}c_{j}f_{j}\\left(x_{i}\\right)\\right|^{2}\n",
    "$$\n",
    "\n",
    "其理论解满足：\n",
    "$$\n",
    "\\frac{\\partial J}{\\partial c_{n}^{*}}=0=\\sum_{i}\\left(y_{i}-\\sum_{j}c_{j}f_{j}\\left(x_{i}\\right)\\right)\\left(-f_{n}^{*}\\left(x_{i}\\right)\\right)\n",
    "$$\n",
    "\n",
    "改写为：\n",
    "$$\n",
    "\\begin{eqnarray*} \\sum_{j}c_{j}\\sum_{i}f_{j}\\left(x_{i}\\right)f_{n}^{*}\\left(x_{i}\\right) & = & \\sum_{i}y_{i}f_{n}^{*}\\left(x_{i}\\right)\\\\ \\mathbf{A}^{H}\\mathbf{Ac} & = & \\mathbf{A}^{H}\\mathbf{y}\\end{eqnarray*}\n",
    "$$\n",
    "\n",
    "其中：\n",
    "$$\n",
    "\\left\\{ \\mathbf{A}\\right\\} _{ij}=f_{j}\\left(x_{i}\\right).\n",
    "$$\n",
    "\n",
    "当 $\\mathbf{A^HA}$ 可逆时，我们有：\n",
    "$$\n",
    "\\mathbf{c}=\\left(\\mathbf{A}^{H}\\mathbf{A}\\right)^{-1}\\mathbf{A}^{H}\\mathbf{y}=\\mathbf{A}^{\\dagger}\\mathbf{y}\n",
    "$$\n",
    "\n",
    "矩阵 $\\mathbf{A}^{\\dagger}$ 叫做 $\\mathbf{A}$ 的伪逆。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 问题求解"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "注意到，我们的模型可以写为：\n",
    "$$\n",
    "\\mathbf{y}=\\mathbf{Ac}+\\boldsymbol{\\epsilon}.\n",
    "$$\n",
    "\n",
    "在给定 $\\mathbf{y}$ 和 $\\mathbf{A}$ 的情况下，我们可以使用 `linalg.lstsq` 求解 $\\mathbf c$。\n",
    "\n",
    "在给定 $\\mathbf{A}$ 的情况下，我们可以使用 `linalg.pinv` 或者 `linalg.pinv2` 求解 $\\mathbf{A}^{\\dagger}$。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####  例子"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "假设我们的数据满足：\n",
    "$$\n",
    "\\begin{align}\n",
    "y_{i} & =c_{1}e^{-x_{i}}+c_{2}x_{i} \\\\\n",
    "z_{i} & = y_i + \\epsilon_i\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "其中 $x_i = \\frac{i}{10},\\ i = 1,\\dots,10$，$c_1 = 5, c_2 = 2$，产生数据"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "c1, c2 = 5.0, 2.0\n",
    "i = np.r_[1:11]\n",
    "xi = 0.1*i\n",
    "yi = c1*np.exp(-xi) + c2*xi\n",
    "zi = yi + 0.05 * np.max(yi) * np.random.randn(len(yi))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "构造矩阵 $\\mathbf A$："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.90483742  0.1       ]\n",
      " [ 0.81873075  0.2       ]\n",
      " [ 0.74081822  0.3       ]\n",
      " [ 0.67032005  0.4       ]\n",
      " [ 0.60653066  0.5       ]\n",
      " [ 0.54881164  0.6       ]\n",
      " [ 0.4965853   0.7       ]\n",
      " [ 0.44932896  0.8       ]\n",
      " [ 0.40656966  0.9       ]\n",
      " [ 0.36787944  1.        ]]\n"
     ]
    }
   ],
   "source": [
    "A = np.c_[np.exp(-xi)[:, np.newaxis], xi[:, np.newaxis]]\n",
    "print A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "求解最小二乘问题："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 4.87016856  2.19081311]\n"
     ]
    }
   ],
   "source": [
    "c, resid, rank, sigma = linalg.lstsq(A, zi)\n",
    "\n",
    "print c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "其中 `c` 的形状与 `zi` 一致，为最小二乘解，`resid` 为 `zi - A c` 每一列差值的二范数，`rank` 为矩阵 `A` 的秩，`sigma` 为矩阵 `A` 的奇异值。\n",
    "\n",
    "查看拟合效果："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAEbCAYAAADDKt+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VPW9//HXJxuEsO+7IAWlgIoLgiim1lrBpa36KNcu\nLrWt16rXLvZqe6+K94qttd5bba11a6XWn9T+1Cp11xoraFGWIAhUVtnDmrAkkO1z/zgnOKSTZJJM\nZnKS9/PxyCMzZ77nnM+ZTN7zne9ZxtwdERGJnox0FyAiIk2jABcRiSgFuIhIRCnARUQiSgEuIhJR\nCnARkYhSgEuzmNlkM1tlZnvN7Atm9qKZXdbIZSwzsyktVWMC6/+RmT1cz+NXmNnbjVjeejM7K7z9\n4/qW3cg6q83s6CbMl29mG5NRg7QuWekuQJrOzNYDfYFKoApYDvweeMgTOMDfzIYBa4Esd69uYhn/\nBdzn7r8M7z8Xs/wrgKvc/YyYaY8BG939lppp7j62ietOCnf/Sc3tJD0nh597d7+zWcWlkJnNAEa4\n+9fTXYskRj3waHPgfHfvCgwFfgrcBDzayOVYM2oYSvDG0dY05zkRSQkFeBvh7vvcfQ4wHbjczMYA\nmNl5ZrbYzErMbIOZ3RYz29/C38Vmts/MTjWzEWb2VzPbaWY7zOwPZtYt3jrNbA1wNDAnHELJMbMC\nM7vKzI4FfgNMCpe9x8y+BXwF+Pdw2nPhcmKHHGaY2VNmNitc5jIzOylmnSeG27M3bPdHM/vvOur7\n2MxODG9/NRyCGB3ev8rMno1Z5+NxnpO9ZjaRsEdtZneb2W4zW2tm5ybyd4ldtpkNC2u4LKxth5n9\nOKbtBDN7N3yutpjZL80su47l9jKzOeHf9T0zuyPRYR4zu8nMNoXbt9LMzgq350fA9PBvszhse4WZ\nrQnbrjWzr4TTM83s5+E2rDGza8NtU6akkJ7sNsbd3wc2AaeHk/YDX3P3bsB5wDVm9oXwsZqhjW7u\n3sXd54f3ZwIDgNHAEGBGHesaAWwg/BTg7uUEYefuvhK4Gng3XHYPd38YeAK4K5xWU0ft4Z4LgCeB\nbsDzwK8AzCwHeBb4LdAjbPPFOPPXKADyw9tnAmvC3zX3C+LME/ucdHX3vxP0xk8FVgK9gJ+R+Kec\neLVNBkYBnwVuNbNjwumVwA3hOiaFj3+njuXeD+wD+gGXA5fVsa4jhOu6Fjg5/OR2DrDe3V8G7gRm\nh3+b8WaWB9wLnBu2nQQUhov6FsHr6QTgZOCSRNYvyaUAb5u2AD0B3P0td/8wvL0UmM0nIfZPwwTu\nvsbd33D3CnffCfxvTPvGqmsYoqHhibfd/eVwHP8PwPHh9IlAprv/0t2r3P1Z4L16lvMWn9R+OvCT\nmPtTwscTre1jd380rOn3wAAz69vAdtS1vNvd/ZC7fwAsIQhB3H2Ru7/n7tXu/jHwEHGeezPLBC4C\nbnP3g+6+AphVT+2xqoAOwBgzy3b3De6+NqbW2suoBsaZWa67F7l7zXDZl4H/dffN7r6HIPw17JRi\nCvC2aRCwGyAcFnnTzLabWTFBr7hXXTOaWT8zmx1+xC4BHq+vfQspirldCnQMP5oPBDbXaruRuoPj\nb8AZZtYfyAT+BEw2s6MIetiFdcwXz7aaG+5eGt7s3Ij54y6LYPvyAMxslJn9xcy2hs/9TOI/930I\nDkCIPbJkUyIrdvfVwHcJPlUVmdmTZjagjrYHCIbk/hXYEtZW82lhQK31b0hk/ZJcCvA2xsxOIQjw\nueGk/wf8GRjs7t0JxqVr/u7xPvLeSdBLGxsOu3ydpr9O4i2/OR+ztxJsW6yhdS0zDKtS4HrgLXff\nRxCe3wZix4u9jtup9gDBDuFPhc/9fxD/ud9BMNwyJGbakDjt4nL3J8Mjg44i2N67ah6K0/ZVdz8H\n6E8whFRzSORWgue+xtDa80rLU4BHnwGYWVczO59gXPjxmmETgl7iHncvN7MJBDsRa/5RdxB8RB4R\ns7zOwAFgr5kNAn7YjNqKgMG1dsQVEez4bIp3gSozu87MssKx/FMamOct4Do+GS4pqHUfjuzBx3tO\nUqUzwbh2abgT+Jp4jdy9CngGmGFmuWHbr5PYGPiocKdlB+AQcJDgDRuCN7dhZlbzmuprwbH9eUAF\nweuipu1TwL+Z2SAz6wHcnMj6JbkU4NE3x8z2EnyE/RFwD3BlzOPfAf4rbHML8MeaB8KhgJnAvPDo\nignA7cCJQAkwB3iapv9jvgF8CGwzs+3htEeBT4dHWjwTZx6Psz4P6y0nGPu9CtgDfBX4C1BeTw1v\nEQTj3+q4f8Q64zwnp9ZXUwJqz1vffDcSvMHuJRj/nl3PvNcR7OTdRjD+/SQxz0N49M6lcebtQLAv\nYAdBL7o3wesGgiEmgF1mtoAgH75HMGy1i2AHb82bysPAKwRj+AsIXicaA08xS+QLHSw4YWQvwbtv\nhbtPqPV4PsEJHDU7Q5529zuSWqlIHGY2H/i1u89Kdy3pZGZ3AX3d/coGG7fM+ofR/BOgpJESPRPT\ngXx3311Pm7fc/cIk1CRSJwtOuf8I2EnQAx8LvJzWotIg3JnYAVhKMIz0DYJPJtKONOZU+oY+Hunj\nk6TCMQTjr3kEx3Vf4u5F9c/SJnUhGDYZSLBf4efu/nx6S9IYeKolOoSylmBMtAp4MDwhI/bxMwl2\nqmwiGC+7MeZ4URERaQGJ9sAnu/tWM+sDvGZmK9099jCsRcAQdy81s6kEh62NSnaxIiLyiYR64EfM\nEFxLY7+731NPm3XASbFj5mamj1ciIk3g7nGHqBs8jNDMOplZl/B2HsG1E5bWatMv5tjRCQRvDP+0\nw9Pd29XPbbfdlvYatM3aZm1ztLe3PokMofQDng3zOQt4wt1fNbOrw1B+kOBCNteYWSXBmW//ksBy\nRUSkGRoMcHdfR3ixnVrTH4y5fT/B1dFERCRFdCZmC8rPz093CSmnbW4f2ts2t9btbfROzCavyMxT\ntS4RkbbCzPCm7sQUEZHWSQEuIhJRCnARkYhSgIuIRJQCXEQkohTgIiIRpQAXEYkoBbiISEQpwEVE\nIkoBLiISUQpwEZGIUoCLiESUAlxEJKIU4CIiEaUAFxGJKAW4iEhEKcBFRCJKAS4iElEKcBGRiFKA\ni4hElAJcRCSiFOAiIhGlABcRiSgFuIhIRCnARUQiSgEuIhJRCnARkYhSgIuIRJQCXEQkohTgIiIR\npQAXEYkoBbiISEQpwEVEIkoBLiISUQkFuJmtN7MPzGyxmb1XR5v7zGyVmS0xs/HJLVNERGrLSrCd\nA/nuvjveg2Y2DfiUu480s1OBB4CJSapRRETiaMwQitXz2IXALAB3nw90N7N+zSlMRETql2iAO/C6\nmS0ws2/FeXwQsDHm/iZgcHOLExGRuiU6hDLZ3beaWR/gNTNb6e5v12pTu4fuzS9PRETqklCAu/vW\n8PcOM3sWmADEBvhmYEjM/cHhtCPMmDHj8O38/Hzy8/MbXbCISFtWUFBAQUFBQm3Nvf6Ospl1AjLd\nfZ+Z5QGvAre7+6sxbaYB17n7NDObCPzC3SfWWo43tC4RETmSmeHucfdBJtID7wc8a2Y17Z9w91fN\n7GoAd3/Q3V80s2lmtho4AFyZpNpFRKQODfbAk7Yi9cBFRBqtvh64zsQUEYkoBbiISEQpwEVEIkoB\nLiISUQpwEZGIUoCLiESUAlxEJKIU4CIiEaUAFxGJKAW4iEhEKcBFRCJKAS4iElEKcBGRiFKAi4hE\nlAJcRCSiFOAiIhGlABcRiSgFuIhIRCnARUQiSgEuIhJRCnARkYhSgIuIRJQCXEQkohTgIiIRpQAX\nEYkoBbiISEQpwEVEIkoBLiISUW06wKu9mlmFs6ioqkh3KSIiSdemA3zvob08sfQJTn3kVBZtXZTu\nckREkqpNB3j3jt155WuvcMOpNzD1ianc/PrNlFaUprssEZGkaNMBDmBmXH7C5Xzwrx+woWQD4x4Y\nx2trXkt3WSIizWbunpoVmXmq1lWfl1a9xDUvXMPpQ0/nfz7/P/TN65vukkRE6mRmuLvFe6zN98Br\nmzpyKh9+50MGdhnI2F+P5aGFD1Ht1ekuq9leeAGKi4+cVlwcTBeRtqndBThAXk4eP/vcz3j9std5\nrPAxJv92MoXbCtNdVrNMngz/8R+fhHhxcXB/8uT01iUiLSehADezTDNbbGZz4jyWb2Yl4eOLzew/\nk19m4yTaGz2u33HM/cZcrhp/FZ//w+f57svfpeRgSeoKTaLu3WHmzCC0168Pfs+cGUwXkbYp0R74\nDcByoK5B7LfcfXz4c0dySmu6xvRGMyyDb574TT78zofsL9/P6PtH8/iSx2kN4/WN1b07/PCHMHx4\n8FvhLdK2NRjgZjYYmAY8AsQdSK9nelo0pTfau1NvHrnwEZ6d/iz3zr+XKY9NidywSnEx3H03rFsX\n/K79KURE2pYGj0Ixsz8BdwJdgRvd/YJaj58JPANsAjaHbZbHWU7Kj0JZvz7oja5bB8OGJT5fVXUV\njy5+lFvfvJUvHfsl7jjrDnp16tVSZSZFzaeMmjeq2vdFJJqafBSKmZ0PbHf3xdTdy14EDHH344Ff\nAn9uTrHJ0pzeaGZGJt8+6dusuHYF2ZnZjL5/NPfNv69Vn5I/b96RYV3zKWTevPTWJSItp94euJnd\nCXwdqAQ6EvTCn3b3y+qZZx1wkrvvrjXdb7vttsP38/Pzyc/Pb1bxdUl2b/TD7R/yvVe+x6a9m7jn\nnHuYOnJq8osWEQEKCgooKCg4fP/222+vswee8Ik84VBJvCGUfgS9dDezCcBT7j4szvwpG0J54YVg\nh2VsWBcXB73R885r2jLdnTkfzeHGV2/k6B5Hc8859zCm75jkFCwiUof6hlAaG+A/cPcLzexqAHd/\n0MyuBa4h6KWXAt9397/Hmb9VnInZXOVV5Tzw/gPMfHsmF42+iNvzb6df537pLktE2qikBHgSimgT\nAV5jT9keZr49k8cKH+OGU2/g+5O+T15OXrrLEpE2RqfSt4AeuT34+Tk/571vvcfyncsZ9atRPLTw\nISqrK9Ndmoi0E+qBJ8mCLQu46fWb2Lx3M3ecdQcXj74Ys1Z1eLyIRJCGUFLE3Xlt7Wvc/PrNZGZk\ncudZd3L20WcryEWkyRTgKVbt1fzpwz9xy5u3MKjrIGaeNZPThpyW7rJEJII0Bp5iGZbB9LHTWX7t\ncr427mtc+vSlTHtiGgu3LEx3aS1Cl7IVSQ8FeAvKysjiqhOv4qPrPuK8kedx4ewL+eLsL0buGisN\n0aVsRdJDQygpVFZRxkMLH+KueXcxcfBEbj3zVk7of0K6y0qKmtD+4Q+DSxfoGiwiyaEx8FamrKKM\n3yz4DXe/czenDDqFW6bcwskDT053Wc3W1IuHiUjdNAbeyuRm5/K9Sd9jzb+t4ezhZ/OlP36JqU9M\nZe6Guekurcl0KVuR1FOAJ1ljdujlZudy/anXs/r61Vx07EVc/ufLOfOxM3l59cuR+kKJ2IuFDRv2\nybXYFeIiLUtDKEnWnCshVlZXMnvZbO6adxdZGVncNPkmLvn0JWRlZKWm+CZqiYuHiUhAY+Ap1twd\neu7Oi6te5KfzfsrmvZv5waQfcOX4K+mU3anlihaRVkkBngbJ2qH37sZ3ufudu5m7YS5Xn3Q11024\nTlc/FGlHtBMzxZK5Q2/SkEk8M/0Z5n5jLjtLdzL6/tFc9dxVLC1amryCRSSS1ANPspb+bsqdpTt5\ncMGD3P/+/YzpO4YbTr2BaSOnkWF6LxZpizSEkkKp2qFXXlXOH5f9kXvn30vxwWKum3AdV5xwBd07\n6uwZkbZEAd6GuTt/3/R37nvvPl5e/TLTx0zn2lOuZVy/cekuTUSSQAHeTmzdt5WHFz3Mgwsf5Oge\nR3PNyddw8eiL6ZDVId2liUgTKcDbmYqqCuZ8NIcHFjzAkm1LuPz4y/n2Sd9mZK+R6S5NRBpJAd6O\nrd69mocXPsxjSx5jTJ8xfPPEb3LR6IvomNUx3aWJSAIU4EJ5VTnPrXyORxY/wsItC7l07KV8Y/w3\nGD9gfLpLE5F6KMDlCOuL1zOrcBa/K/wdPXJ7cMXxV/CVcV+hT16fdJcmIrUowCWuaq/mzXVvMmvJ\nLJ7/x/OcOexMLjvuMs4fdb52fIq0EgpwadC+Q/t4esXTPP7B4xRuK+Ti0Rfz1XFf5YyjztBJQiJp\npACXRtlYspEnlz3JE0ufYE/ZHqaPmc6l4y5lfP/xmMV9HYlIC1GAS5Mt276M2ctm8+SyJ8m0TKaP\nmc6Xx3yZsX3HKsxFUkABLs3m7ry/5X2e+vApnvrwKfJy8rhk9CVc8ulLOK7fcQrzdkDXfU8PBbgk\nVbVX897m93h6+dM8veJpMjMy+dKxX+Ki0RcxYdAEjZm3US19oTaJTwEuLcbdKdxWyDMrnuGZlc+w\np2wPF4y6gC8c+wXOGn6WThhqY5r7ZSXSeApwSZlVu1bx3D+e488r/8zS7Uv57PDPcsGoC5g2cpq+\niKKNSNaXlUhiFOCSFjtLd/LiqheZ89EcXl/7OiN7jmTayGlM/dRUTh54MpkZmekuURpJPfDUU4BL\n2pVXlTN3w1xeWvUSL61+iaIDRZx99Nl8fsTnOWfEOQzsMjDdJUoDNAaeHgpwaXU2lGzg1TWv8sqa\nV3hj7RsM6jqIs4efzedGfI4pR02hc07ndJcotaTrKJT2fvSLAlxatarqKhZuXchra17j9XWv8/7m\n9zm+//GcNewsPjP8M0waPInc7Nx0lylp0t57/gpwiZTSilLe2fgOf133VwrWF/BB0QecOOBEphw1\nhSlHTWHS4El06dAl3WVKCrXnsXcFuETa/vL9vLvxXd76+C3+9vHfWLR1Ecf0PoYzhp7BaUNOY/KQ\nyQzqOijdZUoLa69HvzQ7wM0sE1gAbHL3C+I8fh8wFSgFrnD3xXHaKMAlKQ5VHmLh1oW8/fHbvLPp\nHd7Z+A65WblMGjKJiYMmMnHwRE7of4KGXdoQ9cCbF+DfB04Curj7hbUemwZc5+7TzOxU4F53nxhn\nGQpwaRHuzqrdq5i/aT7vbnqX+Zvns2LHCkb3Gc0pA0/h5IEnc/LAkxnTZwzZmdnpLlcaSWPgzQhw\nMxsMPAbMBL5fuwduZr8B3nT3P4b3VwJnuntRrXYKcEmZsooyCrcVsmDLAt7f8j4LtixgffF6xvQd\nw4n9T2T8gPGc0P8ExvUdR15OXrrLlXroKJTmBfifgDuBrsCNcQJ8DvATd38nvP86cJO7L6zVTgEu\nabW/fD9Lti1h0dZFFG4rZPG2xazcuZIh3YZwXL/jOK7vcYztO5Zx/cYxvPtwnWgkrUJ9AZ7VwIzn\nA9vdfbGZ5dfXtNZ9JbW0Op1zOjN56GQmD518eFpFVQX/2PUPlmxbwrLty/ht4W9ZWrSU7Qe2c0zv\nYxjTZwyje49mdJ/RHNv7WEb0GKFvK5JWo94AB04DLgzHuTsCXc3s9+5+WUybzcCQmPuDw2n/ZMaM\nGYdv5+fnk5+f34SSRZInOzObsX3HMrbv2COm7y/fz4odK1i+Yzkrdq7gscLHWLlzJRtKNjCo6yBG\n9RrFyJ4jGdVrFCN6jGBEzxEM6z6MnMycNG2JtCR3T9klkwsKCigoKEiobcKHEZrZmcQfQondiTkR\n+IV2YkpbVVFVwdo9a1m1exWrdq3io10fsWbPGtbsWcOmvZvo37k/R/c4muHdh3NUt6MY1n0YR3U/\niqHdhjK462AFfCtUWV1J0f4ituzbwuZ9m9m0d9Phnw0lG9i4dyPZGdl8dP1HaakvKceBhwH+A3e/\n0MyuBnD3B8PHfgWcCxwArnT3RXHmV4BLm1ZRVcHGvRtZu2ct6/as4+OSj1lfvJ4NJRvYULKBLfu2\n0DO3J0O6DWFQl0HBT9dBDOg8gAFdBjCg8wD6de5Hn059NP7eTNVeTfHBYrYf2E7R/qLg94EiivYX\nsW3/Nrbu3xr87NvKjtId9MrtxaCugxjYZSCDuww+/Dca2m0oQ7oNYXDXwWm7NLJO5BFpBaqqqyg6\nUHS4d7d572Y279vMln1bDofKtv3b2FO2h565Pemb15c+eX3o06kPvTv1pnen3vTK7UWvTr3omduT\nHh170L1jd7p37E63jt3IzcptU9+M5O4cqjrE3kN7KTlYQvHB4sM/ew7uYU/ZHnaX7WZ32W52le1i\nV9kudpbuZGfpTnaX7aZzTmf65vWlb15f+uX1C34692NA5wH079yf/p37M7DLQPrm9W3Vh5cqwEUi\npLK6kh0HdrCjdMfh3ztLd7KrdBc7SnccEV4lh0rYU7aHkkMlVFZX0rVDV7rkdKFLhy50yelC55zO\ndM7pTF5OHp2yOpGbnUtuVi4dszqSm51Lh8wO5GTmkJOZQ3ZmNtkZ2WRnZpOVkUVWRhaZlkmGZZBh\nGYffHAzDcdwdx6mqrqLaq6nyKiqrK6msrqSiqoKK6grKq8opryrnUOUhDlUd4mDlQcoqyiirLDv8\n+0DFAQ6UH2B/+f7DP/vK97Hv0D4AunXsRrcO3ejWsRvdO3Y//MbVM7fn4Tey3p1606tTL3rl9qJP\nXh965fZq1aHcGApwkXagvKqckoMlh8NvX/m+w8F4oOIAZRVBWB6sPHg4SMuryjlUdYjyqnIqqisO\nB29VdRUV1RVUe3UQztVVAIeD28wwDDMj0zLJzAiCPjsjm8yMzMNvBNkZ2XTI7ECHrA50yOxAbnbw\n5tExqyOdsjuRm5VLXk4eedl55OXkHfGm07VDVx3xgwJcRKTRWssJRPUFuL59VkQkjsmTg1P2i4uD\n+zWn8E+eXP98qaQeuIhIHVrDRbQ0hCIi0kTpvoythlCkTXrhhU8+3tYoLg6miyRDcXHQ8163Lvhd\n+/WWbgpwiawojFFKdMVetnbYsOB37OutNdAQikRaaxijlLYpCkehKMAl8tI9RinSkjQGLm1Wax+j\nFGlJCnCJrHSOUaZrB6p23EosBbhE1rx5R455d+8e3J83r+XXna4dqNpxK7E0Bi7SROnagaodt+2L\ndmKKtJB07UDVjtv2QzsxRVpAunagaset1FCAizRBunagRuHkEkkdDaGINEG6TvJoLSeXSOpoDFxE\nJKI0Bi4i0gYpwEVEIkoBLiISUQpwEZGIUoCLiESUAlxEJKIU4CIiEaUAFxGJKAW4iEhEKcBFRCJK\nAS4iElEKcBGRiFKAi4hElAJcRCSiFOAiIhGlABcRiSgFuIhIRDUY4GbW0czmm1mhmS03s5/EaZNv\nZiVmtjj8+c+WKVdERGpkNdTA3Q+a2WfcvdTMsoC5Zna6u8+t1fQtd7+wZcoUEZHaEhpCcffS8GYO\nkAnsjtMs7ne2iYhIy0gowM0sw8wKgSLgTXdfXquJA6eZ2RIze9HMPp3sQkVE5EiJ9sCr3f0EYDAw\nxczyazVZBAxx9+OBXwJ/TmqVIiLyTxocA4/l7iVm9gJwMlAQM31fzO2XzOzXZtbT3Y8YapkxY8bh\n2/n5+eTn5zetahGRNqqgoICCgoKE2pq719/ArDdQ6e7FZpYLvALc7u5vxLTpB2x3dzezCcBT7j6s\n1nK8oXWJiMiRzAx3j7uPMZEe+ABglpllEAy5PO7ub5jZ1QDu/iBwCXCNmVUCpcC/JKd0ERGpS4M9\n8KStSD1wEZFGq68HrjMxRUQiSgEuIhJRCnARkYhSgIuIRJQCXEQkohTgIiIRpQAXEYkoBbiISEQp\nwEVEIkoBLiISUQpwEZGIUoCLiESUAlxEJKIU4CIiEaUAFxGJKAW4iEhEKcBFRCJKAS4iElEKcBGR\niFKAi4hElAJcRCSiFOAiIhGlABcRiSgFuIhIRCnARUQiSgEuIhJRCnARkYhSgIuIRJQCXEQkohTg\nIiIRpQAXEYkoBbiISEQpwEVEIkoBLiISUQpwEZGIqjfAzayjmc03s0IzW25mP6mj3X1mtsrMlpjZ\n+JYpVUREYtUb4O5+EPiMu58AHAd8xsxOj21jZtOAT7n7SODbwAMtVWzUFBQUpLuElNM2tw/tbZtb\n6/Y2OITi7qXhzRwgE9hdq8mFwKyw7Xygu5n1S2aRUdVa/+gtSdvcPrS3bW6t29tggJtZhpkVAkXA\nm+6+vFaTQcDGmPubgMHJK1FEROJJpAdeHQ6hDAammFl+nGZWe7Yk1CYiIvUw98Sz1sxuAcrc/ecx\n034DFLj77PD+SuBMdy+qNa9CXUSkCdy9dicZgKz6ZjKz3kCluxebWS7wOeD2Ws2eB64DZpvZRKC4\ndnjXV4CIiDRNvQEODABmmVkGwXDL4+7+hpldDeDuD7r7i2Y2zcxWAweAK1u2ZBERgUYOoYiISOuR\n9DMxzexcM1sZnthzUx1t2tSJPw1ts5l9NdzWD8xsnpkdl446kyWRv3HY7hQzqzSzi1JZX0tI8HWd\nb2aLzWyZmRWkuMSkS+B13dvMXg5P9FtmZlekocykMbPfmlmRmS2tp03ryi53T9oPwXHiq4FhQDZQ\nCIyu1WYa8GJ4+1Tg78msIdU/CW7zJKBbePvcKG9zItsb0+6vwF+Ai9Nddwr+xt2BD4HB4f3e6a47\nBds8A/hJzfYCu4CsdNfejG0+AxgPLK3j8VaXXcnugU8AVrv7enevAGYDX6jVpq2d+NPgNrv7u+5e\nEt6dT7SPk0/kbwxwPfD/gR2pLK6FJLLNXwGedvdNAO6+M8U1Jlsi27wV6Bre7grscvfKFNaYVO7+\nNrCnniatLruSHeDxTuoZlECbKAdaItsc6yrgxRatqGU1uL1mNojgn73msgpR39GSyN94JNDTzN40\nswVm9vWUVdcyEtnmh4ExZrYFWALckKLa0qXVZVdDR6E0VqL/qG3pxJ+EazezzwDfACa3XDktLpHt\n/QVws7u7mRn//PeOmkS2ORs4Efgs0Al418z+7u6rWrSylpPINv8YKHT3fDMbAbxmZse7+74Wri2d\nWlV2JTvANwNDYu4PIXiXqq/N4HBaVCWyzYQ7Lh8GznX3+j6mtXaJbO9JBOcFQDA2OtXMKtz9+dSU\nmHSJbPNGYKe7lwFlZvY34HggqgGeyDafBswEcPc1ZrYOOAZYkJIKU6/VZVeyh1AWACPNbJiZ5QDT\nCU70ifU+xBxCAAABwklEQVQ8cBlAfSf+REiD22xmQ4FngK+5++o01JhMDW6vux/t7sPdfTjBOPg1\nEQ5vSOx1/Rxwupllmlkngp1cta8bFCWJbPNK4GyAcCz4GGBtSqtMrVaXXUntgbt7pZldB7xCsBf7\nUXdf0ZZP/Elkm4FbgR7AA2GvtMLdJ6Sr5uZIcHvblARf1yvN7GXgA6AaeNj/+cJvkZHg3/lO4Hdm\ntoSgM/jv7l77aqWRYWZPAmcCvc1sI3AbwdBYq80uncgjIhJR+ko1EZGIUoCLiESUAlxEJKIU4CIi\nEaUAFxGJKAW4iEhEKcBFRCJKAS4iElHJvhaKSCSYWSbB6eFHE1zHZAJwj7u35VPBpY1RD1zaq+OB\npwmu3ZEB/Ing+tYikaEAl3bJ3Re5+yGCb0sqcPcCdy8zs8+nuzaRRCnApV0Kv6+zNzDW3deZ2ekA\n7v5KmksTSZguZiXtkpndAhQBQ4GFwHbgEHC6u/8inbWJJEo7MaVdcvf/rj0t/FaZ4jSUI9IkGkIR\n+cSJKMAlQjSEIiISUeqBi4hElAJcRCSiFOAiIhGlABcRiSgFuIhIRCnARUQiSgEuIhJRCnARkYhS\ngIuIRNT/AXU6OlSJKzNiAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7ff0d40c4f90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "xi2 = np.r_[0.1:1.0:100j]\n",
    "yi2 = c[0]*np.exp(-xi2) + c[1]*xi2\n",
    "\n",
    "plt.plot(xi,zi,'x',xi2,yi2)\n",
    "plt.axis([0,1.1,3.0,5.5])\n",
    "plt.xlabel('$x_i$')\n",
    "plt.title('Data fitting with linalg.lstsq')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 广义逆"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`linalg.pinv` 或 `linalg.pinv2` 可以用来求广义逆，其区别在于前者使用求最小二乘解的算法，后者使用求奇异值的算法求解。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 矩阵分解"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 特征值和特征向量"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 问题描述"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对于给定的 $N \\times N$ 矩阵 $\\mathbf A$，特征值和特征向量问题相当与寻找标量 $\\lambda$ 和对应的向量 $\\mathbf v$ 使得：\n",
    "$$\n",
    "\\mathbf{Av} = \\lambda \\mathbf{v}\n",
    "$$\n",
    "\n",
    "矩阵的 $N$ 个特征值（可能相同）可以通过计算特征方程的根得到：\n",
    "$$\n",
    "\\left|\\mathbf{A} - \\lambda \\mathbf{I}\\right| = 0\n",
    "$$\n",
    "\n",
    "然后利用这些特征值求（归一化的）特征向量。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 问题求解"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- `linalg.eig(A)` \n",
    "    - 返回矩阵的特征值与特征向量\n",
    "- `linalg.eigvals(A)`\n",
    "    - 返回矩阵的特征值\n",
    "- `linalg.eig(A, B)`\n",
    "    - 求解 $\\mathbf{Av} = \\lambda\\mathbf{Bv}$ 的问题"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 例子"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "矩阵为\n",
    "$$\n",
    "\\begin{split}\\mathbf{A}=\\left[\\begin{array}{ccc} 1 & 5 & 2\\\\ 2 & 4 & 1\\\\ 3 & 6 & 2\\end{array}\\right].\\end{split}\n",
    "$$\n",
    "\n",
    "特征多项式为：\n",
    "$$\n",
    "\\begin{eqnarray*} \\left|\\mathbf{A}-\\lambda\\mathbf{I}\\right| & = & \\left(1-\\lambda\\right)\\left[\\left(4-\\lambda\\right)\\left(2-\\lambda\\right)-6\\right]-\\\\  &  & 5\\left[2\\left(2-\\lambda\\right)-3\\right]+2\\left[12-3\\left(4-\\lambda\\right)\\right]\\\\  & = & -\\lambda^{3}+7\\lambda^{2}+8\\lambda-3.\\end{eqnarray*}\n",
    "$$\n",
    "\n",
    "特征根为：\n",
    "$$\n",
    "\\begin{eqnarray*} \\lambda_{1} & = & 7.9579\\\\ \\lambda_{2} & = & -1.2577\\\\ \\lambda_{3} & = & 0.2997.\\end{eqnarray*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 7.95791620+0.j -1.25766471+0.j  0.29974850+0.j]\n",
      "[ 1.  1.  1.]\n",
      "3.23301824835e-15\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[1, 5, 2], \n",
    "              [2, 4, 1],\n",
    "              [3, 6, 2]])\n",
    "\n",
    "la, v = linalg.eig(A)\n",
    "\n",
    "print la\n",
    "\n",
    "\n",
    "# 验证是否归一化\n",
    "print np.sum(abs(v**2),axis=0)\n",
    "\n",
    "# 第一个特征值\n",
    "l1 = la[0]\n",
    "# 对应的特征向量\n",
    "v1 = v[:, 0].T\n",
    "\n",
    "# 验证是否为特征值和特征向量对\n",
    "print linalg.norm(A.dot(v1)-l1*v1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 奇异值分解"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 问题描述"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$M \\times N$ 矩阵 $\\mathbf A$ 的奇异值分解为：\n",
    "$$\n",
    "\\mathbf{A=U}\\boldsymbol{\\Sigma}\\mathbf{V}^{H}\n",
    "$$\n",
    "\n",
    "其中 $\\boldsymbol{\\Sigma}, (M \\times N)$ 只有对角线上的元素不为 0，$\\mathbf U, (M \\times M)$ 和 $\\mathbf V, (N \\times N)$ 为正交矩阵。\n",
    "\n",
    "其具体原理可以查看维基百科：\n",
    "https://en.wikipedia.org/wiki/Singular_value_decomposition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 问题求解"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- `U,s,Vh = linalg.svd(A)` \n",
    "    - 返回 $U$ 矩阵，奇异值 $s$，$V^H$ 矩阵\n",
    "- `Sig = linalg.diagsvd(s,M,N)`\n",
    "    - 从奇异值恢复 $\\boldsymbol{\\Sigma}$ 矩阵"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 例子"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "奇异值分解："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "A = np.array([[1,2,3],[4,5,6]])\n",
    "\n",
    "U, s, Vh = linalg.svd(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\boldsymbol{\\Sigma}$ 矩阵："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 9.508032    0.          0.        ]\n",
      " [ 0.          0.77286964  0.        ]]\n"
     ]
    }
   ],
   "source": [
    "M, N = A.shape\n",
    "Sig = linalg.diagsvd(s,M,N)\n",
    "\n",
    "print Sig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "检查正确性："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1 2 3]\n",
      " [4 5 6]]\n",
      "[[ 1.  2.  3.]\n",
      " [ 4.  5.  6.]]\n"
     ]
    }
   ],
   "source": [
    "print A\n",
    "print U.dot(Sig.dot(Vh))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### LU 分解"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$M \\times N$ 矩阵 $\\mathbf A$ 的 `LU` 分解为：\n",
    "$$\n",
    "\\mathbf{A}=\\mathbf{P}\\,\\mathbf{L}\\,\\mathbf{U}\n",
    "$$\n",
    "\n",
    "$\\mathbf P$ 是 $M \\times M$ 的单位矩阵的一个排列，$\\mathbf L$ 是下三角阵，$\\mathbf U$ 是上三角阵。 \n",
    "\n",
    "可以使用 `linalg.lu` 进行 LU 分解的求解：\n",
    "\n",
    "具体原理可以查看维基百科：\n",
    "https://en.wikipedia.org/wiki/LU_decomposition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.  1.]\n",
      " [ 1.  0.]]\n",
      "[[ 1.    0.  ]\n",
      " [ 0.25  1.  ]]\n",
      "[[ 4.    5.    6.  ]\n",
      " [ 0.    0.75  1.5 ]]\n",
      "[[ 1.  2.  3.]\n",
      " [ 4.  5.  6.]]\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[1,2,3],[4,5,6]])\n",
    "\n",
    "P, L, U = linalg.lu(A)\n",
    "\n",
    "print P\n",
    "print L\n",
    "print U\n",
    "\n",
    "print P.dot(L).dot(U)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cholesky 分解"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Cholesky` 分解是一种特殊的 `LU` 分解，此时要求 $\\mathbf A$ 为 Hermitian 正定矩阵 （$\\mathbf A = \\mathbf{A^H}$）。\n",
    "\n",
    "此时有：\n",
    "$$\n",
    "\\begin{eqnarray*} \\mathbf{A} & = & \\mathbf{U}^{H}\\mathbf{U}\\\\ \\mathbf{A} & = & \\mathbf{L}\\mathbf{L}^{H}\\end{eqnarray*}\n",
    "$$\n",
    "即\n",
    "$$\n",
    "\\mathbf{L}=\\mathbf{U}^{H}.\n",
    "$$\n",
    "\n",
    "可以用 `linalg.cholesky` 求解。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### QR 分解"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$M×N$ 矩阵 $\\mathbf A$ 的 `QR` 分解为：\n",
    "$$\n",
    "\\mathbf{A=QR}\n",
    "$$\n",
    "\n",
    "$\\mathbf R$ 为上三角形矩阵，$\\mathbf Q$ 是正交矩阵。\n",
    "\n",
    "维基链接：\n",
    "https://en.wikipedia.org/wiki/QR_decomposition\n",
    "\n",
    "可以用 `linalg.qr` 求解。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Schur 分解"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对于 $N\\times N$ 方阵 $\\mathbf A$, `Schur` 分解要求找到满足下式的矩阵：\n",
    "$$\n",
    "\\mathbf{A=ZTZ^H}\n",
    "$$\n",
    "\n",
    "其中 $\\mathbf Z$ 是正交矩阵，$\\mathbf T$ 是一个上三角矩阵。\n",
    "\n",
    "维基链接：\n",
    "https://en.wikipedia.org/wiki/Schur_decomposition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1 3 2]\n",
      " [1 4 5]\n",
      " [2 3 6]]\n",
      "[[ 9.90012467  1.78947961 -0.65498528]\n",
      " [ 0.          0.54993766 -1.57754789]\n",
      " [ 0.          0.51260928  0.54993766]] [[ 0.36702395 -0.85002495 -0.37782404]\n",
      " [ 0.63681656 -0.06646488  0.76814522]\n",
      " [ 0.67805463  0.52253231 -0.51691576]]\n",
      "[[ 1.  3.  2.]\n",
      " [ 1.  4.  5.]\n",
      " [ 2.  3.  6.]]\n"
     ]
    }
   ],
   "source": [
    "A = np.mat('[1 3 2; 1 4 5; 2 3 6]')\n",
    "\n",
    "print A\n",
    "\n",
    "T, Z = linalg.schur(A)\n",
    "\n",
    "print T, Z\n",
    "\n",
    "print Z.dot(T).dot(Z.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 矩阵函数"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "考虑函数 $f(x)$ 的泰勒展开：\n",
    "$$\n",
    "f\\left(x\\right)=\\sum_{k=0}^{\\infty}\\frac{f^{\\left(k\\right)}\\left(0\\right)}{k!}x^{k}\n",
    "$$\n",
    "\n",
    "对于方阵，矩阵函数可以定义如下：\n",
    "$$\n",
    "f\\left(\\mathbf{A}\\right)=\\sum_{k=0}^{\\infty}\\frac{f^{\\left(k\\right)}\\left(0\\right)}{k!}\\mathbf{A}^{k}\n",
    "$$\n",
    "\n",
    "这也是计算矩阵函数的最好的方式。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 指数和对数函数"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 指数"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "指数可以定义如下：\n",
    "$$\n",
    "e^{\\mathbf{A}}=\\sum_{k=0}^{\\infty}\\frac{1}{k!}\\mathbf{A}^{k}\n",
    "$$\n",
    "\n",
    "`linalg.expm3` 使用的是泰勒展开的方法计算结果："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  51.96890355   74.73648784]\n",
      " [ 112.10473176  164.07363531]]\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[1, 2], [3, 4]])\n",
    "\n",
    "print linalg.expm3(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "另一种方法先计算 A 的特征值分解：\n",
    "$$\n",
    "\\mathbf{A}=\\mathbf{V}\\boldsymbol{\\Lambda}\\mathbf{V}^{-1}\n",
    "$$\n",
    "\n",
    "然后有（正交矩阵和对角阵的性质）：\n",
    "$$\n",
    "e^{\\mathbf{A}}=\\mathbf{V}e^{\\boldsymbol{\\Lambda}}\\mathbf{V}^{-1}\n",
    "$$\n",
    "\n",
    "`linalg.expm2` 使用的就是这种方法："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  51.9689562    74.73656457]\n",
      " [ 112.10484685  164.07380305]]\n"
     ]
    }
   ],
   "source": [
    "print linalg.expm2(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "最优的方法是用 [`Padé` 近似](https://en.wikipedia.org/wiki/Pad%C3%A9_approximant) 实现，`Padé` 近似往往比截断的泰勒级数准确，而且当泰勒级数不收敛时，`Padé` 近似往往仍可行，所以多用于在计算机数学中。\n",
    "\n",
    "`linalg.expm` 使用的就是这种方法："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  51.9689562    74.73656457]\n",
      " [ 112.10484685  164.07380305]]\n"
     ]
    }
   ],
   "source": [
    "print linalg.expm(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 对数"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "指数的逆运算，可以用 `linalg.logm` 实现："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1 2]\n",
      " [3 4]]\n",
      "[[ 1.  2.]\n",
      " [ 3.  4.]]\n"
     ]
    }
   ],
   "source": [
    "print A\n",
    "print linalg.logm(linalg.expm(A))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 三角函数"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "根据欧拉公式，其定义为：\n",
    "$$\n",
    "\\begin{eqnarray*} \\sin\\left(\\mathbf{A}\\right) & = & \\frac{e^{j\\mathbf{A}}-e^{-j\\mathbf{A}}}{2j}\\\\ \\cos\\left(\\mathbf{A}\\right) & = & \\frac{e^{j\\mathbf{A}}+e^{-j\\mathbf{A}}}{2}.\\end{eqnarray*}\n",
    "$$\n",
    "\n",
    "正切函数定义为：\n",
    "$$\n",
    "\\tan\\left(x\\right)=\\frac{\\sin\\left(x\\right)}{\\cos\\left(x\\right)}=\\left[\\cos\\left(x\\right)\\right]^{-1}\\sin\\left(x\\right)\n",
    "$$\n",
    "\n",
    "因此矩阵的正切函数定义为：\n",
    "$$\n",
    "\\left[\\cos\\left(\\mathbf{A}\\right)\\right]^{-1}\\sin\\left(\\mathbf{A}\\right).\n",
    "$$\n",
    "\n",
    "具体实现：\n",
    "- `linalg.sinm`\n",
    "- `linalg.cosm`\n",
    "- `linalg.tanm`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 双曲三角函数"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{eqnarray*} \\sinh\\left(\\mathbf{A}\\right) & = & \\frac{e^{\\mathbf{A}}-e^{-\\mathbf{A}}}{2}\\\\ \\cosh\\left(\\mathbf{A}\\right) & = & \\frac{e^{\\mathbf{A}}+e^{-\\mathbf{A}}}{2}\\\\ \\tanh\\left(\\mathbf{A}\\right) & = & \\left[\\cosh\\left(\\mathbf{A}\\right)\\right]^{-1}\\sinh\\left(\\mathbf{A}\\right).\\end{eqnarray*}\n",
    "\n",
    "具体实现：\n",
    "- `linalg.sinhm`\n",
    "- `linalg.coshm`\n",
    "- `linalg.tanhm`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 特殊矩阵"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Scipy` 提供了一些特殊矩阵的实现，具体可以参考：\n",
    "\n",
    "http://docs.scipy.org/doc/scipy/reference/tutorial/linalg.html#special-matrices"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}