{ "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": [ "" ] }, "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 }