SciPy で 線形代数を扱う

SciPyTutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SciPy Tutorial\n",
    "https://showa-yojyo.github.io/notebook/python-scipy/tutorial.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 線形方程式を解く\n",
    "https://showa-yojyo.github.io/notebook/python-scipy/linear-equations.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import linalg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[2, -1, 2], [1, -1, -2], [-2, 1, -1]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "b = np.array([8, -1, -6])[:, np.newaxis]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 2, -1,  2],\n",
       "       [ 1, -1, -2],\n",
       "       [-2,  1, -1]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 8],\n",
       "       [-1],\n",
       "       [-6]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.],\n",
       "       [-2.],\n",
       "       [ 2.]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "linalg.solve(A, b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### LU分解"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "A = np.array([[1, 2,2], [2, 5, 6], [3, 8, 12]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1,  2,  2],\n",
       "       [ 2,  5,  6],\n",
       "       [ 3,  8, 12]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "P, L, U = linalg.lu(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.,  1.,  0.],\n",
       "       [ 0.,  0.,  1.],\n",
       "       [ 1.,  0.,  0.]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.        ,  0.        ,  0.        ],\n",
       "       [ 0.33333333,  1.        ,  0.        ],\n",
       "       [ 0.66666667,  0.5       ,  1.        ]])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  3.        ,   8.        ,  12.        ],\n",
       "       [  0.        ,  -0.66666667,  -2.        ],\n",
       "       [  0.        ,   0.        ,  -1.        ]])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  1.,   2.,   2.],\n",
       "       [  2.,   5.,   6.],\n",
       "       [  3.,   8.,  12.]])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P @ L @ U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  1.,   2.,   2.],\n",
       "       [  2.,   5.,   6.],\n",
       "       [  3.,   8.,  12.]])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.dot(P, np.dot(L, U))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  1.,   2.,   2.],\n",
       "       [  2.,   5.,   6.],\n",
       "       [  3.,   8.,  12.]])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.dot(np.dot(P, L), U)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "PLU を手で計算してみよう"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 最小二乗法\n",
    "https://showa-yojyo.github.io/notebook/python-scipy/least-squares.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "xd = np.array([72, 67, 65, 55, 25, 36, 56, 34,\n",
    "               18, 71, 67, 48, 72, 51, 53])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "yd = np.array([202, 186, 187, 180, 156, 169, 174,\n",
    "               172, 153, 199, 193, 174, 198, 183, 178])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[72],\n",
       "       [67],\n",
       "       [65],\n",
       "       [55],\n",
       "       [25],\n",
       "       [36],\n",
       "       [56],\n",
       "       [34],\n",
       "       [18],\n",
       "       [71],\n",
       "       [67],\n",
       "       [48],\n",
       "       [72],\n",
       "       [51],\n",
       "       [53]])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xd[:, np.newaxis]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(15,)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xd.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "A = np.c_[xd[:, np.newaxis], np.ones(xd.shape[0])]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 72.,   1.],\n",
       "       [ 67.,   1.],\n",
       "       [ 65.,   1.],\n",
       "       [ 55.,   1.],\n",
       "       [ 25.,   1.],\n",
       "       [ 36.,   1.],\n",
       "       [ 56.,   1.],\n",
       "       [ 34.,   1.],\n",
       "       [ 18.,   1.],\n",
       "       [ 71.,   1.],\n",
       "       [ 67.,   1.],\n",
       "       [ 48.,   1.],\n",
       "       [ 72.,   1.],\n",
       "       [ 51.,   1.],\n",
       "       [ 53.,   1.]])"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "B = yd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([202, 186, 187, 180, 156, 169, 174, 172, 153, 199, 193, 174, 198,\n",
       "       183, 178])"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "X, residues, rank, s = linalg.lstsq(A, B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([   0.79772656,  138.25306758]),\n",
       " array([ 272.43120523]),\n",
       " 2,\n",
       " array([ 214.24658905,    1.18282701]))"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np. linalg.lstsq(A, B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([   0.79772656,  138.25306758]),\n",
       " 272.43120523201486,\n",
       " 2,\n",
       " array([ 214.24658905,    1.18282701]))"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "linalg.lstsq(A, B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "a = X[0]\n",
    "b = X[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Line: y = 0.798x +138.253\n"
     ]
    }
   ],
   "source": [
    "print(\"Line: y = {:.3f}x {:+.3f}\".format(a, b))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd81dX9x/HXhyFhI1MkYgDZK0DCCiIORNRqXai1Ik5a\n0YrWXRFEqCgUBwoKioAVkapF2vIrqFXZI2CQPQ0SZCN7Jjm/P74XDEhISHLzveP9fDx4cO+533zz\nOSS8882553uOOecQEZHIVcTvAkREJLgU9CIiEU5BLyIS4RT0IiIRTkEvIhLhFPQiIhFOQS8iEuEU\n9CIiEU5BLyIS4Yr5XQBA5cqVXVxcnN9liIiElYULF+5wzlXJ6biQCPq4uDiSk5P9LkNEJKyY2Ybc\nHKehGxGRCKegFxGJcAp6EZEIFxJj9Kdz7Ngx0tLSOHz4sN+lSAGLiYkhNjaW4sWL+12KSFQI2aBP\nS0ujbNmyxMXFYWZ+lyMFxDnHzp07SUtLo1atWn6XIxIVQnbo5vDhw1SqVEkhH2HMjEqVKuk3NZFC\nFLJBDyjkI5S+riKFK6SDXkRE8i9kx+jl9J5//nk6duzIFVdc4XcpInI2Zi6CjMxftxctAh1aBvVT\nK+hzwTmHc44iRfL+C1B6ejrFiuX/n7t///75Pkd+FVRfRKLK6UL+TO0FSEM32UhNTaV+/fp0796d\nJk2asHHjRqZNm0a7du1o2bIlt9xyC/v37wdgypQpNGjQgFatWvGnP/2Ja6+9FoB+/fpx5513kpSU\nxJ133klGRgZPPPEEiYmJNGvWjHfeeQeAzZs307FjR+Lj42nSpAkzZswgIyODHj160KRJE5o2bcqr\nr74KQI8ePfjkk08A+Oqrr2jRogVNmzblnnvu4ciRI4C3pETfvn1p2bIlTZs2ZeXKlb/qX0ZGBo8/\n/jhNmjShWbNmDBs27MTH7tixA4Dk5GQ6dep02r60bduWZcuWnThfp06dSE5O5sCBA9xzzz20bt2a\nFi1a8Pnnnxf0l0ZEzlJ4XJb9ZQYs3V6w52xSBQZefMZD1qxZw9ixY2nbti07duxgwIABfPnll5Qu\nXZqXX36ZoUOH8uSTT9KzZ0+mT59OrVq1uP322086x/Lly5k5cyYlS5Zk5MiRlC9fngULFnDkyBGS\nkpK48sor+eyzz+jSpQt/+ctfyMjI4ODBg6SkpLBp0yaWLl0KwO7du0867+HDh+nRowdfffUV9erV\no3v37owYMYLevXsDULlyZRYtWsTw4cMZMmQI77777kkfP3LkSFJTU0lJSaFYsWLs2rUrx3+yrH15\n9dVXmThxIi+88AKbN29m8+bNJCQk8Oyzz3LZZZcxevRodu/eTevWrbniiisoXbp0jucXkeDQFf0Z\nXHjhhbRt2xaAuXPnsnz5cpKSkoiPj2fs2LFs2LCBlStXUrt27RNzwk8N+uuuu46SJUsCMG3aNMaN\nG0d8fDxt2rRh586drFmzhsTERN5//3369evHkiVLKFu2LLVr12b9+vU8/PDD/Pe//6VcuXInnXfV\nqlXUqlWLevXqAXDXXXcxffr0E6/feOONALRq1YrU1NRf9e3LL7+kZ8+eJ4ZgKlasmOO/R9a+dOvW\n7cRvFhMnTuTmm28+0cdBgwYRHx9Pp06dOHz4MD/++GOO5xaR4AmPK/ocrryDJetVqHOOzp0789FH\nH510TEpKylmdY9iwYXTp0uVXx02fPp3//Oc/9OjRg8cee4zu3buzePFipk6dyttvv83EiRMZPXp0\nrmsvUaIEAEWLFiU9PT3XH1esWDEyM70xw1PnumftS40aNahUqRLff/89H3/8MW+//faJPn766afU\nr18/159TRIJLV/S51LZtW2bNmsXatWsBOHDgAKtXr6Z+/fqsX7/+xFXzxx9/nO05unTpwogRIzh2\n7BgAq1ev5sCBA2zYsIFq1apx//33c99997Fo0SJ27NhBZmYmN910EwMGDGDRokUnnat+/fqkpqae\nqOeDDz7gkksuyXV/OnfuzDvvvHPih8DxoZu4uDgWLlwIwKeffnrGc9x666288sor7Nmzh2bNmp3o\n47Bhw3DOAfDdd9/luiaRiFY0m7jNrr0A5XhFb2YXAOOAaoADRjrnXjezisDHQByQCnRzzv1s3t0w\nrwNXAweBHs65Rac7dzipUqUKY8aM4fbbbz/xpueAAQOoV68ew4cP56qrrqJ06dIkJiZme4777ruP\n1NRUWrZsiXOOKlWqMGnSJL755hsGDx5M8eLFKVOmDOPGjWPTpk3cfffdJ66uX3rppZPOFRMTw/vv\nv88tt9xCeno6iYmJ/OEPf8h1f+677z5Wr15Ns2bNKF68OPfffz8PPfQQffv25d5776VPnz4n3ojN\nzs0338wjjzxCnz59TrT16dOH3r1706xZMzIzM6lVqxb//ve/c12XSMQK8hTKM7HjV17ZHmBWHaju\nnFtkZmWBhcBvgR7ALufcIDN7GjjXOfeUmV0NPIwX9G2A151zbc70ORISEtypG4+sWLGChg0b5rFb\nhWv//v2UKVMG5xy9evWibt26PProo36XFdLC6esrEqrMbKFzLiGn43L8ncE5t/n4Fblzbh+wAqgB\nXA+MDRw2Fi/8CbSPc565QIXAD4uINWrUKOLj42ncuDF79uyhZ8+efpckInLCWb0Za2ZxQAtgHlDN\nObc58NIWvKEd8H4IbMzyYWmBts1EqEcffVRX8CISsnL9LoCZlQE+BXo75/Zmfc154z9nHgP69fke\nMLNkM0vevr2A58iLiMgJubqiN7PieCH/oXPus0DzVjOr7pzbHBia2RZo3wRckOXDYwNtJ3HOjQRG\ngjdGn8f6RURCi49r2mQnxyv6wCya94AVzrmhWV6aDNwVeHwX8HmW9u7maQvsyTLEIyIS2Xxc0yY7\nubmiTwLuBJaY2fG7g54FBgETzexeYAPQLfDaFLwZN2vxplfeXaAVi4jIWckx6J1zM4Hsdoq4/DTH\nO6BXPusKOf369aNMmTI8/vjj2R4zadIk6tWrR6NGjYJSw+7duxk/fjwPPvhgUM4vIpFJd8YWoEmT\nJrF8+fKgnX/37t0MHz48aOcXkSD7NvnXf2YG/37SyAj6mYuC8g84cOBA6tWrR4cOHVi1atWJ9lGj\nRpGYmEjz5s256aabOHjwILNnz2by5Mk88cQTxMfHs27dutMed6pvv/2W+Ph44uPjadGiBfv27QNg\n8ODBJ5Yz7tu3LwBPP/0069atIz4+nieeeCJffROREKH16HMpCG9+LFy4kAkTJpCSksKUKVNYsGDB\nidduvPFGFixYwOLFi2nYsCHvvfce7du357rrrmPw4MGkpKRQp06d0x53qiFDhvDWW2+RkpLCjBkz\nKFmyJNOmTWPNmjXMnz+flJQUFi5cyPTp0xk0aBB16tQhJSWFwYMH57lvIhJEhbB2zdkKj9UrfTBj\nxgxuuOEGSpUqBXhL9B63dOlSnnvuOXbv3s3+/ftPuxplbo9LSkriscce44477uDGG28kNjaWadOm\nMW3aNFq0aAF4SyysWbOGmjVrBqGnIlKgsptC+W3y6dsLgYI+D3r06MGkSZNo3rw5Y8aM4Ztvvsnz\ncU8//TTXXHMNU6ZMISkpialTp+Kc45lnnvnVUgqnW1deRCQnofc7Rojo2LEjkyZN4tChQ+zbt49/\n/etfJ17bt28f1atX59ixY3z44Ycn2suWLXtijP1Mx2W1bt06mjZtylNPPUViYiIrV66kS5cujB49\n+sRWhZs2bWLbtm2/Or+ISG7oij4bLVu25NZbb6V58+ZUrVr1pOWHX3zxRdq0aUOVKlVo06bNifC9\n7bbbuP/++3njjTf45JNPsj0uq9dee42vv/6aIkWK0LhxY7p27UqJEiVYsWIF7dq1A6BMmTL8/e9/\np06dOiQlJdGkSRO6du2qcXqRcFK0SPZ3zAZZjssUF4Z8L1Mcgrccy5lpmWKR/MvtMsWRcUWvMBcR\nyZbG6EVEIlxIB30oDCtJwdPXVaLaT/thy4FC/ZQhG/QxMTHs3LlToRBhnHPs3LmTmJgYv0sRKVwH\njsEr86Dt3+HF2YX6qUN2jD42Npa0tDS0KUnkiYmJITY21u8yRApHpoNPVsGAObD5AFx3ETx5xm20\nC1zIBn3x4sWpVauW32WIiOTd3J+gz0xI2QbxVWFkF2h7fqGXEbJBLyIStjbshf6zYfJaqF4a3roC\nbq4PRbJb8T24FPQiIgVl31F4LRneWQxFDZ5sDQ+2gNLFfS1LQS8ikl8ZmTB+Bbw0F7Yfgm714S/t\n4PwyflcGKOhFRPJn+kZ4fiYs2wmtq8OH10KLan5XdRIFvYhIXqz7GfrOgqmpULMsjOoC118E5s84\n/Jko6EUkOuV1jazdh2HIAnhvCcQUhefaQc/mEBO6cRq6lYmIBNPZ7kx3LAPGLIXB82HPUbijITzd\nFqqWCl6NBURBLyJyJs7BFxug3yxY8zN0jIX+HaBxZb8ryzUFvYhIdlbs9N5o/WYj1C4PH1wDXeJC\nchz+TBT0IiKn2nEIXp4H45ZB2XPgxQ5wT1M4p6jfleWJgl5E5LhjmTBlC9y9CA4eg7ubeDc9VSzp\nd2X5oqAXkeiUdWs/52DeLvjgR9h6BK64EF5IgnoV/a2xgCjoRSQ6HZ9CuXibt/DYnJ+gQUUY1gUu\nrelvbQVMQS8i0WnLfhg4Fz5eCRVjYHAn+H0jKBay23TkmYJeRKLLwWMwIgXeWOTNje/VAh5NgHIl\n/K4saBT0IhIdnIPPVsOLc2DTfrimNvRNglrl/a4s6HL8HcXMRpvZNjNbmqWtuZnNMbMlZvYvMyuX\n5bVnzGytma0ysy7BKlxEJNcWbIaun8AfvvBm0Ey6AcZcHRUhD7nbM3YMcNUpbe8CTzvnmgL/BJ4A\nMLNGwG1A48DHDDez8Jx4KiLhb+NeeGAqXP0pbNwHb1wOX3aDpBp+V1aochy6cc5NN7O4U5rrAdMD\nj78ApgJ9gOuBCc65I8APZrYWaA3MKaiCRURytP+oNwY/4jvv+WMJ8HBLKHOOv3X5JK9j9MvwQn0S\ncAtwQaC9BjA3y3FpgTYRkeDLyIQJK+Gvc2HbQbipnre6ZGxZvyvzVV6D/h7gDTPrA0wGjp7tCczs\nAeABgJo1I2vOqoj4YFYaPDcTlu6AhGow9mpIOM/vqkJCnoLeObcSuBLAzOoB1wRe2sQvV/cAsYG2\n051jJDASICEhweWlDhER1u/2NuL+z3qoUQZGXgm/rRt2C48FU56C3syqOue2mVkR4Dng7cBLk4Hx\nZjYUOB+oC8wvkEpFRLLacwSGLoBR30PxovBMW/hjPJTUrPFT5fgvYmYfAZ2AymaWBvQFyphZr8Ah\nnwHvAzjnlpnZRGA5kA70cs5lBKNwEYlS6ZneqpKvzINdh+H2hl7In1fa78pCljnn/6hJQkKCS05O\n9rsMEQl1/9sAz8+CVbu8KZL9O0CzKn5X5RszW+icS8jpOP2OIyKhb/UuL+C/2gBx5b03WrvW0jh8\nLinoRSR07TwEr8yHsUuhdHFv6eB7m0EJ3Yd5NhT0IhJ6jmbA6CUwZAHsOwp3NYYn20Dl8N4AxC8K\nehEJHc7Bf3/wNuJev8dbF75/EjSo5HdlYU1BLyKhYekO6DMDZm6CeufCR7/xdno6GzMX/bJrVFZF\ni/yy0UgUUtCLiL+2HoBB8+DD5XBuDAzqCN0be3Pjz9bpQv5M7VFCQS8i/jicDm+nwGsLvTH5ns3h\nz4lQIcbvyiKOgl5ECpdzMGktvDjbWzq4ay1vA5A6FfyuLGIp6EWk8Cza6m3EPX8zNK4Mn10OF8f6\nU0sUjecr6EUk+H7aDwPmwD9WQZVSMPRS+F1DL1T9EkXj+Qp6EQmeA8fgzUXw1neQ6eCRVtC7VfA2\nAClaJPur9CimoBeRgpfpYOJKGDgXthzwlg3u0w5qlsv5Y/MjwoZcCoqCXkQK1tyfvHH4lG3Qshq8\ndxW0ru53VVFNQS8iBWPDXnhhFvxrHZxfBoZ39rbyK6KFx/ymoBeR/Nl3FF5NhndSoFgReKo1PNgC\nShX3u7Izi6LxfAW9iORNRqZ3N+ugebD9ENzaAP7SFqqX8buy3Imi8XwFvYicvW83wvMzYflOaFMd\nxl8L8dX8rkqyoaAXkdxb+zP0nQXTUuHCct4brb+pow1AQpyCXkRy9vNhb2340Usgpig83x7ubwYx\nipBwoK+SSKgKhVv0j2XA+0th8HzYexTubARPtfHubpWwoaAXCVV+3qLvHHyR6g3TrN0Nl1zgbQDS\nqHLwP7cUOAW9iJxs+Q5vI+5vN3orSn54DXSO0zh8GFPQi4hn+0F4eR58sBzKnQMDL4a7m+RtAxAJ\nKQp6kWh3JANGLvZuejqUDvc2hSdae7s9SURQ0ItEK+e85Qr6z/aWL7gyDl5IgovO9bsyKWAKepFQ\nFcxb9Bdvg+dmeguQNawIn1zvveEqEUlBLxKqgjGFcst+GDDXW0K4UkkY0gnuaOStUSMRS0EvEg0O\nHoPh38GwRZCeCQ+19DYAKVfC78qkECjoRSJZpoNPV3vb+P2031uu4Pn2EFfe78qkECnoRSLV/M3e\nBiCLtkLzKvD2ldDufL+rEh8o6EUizca90H8OTFoD55WGYZdDtwbaACSK5Rj0ZjYauBbY5pxrEmiL\nB94GYoB04EHn3HwzM+B14GrgINDDObcoWMWLSBb7j8JrC+HtFC/UH0/0xuJLn7IBSCisoSOFKjdX\n9GOAN4FxWdpeAV5wzv2fmV0deN4J6ArUDfxpA4wI/C0iwZKRCR+thL/O9e5uvaW+twFIjbLZH382\n7RL2cgx659x0M4s7tRk4vp17eeCnwOPrgXHOOQfMNbMKZlbdObe5gOoVkaxmpnnz4ZftgMTz4O/X\neBtyi2SR1zH63sBUMxsCFAHaB9prABuzHJcWaPtV0JvZA8ADADVr1sxjGSJRat1ubyPu//sBLigL\nI7vAby/SwmNyWnkN+j8CjzrnPjWzbsB7wBVncwLn3EhgJEBCQoLLYx0i0WXPEfjbAnj3ezinqDdE\n0zMeSmpehWQvr98ddwGPBB7/A3g38HgTkPU+6thAm4jkR3omjF0Kr8z3dnu6oxE83Qaqlfa7MgkD\neQ36n4BLgG+Ay4A1gfbJwENmNgHvTdg9Gp8XyacvN0DfmbD6Z+hQA/p3gKZV8n6+YK6hIyEpN9Mr\nP8KbUVPZzNKAvsD9wOtmVgw4TGCsHZiCN7VyLd70yruDULNIdFi509sA5OsfoVZ5GHc1XFUr/+Pw\nmkIZdXIz6+b2bF5qdZpjHdArv0WJRLWdh+Dl+TBuKZQ5B17sAPc09cbks6O58XIGegdHJFQczfDe\nZP3bAjhwDHo08TYAqVQy54/V3Hg5AwW9iN+cgynrod9sSN0Dl1/obQBSv6LflUmEUNCL+On77fD8\nTJi1yQv2j38Dl13od1USYRT0In7YcgBemgsfrYCKMfDKJXBnY20AIkGhoBcpTIfSYUQKvL4QjmXA\ngy3g0QQorw1AJHgU9CKFwTlv2eD+cyBtH1xT29sApHaFgjm/5sbLGSjoRYJt4RZvA5AFW7wbnd68\nHJJiC/ZzaAqlnIGCXiRYNu2DF+d4W/lVLQWvXwa3NtBVthQ6Bb1IQdt/1NuEe/h33vNHE+BPLb2b\nn0R8oKAXKSiZDj5eCQPnwNaDcGNd6NMeYrPZAESkkCjoRQrC7E3eOPz326FVNXi/KyRW97sqEUBB\nL5I/P+zxNgD5z3qoUQbe7gw31tMGIBJSFPQiebH3CAxNhlGLoVhRb234P8ZDqeI5f6xIIVPQi5yN\n9Ez4+3J4eZ63yuRtDeHZNnBeGb8rE8mWgl4kt77+0VuXZuUuaHc+TPgNNK/qd1UiOVLQi+Rkzc9e\nwH+5AeLKeW+0XlP7zOPw4bY+fLjVK2dFQS+SnV2HYPACeH+JN/berz3c1xxKnGEDkOPCbX34cKtX\nzoqCXuRURzNg9BJvA5C9R6F7Y3iyNVQp5XdlInmioBc5zjmYluoN06zfA50u8DbibljJ78pE8kVB\nLwKwbIcX8NPToO65MP5auOJCzYeXiKCgl+i27SAMmgsfroDy58BLHeGuxlA8F+PwImFCQS/R6XA6\njFwMrybD4Qy4vxk8nggVYgrm/OG2Pny41StnRUEv0cU5mLwO+s+CH/fBVbW82TR1zi3YzxNuUxLD\nrV45Kwp6yVmkzLFO2QrPzYR5m6FxJfj0euh4gd9ViQSdgl5yFu5zrH/a7y0dPHEVVCkJQy+F3zXU\nsIREDQW9RK4Dx7zNP4Yt8taK/1NL6J0AZbUBiEQXBb1EnkwHn6yCAXNg8wG47iJvI+4Ly/ldmYgv\nFPQSWeZthj4z4LttEF8VRnaBtuf7XZWIrxT0Ehl+3Av9Z8Pna6F6aXjrCri5PhTRDU8iCnrJWSjP\nsd53FF5LhncWe6H+RCL0agmltQGIyHE5Br2ZjQauBbY555oE2j4G6gcOqQDsds7FB157BrgXyAD+\n5JybGozCpRCF4hTKjEwYvwJemgvbD0G3+vCXdnC+NgAROVVurujHAG8C4443OOduPf7YzP4G7Ak8\nbgTcBjQGzge+NLN6zrmMAqxZot2MNG8cftlOaF0dPrwWWlTzuyqRkJVj0DvnpptZ3OleMzMDugGX\nBZquByY4544AP5jZWqA1MKdAqpXotu5n6Dcb/vsD1CwLo7rA9Rdp4TGRHOR3jP5iYKtzbk3geQ1g\nbpbX0wJtInm3+zAMWQDvLYGYovBcO+jZHGL0FpNIbuT3f8rtwEd5+UAzewB4AKBmzZr5LENCRkEu\nl3AsA8Yug1fmwe4j8PtG8HRbqKoNQETORp6D3syKATcCrbI0bwKyLh4SG2j7FefcSGAkQEJCgstr\nHRJiCmK5BOe8/Vn7zvL2a7041tsApEnlgqlRJMrk54r+CmClcy4tS9tkYLyZDcV7M7YuMD8fn0Oi\nzcqd0GcmfLMRapeHD66BLnEahxfJh9xMr/wI6ARUNrM0oK9z7j282TUnDds455aZ2URgOZAO9NKM\nG8mVHYe8IZqxy7y1aF7sAPc0hXO0AYhIfuVm1s3t2bT3yKZ9IDAwf2VJ1DiSAe8uhqHJ3iJk9zT1\nbnqqWNLvykQihqYtiD+cgynrvemSqXu8/VlfSIJ6Ff2uTCTiKOilYOVmuYTvt3s3PM3+CRpUhInX\nwaWaeSUSLAp6KVhnmkK55QD8dS5MWAEVY2BwJ2/KZLEQWDNHJIIp6CX4DqXDiBR4faE3N75XC3g0\nAcqV8LsykaigoJfgcQ7+uQZenANp++Ca2tA3CWqV97sykaiioJfgSN7ijcMnb4WmVeDNKyBJq2GI\n+EFBLwUrbR+8OBs+W+MtVfD6ZXBbQ20AIuIjBb0UjP1HvU24h3/nPX8sAR5uCWW0EbeI3xT0kj+Z\nzptF89e5sPUg3FTPW10ytqzflYlIgIJe8m7WJm9dmiXbIaEajLkaEs7zuyoROYWCXs7eD3vghVnw\nn/VQowy8cyXcUFcLj4mEKAW95N7eI96aNKMWQ7Gi8Exb+GM8lNS3kUgo0/9QyVl6JnywDF6eB7sO\ne7Nonm0L55X2uzIRyQUFvZzZ/wIbgKzcBe3P9zYAaV7V76pE5Cwo6OX0Vu/yAv7LDRBXHsZ0hatr\naxxeJAwp6OVkuw7BK/NhzFIoVRz6tYf7mkMJbQAiEq4U9OI5mgGjl8CQBbDvKNzVGJ5sA5W1AYhI\nuFPQRzvnYGoq9J0J6/dApwu8bfwaVPK7MhEpIAr6aLZ0Bzw/E2akQd1zYfy13k5PGocXiSgK+mi0\n7SC8NBc+XA4VSsBLHb2hmuIahxeJRAr6aHI4Hd5ZDK8lw+EM6Nkc/pwIFWL8rkxEgkhBHw2cg8/X\nessH/7gPutbyNgCpU8HvykSkECjoI913W+G5mTB/MzSuBJ/9Fi6O9bsqESlECvpI9dN+GDAH/rEK\nqpSCoZfC7xpCUW3ELRJtFPSR5sAxeGsRvPmdt1b8I62gdyttACISxRT0kSLTeVfvA+fA5gNw/UXw\nfHuoWc7vykTEZwr6SDD3J28DkJRt0KIqjLoK2lT3uyoRCREK+nC2YS/0nw2T10L10jC8s7eVnzbi\nFpEsFPThaN9ReDUZ3kmBYkXgqdbwYAtvETIRkVMo6MNJRqZ3N+ugebD9EHSr723EXb2M35WJSAjL\nca6dmY02s21mtvSU9ofNbKWZLTOzV7K0P2Nma81slZl1CUbRUWn6RrjsY/jzN1C7AnxxC7zVWSEv\nIjnKzRX9GOBNYNzxBjO7FLgeaO6cO2JmVQPtjYDbgMbA+cCXZlbPOZdR0IVHjXU/exuATE2FC8vB\ne1fBb+po4TERybUcg945N93M4k5p/iMwyDl3JHDMtkD79cCEQPsPZrYWaA3MKbCKo8XPh7214Ucv\ngZii3lTJ+5tBjEbbROTs5DU16gEXm9lA4DDwuHNuAVADmJvluLRAm+TWsQx4fykMmQ97jsLvG8FT\nbaBqKb8rE5EwldegLwZUBNoCicBEM6t9NicwsweABwBq1qyZxzIiiHPwxQZvA5C1u+GSC6B/EjSq\n7HdlIhLm8hr0acBnzjkHzDezTKAysAm4IMtxsYG2X3HOjQRGAiQkJLg81hEZVuz0NgD5ZqO3ouSH\n10DnOI3Di0iByGvQTwIuBb42s3rAOcAOYDIw3syG4r0ZWxeYXxCFRqTtB+HlefDBcih3Dgy4GO5u\nAudoAxARKTg5Br2ZfQR0AiqbWRrQFxgNjA5MuTwK3BW4ul9mZhOB5UA60Eszbk7jSAaMWgxDk+FQ\nOtzbFJ5oDedqAxARKXjm5bO/EhISXHJyst9lBJ9z8O913rIFqXuh84XwQgdvv1YRkbNkZgudcwk5\nHae5eoVl8TZv4bE5P0HDivCP66CT3oQWkeBT0Afblv0wYC5MXAmVSsKQTnBHI2+NGhGRQqCgD5aD\nx2BECryxENIz4aGW3gYg5Ur4XZmIRBkFfUHLdPDZanhxjred37V1oG97iCvvd2UiEqUU9AVpwWZv\nI+5FW6FZFRjRGdrrxmAR8ZeCviBs3Av958CkNXBeaRh2OXRroA1ARCQkKOjzY/9ReH2hNxZfxODx\nRG8svrTPPI9HAAAHIUlEQVQ2ABGR0KGgz4uMTJiwEv46F7YdhJvreRuA1Cjrd2UiIr+ioD9bs9K8\ncfilOyDxPBh3NbQ67/THzlzk/VA4VdEi0KFlcOsUEQlQ0OfW+t3wwmyYsh5iy8LILvDbi8688Njp\nQv5M7SIiQaCgz8meI/C3BfDu995iY8+2hT/EQ0n904lIeFBaZSc9E8Yt81aX/Pkw/K4RPNMGqpX2\nuzIRkbOioD+drzZ4+7Su2gUdakD/DtC0it9ViYjkiYI+q1W7vA1A/vcj1CrvvdF6VS1tACIiYU1B\nD7DzELw8H8Yt9ebA90+Ce5vlfwOQokWyn3UjIlJIojvoj2bAe9/DkAVw4Bh0bwJPtfZWmSwImkIp\nIiEgOoPeOfi/H6DfLPhhD1xW0xuHr1/R78pERApc9AX9ku3eBiCzNnnBPuE3cPmFflclIhI00RP0\nWw/AS/Ng/HJvb9aXL4HujbUBiIhEvPAP+pyWGTiUDm+neIuPHc3wbnb6cyKU1wYgIhIdwj/os1tO\nID0D/rnaWz44bR9cXRuebw91KhRufSIiPgv/oD+d1ftgzAZYvR+aVPbWh+8Q63dVIiK+iKyg33EE\nPtwIM3ZAheLw6mVwewPNWxeRqBY5QT9jB4xY702dvPF8uKEGXNXI76pERHwXOUEfVwpanwt31IQq\neqNVROS48A/648sMXFAKetc9uV1ERCIg6LXMgIjIGemyV0QkwinoRUQinIJeRCTC5Rj0ZjbazLaZ\n2dIsbf3MbJOZpQT+XJ3ltWfMbK2ZrTKzLsEqXEREcic3V/RjgKtO0/6qcy4+8GcKgJk1Am4DGgc+\nZriZ5XP3DhERyY8cg945Nx3YlcvzXQ9McM4dcc79AKwFWuejPhERyaf8jNE/ZGbfB4Z2zg201QA2\nZjkmLdAmIiI+yWvQjwDqAPHAZuBvZ3sCM3vAzJLNLHn79u15LENERHKSpxumnHNbjz82s1HAvwNP\nNwEXZDk0NtB2unOMBEYGzrHdzDbkpZZ8qAzsKOTPWZgiuX/qW/iK5P750bdcbY+Xp6A3s+rOuc2B\npzcAx2fkTAbGm9lQ4HygLjA/p/M556rkpY78MLNk51xCYX/ewhLJ/VPfwlck9y+U+5Zj0JvZR0An\noLKZpQF9gU5mFg84IBXoCeCcW2ZmE4HlQDrQyzmXEZzSRUQkN3IMeufc7adpfu8Mxw8EBuanKBER\nKTjRfGfsSL8LCLJI7p/6Fr4iuX8h2zdzzvldg4iIBFE0X9GLiESFqAh6M7vAzL42s+VmtszMHgm0\nVzSzL8xsTeDvc3M6V6gxsxgzm29miwN9eyHQXsvM5gXWHfrYzM7xu9a8MrOiZvadmf078DyS+pZq\nZksCa0YlB9rC/vsSwMwqmNknZrbSzFaYWbtI6JuZ1c+yzleKme01s96h3LeoCHq8GUB/ds41AtoC\nvQLr8jwNfOWcqwt8FXgebo4AlznnmuPdwHaVmbUFXsZbj+gi4GfgXh9rzK9HgBVZnkdS3wAuDawZ\ndXxqXiR8XwK8DvzXOdcAaI73NQz7vjnnVh1f5wtoBRwE/kko9805F3V/gM+BzsAqoHqgrTqwyu/a\n8tmvUsAioA3ejRvFAu3tgKl+15fHPsXi/ae5DO/GPIuUvgXqTwUqn9IW9t+XQHngBwLvA0ZS307p\nz5XArFDvW7Rc0Z9gZnFAC2AeUM39cuPXFqCaT2XlS2BoIwXYBnwBrAN2O+fSA4eE85pDrwFPApmB\n55WInL6Bdy/KNDNbaGYPBNoi4fuyFrAdeD8w7PaumZUmMvqW1W3AR4HHIdu3qAp6MysDfAr0ds7t\nzfqa834Mh+UUJOdchvN+jYzFWy20gc8lFQgzuxbY5pxb6HctQdTBOdcS6Io3pNgx64th/H1ZDGgJ\njHDOtQAOcMpQRhj3DYDAe0PXAf849bVQ61vUBL2ZFccL+Q+dc58FmreaWfXA69XxrojDlnNuN/A1\n3nBGBTM7fkNctmsOhbgk4DozSwUm4A3fvE5k9A0A59ymwN/b8MZ5WxMZ35dpQJpzbl7g+Sd4wR8J\nfTuuK7DI/bL2V8j2LSqC3swM727eFc65oVlemgzcFXh8F97YfVgxsypmViHwuCTeew8r8AL/5sBh\nYdk359wzzrlY51wc3q/I/3PO3UEE9A3AzEqbWdnjj/HGe5cSAd+XzrktwEYzqx9ouhxvaZSw71sW\nt/PLsA2EcN+i4oYpM+sAzACW8MtY77N44/QTgZrABqCbcy63m6yEBDNrBowFiuL94J7onOtvZrXx\nroIrAt8Bv3fOHfGv0vwxs07A4865ayOlb4F+/DPwtBgw3jk30MwqEebflwCB9bDeBc4B1gN3E/ge\nJfz7Vhr4EajtnNsTaAvZr1tUBL2ISDSLiqEbEZFopqAXEYlwCnoRkQinoBcRiXAKehGRCKegFxGJ\ncAp6EZEIp6AXEYlw/w9af6UxYIrGWgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10868b780>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "xs = np.r_[min(xd):max(xd):15j]\n",
    "ys = a * xs + b\n",
    "\n",
    "plt.plot(xs, ys, color='deeppink', label='regression curve')\n",
    "plt.scatter(xd, yd, color='pink', marker='s', label='data set')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 正規直交化\n",
    "http://www.python.ambitious-engineer.com/archives/1279\n",
    "\n",
    "QR分解\n",
    "正則行列Aが直交行列Qと上三角行列R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "a = np.array([[1, 1], [1, -1], [0, 0]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1,  1],\n",
       "       [ 1, -1],\n",
       "       [ 0,  0]])"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "q, r = np.linalg.qr(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.70710678, -0.70710678],\n",
       "       [-0.70710678,  0.70710678],\n",
       "       [-0.        , -0.        ]])"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ -1.41421356e+00,   3.33066907e-16],\n",
       "       [  0.00000000e+00,  -1.41421356e+00]])"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "q, r = linalg.qr(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.70710678, -0.70710678,  0.        ],\n",
       "       [-0.70710678,  0.70710678,  0.        ],\n",
       "       [-0.        , -0.        ,  1.        ]])"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ -1.41421356e+00,   3.33066907e-16],\n",
       "       [  0.00000000e+00,  -1.41421356e+00],\n",
       "       [  0.00000000e+00,   0.00000000e+00]])"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}