{ "cells": [ { "cell_type": "markdown", "metadata": { "toc": "true" }, "source": [ "# Table of Contents\n", "

1  固有値
2  Googleのページランク
2.1  課題
3  累乗(べき乗)法により最大固有値が求まる原理
4  Jacobi回転による固有値の求め方
4.1  Mapleでみる回転行列
4.2  C 言語での実装
5  数値計算ライブラリーについて
5.1  python謹製 lapack利用逆行列
5.2  計算速度競争
6  課題
6.1  行列のpprint設定イロイロ
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "
\n", "線形代数-固有値(LAEigen)\n", "
\n", "
\n", "
\n", "file:/Users/bob/Github/TeamNishitani/jupyter_num_calc/laeigenvectors\n", "
\n", "https://github.com/daddygongon/jupyter_num_calc/tree/master/notebooks_python\n", "
\n", "cc by Shigeto R. Nishitani 2017 \n", "
\n", "\n", "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 固有値\n", "\n", "$A$を対称正方行列，$x$をベクトルとしたときに，\n", "$$\n", "Ax = \\lambda x\n", "$$\n", "の解，$\\lambda$を固有値，$x$を固有ベクトルという．$x$がゼロベクトルではない意味のある解は特性方程式det$(A-\\lambda E)=0$が成り立つときにのみ得られる．\n", "\n", "この特性方程式をpythonで特には，sympyってのを使わないといけない．一応示しておくと，\n", "\n", "$$\n", "\\begin{align}\n", " A\\, &= \\, \\left( \\begin {array}{cc} 3&2/3\\\\2/3&2\\end {array} \\right) \\\\\n", " \\det( A\\, - t E) &= \n", "\\, \\det \\left( \\begin {array}{cc} 3-t&2/3\\\\2/3& 2-t \\end {array} \\right) \\\\\n", "& = t^2 -5t + 6 - 4/9\n", "\\\\\n", "& = t^2 -5t + 50/9 \\\\\n", "t (\\lambda) &= 10/3,\\,5/3\n", "\\end{align}\n", "$$\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "⎡t 0⎤\n", "⎢ ⎥\n", "⎣0 t⎦\n", " 2 \n", "t - 5⋅t + 5.55555555555556\n" ] }, { "data": { "text/plain": [ "[1.66666666666667, 3.33333333333333]" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from sympy import *\n", "t = symbols('t')\n", "A = Matrix([[3, 2/3],[2/3,2]]\n", ")\n", "xx = diag(t,t)\n", "pprint(xx)\n", "eq1 = simplify((A-xx).det())\n", "pprint(eq1)\n", "solve(eq1,t)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "固有値を求めるコマンドEigenvectorsを適用すると，固有値と固有ベクトルが求まる．ここで，固有ベクトルは行列の列(Column)ベクトルに入っている．" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "array([ 3.33333333+0.j, 1.66666667+0.j])\n", "array([[ 0.89442719, -0.4472136 ],\n", " [ 0.4472136 , 0.89442719]])\n", "array([ 0.89442719, 0.4472136 ])\n" ] } ], "source": [ "import pprint\n", "import scipy.linalg # SciPy Linear Algebra Library\n", "\n", "A = scipy.array([[3, 2/3],[2/3,2]])\n", "l, V = scipy.linalg.eig(A)\n", "pprint.pprint(l)\n", "pprint.pprint(V)\n", "pprint.pprint(V[:,0]) # columnの取り出し方" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "固有ベクトルは規格化されている．" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Googleのページランク\n", "\n", ">多くの良質なページからリンクされているページはやはり良質なページである\n", "\n", "Googleのpage rankは上のような非常に単純な仮定から成り立っている．ページランクを実際に求めよう．つぎのようなリンクが張られたページを考える．\n", "\n", "![link structure](figs/linkstruct.png)\n", "\n", "\n", "計算手順は以下の通り\n", "1. リンクを再現する**隣接行列**を作る．ページに番号をつけて，その間が結ばれているi-j要素を1，そうでない要素を0とする．\n", "1. 隣接行列を**転置**する\n", "1. 列ベクトルの総和が1となるように**規格化**する．\n", "1. こうして得られた**推移確率行列**の最大固有値に属する固有ベクトルを求め，適当に規格化する．\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 課題 \n", "\n", "1. 上記手順を参考にして，pythonでページランクを求めよ．\n", "1. このような問題ではすべての固有値・固有ベクトルを求める必要はなく，最大の固有値を示す固有ベクトルを求めるだけでよい．初期ベクトルを適当に決めて，何度も推移確率行列を掛ける反復法でページランクを求めよ．\n", "
\n", "
隣接行列
\n", "
\n", "\n", "$$\n", "{\\it A1}\\, := \\, \\left[ \\begin {array}{c|c|c|c|c|c|c|c} \n", "&1&2&3&4&5&6&7\\\\\n", "1&0&1&1&1&1&0&1\\\\\n", "2&1&0&0&0&0&0&0\\\\\n", "3& & & & & & & \\\\\n", "4& & & & & & & \\\\\n", "5& & & & & & & \\\\\n", "6& & & & & & & \\\\\n", "7& & & & & & & \n", "\\end {array} \\right] \n", "$$\n", "
\n", "
転置行列
\n", "
\n", "\n", "$$\n", "{Transpose}({\\it A1})\\, := \\, \\left[ \\begin {array}{c|c|c|c|c|c|c} \n", "\\, \\, &\\, \\, &\\, \\, &\\, \\, &\\, \\, &\\, \\, &\\, \\, \\\\\n", "& & & & & & \\\\\n", "& & & & & & \\\\\n", "& & & & & & \\\\\n", "& & & & & & \\\\\n", "& & & & & & \\\\\n", "& & & & & & \n", "\\end {array} \\right] \n", "$$\n", "
\n", "
規格化
\n", "
\n", "\n", "$$\n", "\\left[ \\begin {array}{c|c|c|c|c|c|c} \n", "\\, \\, &\\, \\, &\\, \\, &\\, \\, &\\, \\, &\\, \\, &\\, \\, \\\\\n", "& & & & & & \\\\\n", "& & & & & & \\\\\n", "& & & & & & \\\\\n", "& & & & & & \\\\\n", "& & & & & & \\\\\n", "& & & & & & \n", "\\end {array} \\right] \n", "$$\n", "
\n", "
遷移
\n", "
\n", "\n", "$$\n", "\\left( \\begin {array}{ccccccc} \n", "0 &1 &1/2 &0 &1/4 &1/2 &0 \\\\\n", "1/5 &0 &1/2 &1/3 &0 &0 &0 \\\\\n", "1/5 &0 &0 &1/3 &1/4 &0 &0 \\\\\n", "1/5 &0 &0 &0 &1/4 &0 &0 \\\\\n", "1/5 &0 &0 &1/3 &0 &1/2 &1 \\\\\n", "0 &0 &0 &0 &1/4 &0 &0 \\\\\n", "1/5 &0 &0 &0 &0 &0 &0 \n", "\\end {array} \\right) \n", "\\left( \\begin {array}{c} \n", "1/7\\\\ \n", "1/7\\\\ \n", "1/7\\\\ \n", "1/7\\\\ \n", "1/7\\\\ \n", "1/7\\\\ \n", "1/7 \n", "\\end {array} \\right) \\, = \\, \n", "\\left( \\begin {array}{ccccccc} \n", "\\, \\, &\\, \\, &\\, \\, &\\, \\, &\\, \\, &\\, \\, &\\, \\, \\\\\n", "& & & & & & \\\\\n", "& & & & & & \\\\\n", "& & & & & & \\\\\n", "& & & & & & \\\\\n", "& & & & & & \\\\\n", "& & & & & & \n", "\\end {array} \\right)\n", "\\, = \\, \\left( \\begin {array}{c} \n", "0.32\\\\ \n", "0.15\\\\ \n", "0.11\\\\ \n", "0.06\\\\ \n", "0.29\\\\ \n", "0.04\\\\ \n", "0.03 \n", "\\end {array} \\right) \n", "$$\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "array([ 0.35 , 0.16666667, 0.11666667, 0.05 , 0.26666667,\n", " 0.05 , 0. ])\n", "array([ 0.69945653+0.j, 0.38286042+0.j, 0.32395882+0.j, 0.24296911+0.j,\n", " 0.41231122+0.j, 0.10307780+0.j, 0.13989131+0.j])\n" ] } ], "source": [ "from pprint import pprint\n", "from numpy import array, zeros, diagflat, dot, transpose\n", "from scipy.linalg import eig\n", "\n", "A = array([[0,1,1,1,1,0,1],\n", " [1,0,0,0,0,0,0],\n", " [1,1,0,0,0,0,0],\n", " [0,1,1,0,1,0,0],\n", " [1,0,1,1,0,1,0],\n", " [1,0,0,0,1,0,0],\n", " [0,0,0,0,1,0,0]])\n", "\n", "diag = []\n", "for i in range(0,7):\n", " tmp = 0.0\n", " for j in range(0,7):\n", " tmp += A[i,j]\n", " diag.append(1.0/tmp)\n", "\n", "D = diagflat(diag)\n", "tA = dot(transpose(A),D)\n", "\n", "x = array([1,0,0,0,0,0,0])\n", "pprint(dot(tA,dot(tA,x)))\n", "\n", "l, V = eig(tA)\n", "v0 = V[:,0]\n", "pprint(v0)\n" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "array([[ 0.200, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000],\n", " [ 0.000, 1.000, 0.000, 0.000, 0.000, 0.000, 0.000],\n", " [ 0.000, 0.000, 0.500, 0.000, 0.000, 0.000, 0.000],\n", " [ 0.000, 0.000, 0.000, 0.333, 0.000, 0.000, 0.000],\n", " [ 0.000, 0.000, 0.000, 0.000, 0.250, 0.000, 0.000],\n", " [ 0.000, 0.000, 0.000, 0.000, 0.000, 0.500, 0.000],\n", " [ 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 1.000]])\n", "array([[ 0.000, 1.000, 0.500, 0.000, 0.250, 0.500, 0.000],\n", " [ 0.200, 0.000, 0.500, 0.333, 0.000, 0.000, 0.000],\n", " [ 0.200, 0.000, 0.000, 0.333, 0.250, 0.000, 0.000],\n", " [ 0.200, 0.000, 0.000, 0.000, 0.250, 0.000, 0.000],\n", " [ 0.200, 0.000, 0.000, 0.333, 0.000, 0.500, 1.000],\n", " [ 0.000, 0.000, 0.000, 0.000, 0.250, 0.000, 0.000],\n", " [ 0.200, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000]])\n" ] } ], "source": [ "import numpy as np\n", "from pprint import pprint\n", "\n", "np.set_printoptions(formatter={'float': '{: 0.3f}'.format}) \n", "#pprint(A)\n", "pprint(D)\n", "#pprint(transpose(A))\n", "pprint(dot(transpose(A),D))\n", "x = array([1,0,0,0,0,0,0])\n", "#pprint(x)\n", "#pprint(dot(tA,dot(tA,dot(tA,x))))\n", "l, V = eig(tA)\n", "v0 = V[:,0]\n", "#pprint(v0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 累乗(べき乗)法により最大固有値が求まる原理\n", "\n", "累乗(べき乗)法は，最大固有値とその固有ベクトルを効率的に見つける算法である．すこし，固有値について復習しておく．正方行列$A$に対して，\n", "$$\n", "A x = \\lambda x\n", "$$\n", "の解$\\lambda$を固有値，$x$を固有ベクトルという．$\\lambda$は，\n", "$$\n", "\\det( A - \\lambda E) =0\n", "$$\n", "として求まる永年方程式の解である．\n", "\n", "では，なぜ適当な初期ベクトル$x_0$から始めて，反復\n", "$$\n", "x_{k+1} = A x_k\n", "$$\n", "を繰り返すと，$A$の絶対値最大の固有値に属する固有ベクトルに近づいていくのかを見ておこう．\n", "\n", "すべての固有値がお互いに異なる場合を考える．今，行列の固有値を絶対値の大きなもの順に並べて，$|\\lambda_1|>|\\lambda_2|>\\cdots>|\\lambda_n|$とし，対応する長さを1に規格化した固有ベクトルを$x_1, x_2, \\ldots, x_n$とする．初期ベクトルは固有ベクトルの線形結合で表わせて，\n", "$$\n", "X_0 = c_1x_1+c_2x_2+\\cdots+c_nx_n\n", "$$\n", "となるとする．これに行列$A$を$N$回掛けると，\n", "$$\n", "A^N X_0 = c_1 \\lambda_1^N x_1+\n", "c_2 \\lambda_2^N x_2+\\cdots+\n", "c_n \\lambda_n^N x_n\n", "$$\n", "となる．これを変形すると，\n", "$$\n", "A^NX_0 = X_{N}\n", "= c_1 \\lambda_1^N \\left\\{ x_1+\n", "\\frac{c_2}{c_1}\\left(\\frac{\\lambda_2}{\\lambda_1}\\right)^N x_2+\\cdots+\n", "\\frac{c_n}{c_1}\\left(\\frac{\\lambda_n}{\\lambda_1}\\right)^N x_n \\right\\}\n", "$$\n", "となる．$|\\lambda_1|>|\\lambda_i|(i\\ge2)$だから括弧の中は$x_1$だけが生き残る．\n", "\n", "こうして最大固有値に属する固有ベクトルが，反復計算を繰り返すだけで求められる．\n", "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Jacobi回転による固有値の求め方\n", "\n", "固有値を求める手法として，永年方程式を解くというやり方は回りくどすぎる．少し古めかしいが非対角要素を0にする回転行列を反復的に作用させるJacobi(ヤコビ)法を紹介する．現在認められている最適の方策は，ハウスホルダー(Householder)変換で行列を単純な三重対角化行列に変形してから，反復法で解を追い込んでいくやり方である．Jacobi法は，Householder法ほど万能ではないが，10次程度までの行列には今でも役に立つ．\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Mapleでみる回転行列 \n", "\n", "行列の軸回転の復習をする．対称行列$B$に回転行列$U$を作用すると\n", "$$\n", "B.U = \n", "\\left(\n", "\\begin{array}{cc}\n", "{a_{1\\,1}} & {a_{1\\,2}}\\\\\n", "{a_{2\\,1}(={a_{1\\,2}})} & {a_{2\\,2}}\n", "\\end{array}\n", "\\right)\n", "\\left( \n", "\\begin{array}{cc}\n", "\\cos(\\theta) & -\\sin(\\theta)\\\\\n", "\\sin(\\theta) & \\cos(\\theta)\n", "\\end{array}\n", "\\right) \n", "$$\n", "\n", "となる．回転行列を4x4の行列に\n", "$$\n", "U^t B U\n", "$$\n", "と作用させたときの各要素の様子を以下に示した．\n", "```maple\n", "> restart:\n", "> n:=4:\n", "> with(LinearAlgebra):\n", "> B:=Matrix(n,n,shape=symmetric,symbol=a);\n", "```\n", "\n", "$$\n", "B := \\left[{\n", "\\begin{array}{cccc}\n", "{a_{1, \\,1}} & {a_{1, \\,2}} & {a_{1, \\,3}} & {a_{1, \\,4}} \\\\\n", "{a_{1, \\,2}} & {a_{2, \\,2}} & {a_{2, \\,3}} & {a_{2, \\,4}} \\\\\n", "{a_{1, \\,3}} & {a_{2, \\,3}} & {a_{3, \\,3}} & {a_{3, \\,4}} \\\\\n", "{a_{1, \\,4}} & {a_{2, \\,4}} & {a_{3, \\,4}} & {a_{4, \\,4}}\n", "\\end{array}}\n", "\\right] \n", "$$\n", "\n", "```maple\n", "> U:=Matrix(n,n,[[c,-s,0,0],[s,c,0,0],[0,0,1,0],[0,0,0,1]]);\n", "\n", "#U:=Matrix(n,n,[[c,-s],[s,c]]);\n", "```\n", "$$\n", "U := \\left[ \n", "{\\begin{array}{ccrr}\n", "c & - s & 0 & 0 \\\\\n", "s & c & 0 & 0 \\\\\n", "0 & 0 & 1 & 0 \\\\\n", "0 & 0 & 0 & 1\n", "\\end{array}}\n", "\\right] \n", "$$\n", "\n", "```maple\n", ">TT:=Transpose(U).B.U;\n", "```\n", "\n", "$$\n", "\\mathit{TT} := \\\\ \\notag\n", "{\\begin{array}{c}\n", "\\left[ \\right. \n", "(c\\,{a_{1, \\,1}} + s\\,{a_{1, \\,2}})\\,c + (c\\,{a_{1, \\,2}} + s\n", "\\,{a_{2, \\,2}})\\,s\\,, \\, - (c\\,{a_{1, \\,1}} + s\\,{a_{1, \\,2}})\\,s\n", "+ (c\\,{a_{1, \\,2}} + s\\,{a_{2, \\,2}})\\,c\\,, \\\\\n", "c\\,{a_{1, \\,3}} + s\\,{a_{2, \\,3}}\\,, \\,c\\,{a_{1, \\,4}} + s\\,{a_{2\n", ", \\,4}} \\left. \\right] \\\\\n", "\\left[ \\right. \n", "( - s\\,{a_{1, \\,1}} + c\\,{a_{1, \\,2}})\\,c + ( - s\\,{a_{1, \\,2\n", "}} + c\\,{a_{2, \\,2}})\\,s\\,, \\, - ( - s\\,{a_{1, \\,1}} + c\\,{a_{1, \n", "\\,2}})\\,s + ( - s\\,{a_{1, \\,2}} + c\\,{a_{2, \\,2}})\\,c\\,, \\\\\n", "- s\\,{a_{1, \\,3}} + c\\,{a_{2, \\,3}}\\,, \\, - s\\,{a_{1, \\,4}} + c\n", "\\,{a_{2, \\,4}} \\left. \\right] \\\\\n", "\\left[ c\\,{a_{1, \\,3}} + s\\,{a_{2, \\,3}}\\,, \\, - s\\,{a_{1, \n", "\\,3}} + c\\,{a_{2, \\,3}}\\,, \\,{a_{3, \\,3}}\\,, \\,{a_{3, \\,4}} \n", "\\right] \\\\\n", "\\left[ c\\,{a_{1, \\,4}} + s\\,{a_{2, \\,4}}\\,, \\, - s\\,{a_{1, \n", "\\,4}} + c\\,{a_{2, \\,4}}\\,, \\,{a_{3, \\,4}}\\,, \\,{a_{4, \\,4}} \n", "\\right] \n", "\\end{array}}\n", "$$\n", "\n", "```maple\n", ">expand(TT[1,1]);\n", "expand(TT[2,2]);\n", "expand(TT[1,2]);\n", "expand(TT[2,1]);\n", "```\n", "$$\n", "c^{2}\\,{a_{1, \\,1}} + 2\\,c\\,s\\,{a_{1, \\,2}} + s^{2}\\,{a_{2, \\,2}}\n", "$$\n", "$$\n", "s^{2}\\,{a_{1, \\,1}} - 2\\,c\\,s\\,{a_{1, \\,2}} + c^{2}\\,{a_{2, \\,2}}\n", "$$\n", "$$\n", "- s\\,c\\,{a_{1, \\,1}} - s^{2}\\,{a_{1, \\,2}} + c^{2}\\,{a_{1, \\,2}}\n", "+ c\\,s\\,{a_{2, \\,2}}\n", "$$\n", "$$\n", "- s\\,c\\,{a_{1, \\,1}} - s^{2}\\,{a_{1, \\,2}} + c^{2}\\,{a_{1, \\,2}}\n", "+ c\\,s\\,{a_{2, \\,2}}\n", "$$\n", "この非対角要素を0にする$\\theta$は以下のように求まる．\n", "\n", "|　　　　　　　　 |\n", "|:----|\n", "\n", "このとき注目している$i,j=1,2$以外の要素も変化する．\n", "```maple\n", ">expand(TT[3,1]);\n", "expand(TT[3,2]);\n", "```\n", "$$\n", "c\\,{a_{1, \\,3}} + s\\,{a_{2, \\,3}}\n", "$$\n", "$$\n", "- s\\,{a_{1, \\,3}} + c\\,{a_{2, \\,3}}\n", "$$\n", "これによって一旦0になった要素も値を持つが，なんども繰り返すことによって，徐々に0へ近づいていく．\n", "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## C 言語での実装\n", "\n", "lapack_codes/c_versions/Jacobi2.c にコードを入れている．\n", "```bash\n", "$ gcc Jacobi2.c\n", "$ a.out < input.txt\n", "```\n", "\n", "で固有値を求めていく過程が表示される．" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 数値計算ライブラリーについて\n", "\n", "- この辺りに記述しているC, Python, Maple, ruby言語のコードは，numerical_calc/lapack_codesにおいている\n", "\n", "\n", "一般の数値計算ライブラリーについては，時間の関係で講義ではその能力を紹介するにとどめる．昔の演習で詳しく取り上げていたので，研究や今後のために必要と思うときは，テキストを取りにおいで．\n", "\n", "行列の計算は，数値計算の中でも特に利用する機会が多く，また，律速ルーチンとなる可能性が高い．そこで，古くから行列計算の高速ルーチンが開発されてきた．なかでもBLASとLAPACKはフリーながら非常に高速である． \n", "\n", "前回に示した，逆行列を求める単純なLU分解法をC言語でコーディングしたものと，LAPACKのルーチンを比べた場合，1000次元の行列で計測すると\n", "```maple\n", "> 1000 [dim] 2.5200 [sec] #BOB\n", "> 1000 [dim] 0.4700 [sec] #LAPACK\n", "```\n", "となった．2006年に初めてこの計算に用いたPCはMacBook(2GHz, Intel Core Duo)であるが，\n", "この計算での0.47秒は1.4GFLOPに相当する．\n", "07年のMacBook(2GHz, Intel Core 2 Duo)ではさらに速くなって\n", "```maple\n", "bob% gcc -O3 bob.c -o bob\n", "bob% ./bob\n", "1000\n", " 1000 [dim] 1.7543 [sec] #BOB\n", "bob% gcc -O3 lapack.c -llapack -lblas -o lapack\n", "bob% ./lapack\n", "1000\n", " 1000 [dim] 0.1893 [sec] #LAPACK\n", "```\n", "で，3.5GFLOPSが出ている．今(2016年)は，MacBookAir(2.2GHz, Intel Core i7)で...\n", "top500.orgが毎年２回High Performance Computerのランクを発表している．\n", "今は，Top1は100PFlopsであるが，\n", "初回の1994年6月の500位は0.4GFlopsで，今のlaptopがはるかに凌いでいる．\n", "まさにlaptopスパコンの時代なんですよ．\n", "\n", "ライブラリーは世界中の計算機屋さんがよってたかって検証しているので，バグがほとんど無く，また，高速である．\n", "初学者はライブラリーを使うべきである．ただし，下のサンプルプログラムの行列生成の違いのように，ブラックボックス化すると思わぬ間違い（ここではFortranとCでの行列の並び順の違いが原因)をしでかすことがあるので，プログラムに組み込む前に必ず小さい次元(サンプルコード)で検証しておくこと\\footnote{少し前(2002年ごろ)GotoBLASが開発されて，性能が10%ほども上がった}．\n", "\n", "添付のコードはちょっと長いが時間があればフォローせよ．コンパイルは，OSXでは\n", "```maple\n", "> gcc -O3 -UPRINT lapack.c -llapack -lblas\n", "```\n", "とすればできる．linuxではLAPACK, BLASがインストールされていれば，\n", "```maple\n", "> #include \n", "```\n", "をコメントアウトして，\n", "```maple\n", "> gcc -O3 -DPRINT lapack.c -L/usr/local/lib64 -llapack -lblas -lg2c\n", "```\n", "などとすればコンパイルできるはず．\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## python謹製 lapack利用逆行列" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1000 [dim], 0.10614 [sec] # python\n", "2000 [dim], 0.50773 [sec] # python\n", "4000 [dim], 2.25623 [sec] # python\n" ] } ], "source": [ "import numpy as np\n", "import scipy.linalg as la\n", "import time\n", "\n", "sizes = [1000, 2000, 4000]\n", "for n in sizes:\n", " A = np.random.random((n, n))\n", " b = np.random.random((n))\n", " start = time.time()\n", " la.lapack.dgesv(A, b)\n", " print('%s [dim], %10.5f [sec] # python' % (n,time.time()-start))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "|size|elapsed_time[sec]|\n", "|:----|----:|\n", "|100 |0.004878|\n", "|1000 |0.085786|\n", "|10000|46.031718|" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 計算速度競争\n", "\n", "\n", "Mapleと，C, pythonで競争させてみました．Cはそのままです．Mapleは\n", "``` octave\n", "with(LinearAlgebra):\n", "data:=[1000,2000,4000];\n", "for n in data do\n", " A:=RandomMatrix(n,n,generator = 0..1.0):\n", " b:=RandomVector(n,generator = 0..1.0):\n", " st:=time();\n", " LUDecomposition():\n", " print(time()-st);\n", "end:\n", "```\n", "\n", "です．pythonのcodeは，lapackにあるdgesvを指定して\n", "呼び出すようにしています．scipy.laのsolveでは何を使っているかはよくわからないので．\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "結果は，次の通りです．\n", "\n", "|size|Maple[sec]|C[sec]|python[sec]|\n", "|:----|-----:|-----:|------:|\n", "|1000 | 0.339|0.0589|0.04101\n", "|2000 | 1.384|0.3826|0.12682\n", "|4000 |45.482|3.0111|0.81373" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "pythonの圧勝ですね．絶対値はあまり意味がありません．MacBook Air(13-inch, Early 2015/2.2GHz Core i7, ElCapitan 10.11.6)ですが，softのversionにもよります．\n", "\n", "ただし，どんなrandomを生成しているかで結果は大きく変わるので，ちょっと怪しいです．\n", "たとえば，MapleのRandomMatrixでgeneratorを指定しないと\n", "size=1000でも9.773[sec]ととんでもないぐらい時間がかかります．\n", "\n", "MapleではさらにRandomMatrixの生成にも時間がかかります．\n", "MapleではNAGのライブラリを使っているので．．．" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 課題\n", "\n", "1. 次の行列$A$　の固有値をpythonで求めよ．\n", "また，対角化行列$P$ を求めて，対角化せよ．\n", "$$ A\\, = \\, \\left[ \\begin {array}{ccc} 0&1&-2\\\\\n", "-3&7&-3 \\\\\n", "3&-5&5\n", "\\end {array} \\right]\n", "$$\n" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 1.000 9.000 2.000]\n", "[[-0.905 0.229 -0.707]\n", " [-0.302 0.688 0.000]\n", " [ 0.302 -0.688 0.707]]\n" ] } ], "source": [ "import numpy as np\n", "A = np.array([[0,1,-2],[-3,7,-3],[3,-5,5]])\n", "l, P = np.linalg.eig( A )\n", "print(l)\n", "print(P)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "'[[ 1. 0. -0.]\\n [-0. 9. -0.]\\n [ 0. 0. 2.]]'\n" ] } ], "source": [ "from pprint import pprint\n", "dA = np.dot(np.linalg.inv(P),np.dot(A,P))\n", "np.set_printoptions(precision=3, suppress=True)\n", "pprint(format(dA))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 行列のpprint設定イロイロ" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "'[[ 1. 0. -0.]\\n [-0. 9. -0.]\\n [ 0. 0. 2.]]'\n" ] } ], "source": [ "from pprint import pprint\n", "dA = np.dot(np.linalg.inv(P),np.dot(A,P))\n", "pprint(format(dA))" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "array([[ 1., 0., -0.],\n", " [-0., 9., -0.],\n", " [ 0., 0., 2.]])\n" ] } ], "source": [ "np.set_printoptions(precision=3)\n", "pprint(dA)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "array([[ 1., 0., -0.],\n", " [-0., 9., -0.],\n", " [ 0., 0., 2.]])\n" ] } ], "source": [ "np.set_printoptions(precision=3, suppress=True)\n", "pprint(dA)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "array([[ 1.000, 0.000, -0.000],\n", " [-0.000, 9.000, -0.000],\n", " [ 0.000, 0.000, 2.000]])\n" ] } ], "source": [ "np.set_printoptions(formatter={'float': '{: 0.3f}'.format}) \n", "pprint(dA)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "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.1" }, "latex_envs": { "LaTeX_envs_menu_present": true, "autocomplete": true, "bibliofile": "biblio.bib", "cite_by": "apalike", "current_citInitial": 1, "eqLabelWithNumbers": true, "eqNumInitial": 1, "hotkeys": { "equation": "Ctrl-E", "itemize": "Ctrl-I" }, "labels_anchors": false, "latex_user_defs": false, "report_style_numbering": false, "user_envs_cfg": false }, "toc": { "colors": { "hover_highlight": "#DAA520", "navigate_num": "#000000", "navigate_text": "#333333", "running_highlight": "#FF0000", "selected_highlight": "#FFD700", "sidebar_border": "#EEEEEE", "wrapper_background": "#FFFFFF" }, "moveMenuLeft": true, "nav_menu": { "height": "12px", "width": "252px" }, "navigate_menu": true, "number_sections": true, "sideBar": true, "threshold": 4, "toc_cell": true, "toc_section_display": "block", "toc_window_display": true, "widenNotebook": false } }, "nbformat": 4, "nbformat_minor": 2 }