Cartman0 · June 1, 2019 20:19
diff --git a/共役勾配CG法.ipynb b/共役勾配CG法.ipynb
 {
  "cells": [
    {
      "metadata": {
        "toc": "true"
      },
      "cell_type": "markdown",
      "source": "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#共役勾配法-Conjugate-Gradient-Method\" data-toc-modified-id=\"共役勾配法-Conjugate-Gradient-Method-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>共役勾配法 Conjugate Gradient Method</a></span><ul class=\"toc-item\"><li><span><a href=\"#共役\" data-toc-modified-id=\"共役-1.1\"><span class=\"toc-item-num\">1.1&nbsp;&nbsp;</span>共役</a></span><ul class=\"toc-item\"><li><span><a href=\"#コレスキー分解\" data-toc-modified-id=\"コレスキー分解-1.1.1\"><span class=\"toc-item-num\">1.1.1&nbsp;&nbsp;</span>コレスキー分解</a></span></li></ul></li><li><span><a href=\"#wikipedia解説\" data-toc-modified-id=\"wikipedia解説-1.2\"><span class=\"toc-item-num\">1.2&nbsp;&nbsp;</span>wikipedia解説</a></span><ul class=\"toc-item\"><li><span><a href=\"#直接法\" data-toc-modified-id=\"直接法-1.2.1\"><span class=\"toc-item-num\">1.2.1&nbsp;&nbsp;</span>直接法</a></span></li><li><span><a href=\"#反復法\" data-toc-modified-id=\"反復法-1.2.2\"><span class=\"toc-item-num\">1.2.2&nbsp;&nbsp;</span>反復法</a></span><ul class=\"toc-item\"><li><span><a href=\"#グラミシュミット法による共役ベクトルの作り方\" data-toc-modified-id=\"グラミシュミット法による共役ベクトルの作り方-1.2.2.1\"><span class=\"toc-item-num\">1.2.2.1&nbsp;&nbsp;</span>グラミシュミット法による共役ベクトルの作り方</a></span></li><li><span><a href=\"#CGにおける共役ベクトルの作成\" data-toc-modified-id=\"CGにおける共役ベクトルの作成-1.2.2.2\"><span class=\"toc-item-num\">1.2.2.2&nbsp;&nbsp;</span>CGにおける共役ベクトルの作成</a></span></li><li><span><a href=\"#各パラメータの求め方\" data-toc-modified-id=\"各パラメータの求め方-1.2.2.3\"><span class=\"toc-item-num\">1.2.2.3&nbsp;&nbsp;</span>各パラメータの求め方</a></span></li></ul></li></ul></li></ul></li></ul></div>"
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "# 共役勾配法 Conjugate Gradient Method"
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "資料：\n\n- https://ja.wikipedia.org/wiki/%E5%85%B1%E5%BD%B9%E5%8B%BE%E9%85%8D%E6%B3%95\n\n- [システムの最適化 － ２．非線形計画法（NP : Nonlinear Programming） －, ２．５. 共役勾配法（Conjugate Gradient Method）,sist.ac.jp](http://www.sist.ac.jp/~suganuma/kougi/other_lecture/SE/opt/nonlinear/nonlinear.htm#2.5)\n- [共役勾配法, 矢部博,  日本オペレーションズ・リサーチ学会, 1987年6月号](http://www.orsj.or.jp/~archive/pdf/bul/Vol.32_06_363.pdf])\n- NumericalRecipe 3rd edition, 10.8 Conjugate Gradient Methods in Multidimensions, p.539"
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "## 共役"
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "- $A$: 正定値対称行列\n- 非零ベクトル$\\vec{u}、\\vec{v}$\n\n$$\n\\mathbf {u} ^{ \\mathrm {T} }\\mathbf {A} \\mathbf {v} \n= 0\n$$\n\nを満たすとき、u、vはAに関して共役であるという。\n\n対称行列$A$を単位行列にすれば内積0，つまり直交を表す．つまり，共役は直交関係を拡張したもの\n\nAは対称正定値より，左辺を展開すると，\n\n$$\n\\vec{u} ^{ \\mathrm {T} }\\mathbf {A} \\vec{v} \n= \\vec{u}^{T} \\cdot \\mathbf{A}\\vec{v}\n= \\mathbf{A}^{T} \\vec{u} \\cdot \\vec{v}\n=  \\mathbf{A} \\vec{u} \\cdot \\vec{v}\n= 0\n$$\n\n片方を行列Aで変換したベクトルとの内積が0，つまり直交すれば共役となる．\n\n"
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "### コレスキー分解"
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "http://www.orsj.or.jp/~archive/pdf/bul/Vol.32_06_363.pdf\n\n\n対称行列$A$をコレスキー分解して\n\n$$\nA = P^{T}P\n$$\n\n$$\n\\vec{u}^{ \\mathrm {T} }\\mathbf {A} \\vec{v} \n= \\vec{u}^{\\mathrm {T}} P^{T} P \\vec{v}\n= (P\\vec{u})^{T}(P\\vec{v})\n= 0\n$$\n\n分解した行列で写像したベクトルが直交していることを表している．"
    },
    {
      "metadata": {
        "ExecuteTime": {
          "start_time": "2019-06-01T12:26:32.184466Z",
          "end_time": "2019-06-01T12:26:32.190415Z"
        },
        "trusted": true,
        "scrolled": true
      },
      "cell_type": "code",
      "source": "import scipy as sp\nfrom scipy import linalg\n\n# P = sp.matrix([[1, 0],\n#                [2, 1]])\n\n# P.T.dot(P)\nA = sp.matrix([[2, 1],\n               [1, 2]])\nP = linalg.cholesky(A)\nP",
      "execution_count": 9,
      "outputs": [
        {
          "output_type": "execute_result",
          "execution_count": 9,
          "data": {
            "text/plain": "array([[1.41421356, 0.70710678],\n       [0.        , 1.22474487]])"
          },
          "metadata": {}
        }
      ]
    },
    {
      "metadata": {
        "ExecuteTime": {
          "start_time": "2019-06-01T12:26:42.320751Z",
          "end_time": "2019-06-01T12:26:42.324510Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "P.T.dot(P)",
      "execution_count": 10,
      "outputs": [
        {
          "output_type": "execute_result",
          "execution_count": 10,
          "data": {
            "text/plain": "array([[2., 1.],\n       [1., 2.]])"
          },
          "metadata": {}
        }
      ]
    },
    {
      "metadata": {
        "collapsed": true,
        "trusted": false
      },
      "cell_type": "markdown",
      "source": "## wikipedia解説"
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "### 直接法"
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "$p_{k}$を$n$個の互いに共役なベクトル列とする。**$p_{k}$ は基底Rnを構成する？(ここの理論がよくわからない)**\n\n> 定数ベクトル$\\vec{b}$との直交基底を考えている？\n$A\\vec{x}^{*}$ より，定数ベクトルと直交したベクトルは共役関係にあるので．\n\n共役より，\n\n$$\n\\vec{p}_{i}^{T} A \\vec{p}_{j} = 0, (i \\neq j)\n$$\n\nので、$A\\vec{x} = \\vec{b}$の解$x^{*}$をこの基底で展開すると、線形結合で表されるので，\n\n$$\n\\vec{x}^{*} = \\sum_{i=1}^{n} \\alpha_{i} \\vec{p} _{i}\n$$\n\nと書ける。\n\nこの解を代入すると，\n\n$$\nA\\vec{x}^{*} = A \\sum_{i=1}^{n} \\alpha_{i} \\vec{p} _{i}\n= \\sum_{i=1}^{n} \\alpha_{i} A \\vec{p}_{i}\n= \\vec{b}\n$$\n\n左から任意の共役ベクトルをかけて\n\n$$\n\\vec{p}_{k}^{T} A\\vec{x}^{*} \n= \\sum_{i=1}^{n} \\alpha_{i} \\vec{p}_{k}^{T} A \\vec{p} _{i}\n= \\vec{p}_{k}^{T} \\vec{b}\n$$\n\n\n共役により\n\n$$\n\\vec{p}_{k}^{T} A \\vec{p} _{i} = 0, (i \\neq k)\n$$\n\n残るのは二次形式の項のみ\n\n$$\n\\vec{p}_{k}^{T} A \\vec{p} _{k} > 0, (i==k)\n$$\n\n$$\n\\vec{p}_{k}^{T} A\\vec{x}^{*} \n= \\alpha_{k} \\vec{p}_{k}^{T} A \\vec{p} _{k}\n= \\vec{p}_{k}^{T} \\vec{b}\n$$\n\n\n$$\n\\alpha_{k} = \\frac{\\vec{p}_{k}^{T} \\vec{b}}{ \\vec{p}_{k}^{T} A \\vec{p} _{k} }\n$$\n\n正定値より，分母$\\vec{p}_{k}^{T} A \\vec{p} _{k} > 0$\n\nよって，解は\n\n$$\n\\vec{x}^{*} = \\sum_{i=1}^{n} \\alpha_{i} \\vec{p} _{i}\n=  \\sum_{i=1}^{n} \\frac{\\vec{p}_{i}^{T} \\vec{b}}{ \\vec{p}_{i}^{T} A \\vec{p} _{i} } \\vec{p} _{i}\n$$\n\n以上から、Ax = bを解くための方法が得られる。すなわち、まずn個の共役な方向を見つけ、それから係数αkを計算すればよい。"
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "### 反復法"
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "共役なベクトル列pkを注意深く選ぶことにより、一部のベクトルからx*の良い近似を得られる可能性がある。そこで、共役勾配法を反復法として利用することを考える。こうすることで、nが非常に大きく、直接法では解くのに時間がかかりすぎるような問題にも適用することができる。\n\nx*の初期値をx0 = 0 とする。x*が二次形式\n\n\n目的関数をいかにする．\n\n$$\nf(\\vec{x}) = \\frac{1}{2} \\vec{x}^{T} A \\vec{x} - \\vec{b}^{T}\\vec{x}, \n$$\n\n微分すると，\n\n$$\n\\nabla f \n= \\frac{d f(\\vec{x})}{d\\vec{x}} \n= A \\vec{x} - \\vec{b}\n$$\n\n勾配が元の連立方程式になっている．\n最小化問題なので，この勾配の負の方向に最適化していく．\nまた，最適解に近づくと，勾配自体も小さくなっていく．\n\n$$\n- \\nabla f  \n=  \\vec{b} - A \\vec{x} \n$$\n\n最初の基底ベクトルp1を$\\vec{x} = \\vec{x}_{0}$でのfの勾配$A \\vec{x}_{0}-\\vec{b}=−\\vec{b}$となるように取る。つまり，$\\vec{x}_{0} = -\\nabla f = \\vec{b}$.\nこのとき、基底の他のベクトルは勾配に共役である。そこで、この方法を共役勾配法と呼ぶ。\n\n$\\vec{r_{k}} =\n= - \\nabla f = \\vec{b} - A \\vec{x_{k}}$ を勾配かつ残差ベクトルとし，\n\n定数ベクトル$\\vec{b}$は，反復回数を重ねるごとにA_{x_k}で補完されていくので勾配・残差ベクトルは小さくなっていく．\n\n係数は，直接法の式の定数ベクトルを$\\vec{b} = $\n\n$$\n\\vec{r_{k}} = \\vec{b} - A \\vec{x_{k}} \\\\\nA \\vec{x_{k}} = \\vec{b} - \\vec{r_{k}} \\\\\nA \\vec{x_{k}} \n= \\sum_{i=1}^{k}\\alpha_i A \\vec{p_i}\n= \\vec{b} - \\vec{r_{k}} \\\\\n$$\n\n"
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "#### グラミシュミット法による共役ベクトルの作り方"
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "グラミシュミットの正規直交基底の共役ver 内積の式を共役に変えればいい．\n\nhttp://www.orsj.or.jp/~archive/pdf/bul/Vol.32_06_363.pdf\n\n\n1次独立な2つのベクトル $\\vec{v}_{0}$ と$\\vec{v}_{1}$が与えら\nれている場合を考える.このとき$\\vec{v}_{1}$を $\\vec{v}_{0}$ 共役の意味で正射影して得られるベクトルは $\\gamma_{10} \\vec{v}_{0}$ と表わされるから $\\vec{d}_{1} = \\vec{v}_{1} - \\gamma_{10}\\vec{v}_{0}$ が $\\vec{d_{1}, \\vec{v_0}}_{A} = \\vec{d_{1}}^{T}A\\vec{v}_{0}  = 0$ を満たすように係\n$\\gamma_{10}$を決めれば\n\n$$\n\\vec{d_{1}} \n= \\vec{v}_{1} - (\\frac{(\\vec{v_0}, \\vec{v_1})_{A}}{(\\vec{v_0}, \\vec{v_0})_{A}})\\vec{v_{0}} \n= \n\\vec{v}_{1} - \\left(\\frac{\\vec{v_0}^{T} A \\vec{v_1}}{\\vec{v_0}^{T} A \\vec{v_0} }\\right) \\vec{v_{0}}\n$$"
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "#### CGにおける共役ベクトルの作成"
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "新しく作る$p_{k+1}$ベクトルは,今までの基底ベクトルとも共役関係になる．\n\n$$\n\\vec{p}_{k+1}^{T} A \\vec{p}_{k} = 0\n$$"
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "勾配ベクトル$\\vec{r_{k}}$ に近いほうに共役である基底ベクトルを作る．\n$\\vec{r_{k}}$　から $p_{k}$との基底をつくる．\n\n$$\n\\vec{p_{k+1}} \n= \n\\vec{r}_{k} - \\left(\\frac{\\vec{p_{k}}^{T} A \\vec{r}_{k}}{\\vec{p_{k}}^{T} A \\vec{p_{k}} }\\right) \\vec{p_{k}}\n$$\n"
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "#### 各パラメータの求め方"
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "誤差・勾配ベクトルは\n\n$$\n\\vec{b} - A\\vec{x}_{k} = \\vec{r}_{k} \\\\\nA\\vec{x}^{*} - A\\vec{x}_{k} = \\vec{r}_{k}\\\\\nA (\\vec{x}^{*} - \\vec{x}_{k}) = \\vec{r}_{k}\\\\\nA (\\sum_{i=1}^{n}\\alpha_{i}\\vec{p}_{i} - \\sum_{i=1}^{k}\\alpha_{i}\\vec{p}_{i}) = \\vec{r}_{k}\\\\\nA \\sum_{i=k+1}^{n}\\alpha_{i}\\vec{p}_{i} = \\vec{r}_{k}\\\\\n\\sum_{i=k+1}^{n}\\alpha_{i} A\\vec{p}_{i} = \\vec{r}_{k}\\\\\\\n\\sum_{i=k+1}^{n}\\alpha_{i} \\vec{p}_{k+1}^{T} A \\vec{p}_{i} =  \\vec{p}_{k+1}^{T} \\vec{r}_{k}\\\\\\\n$$\n\n共役により$\\vec{p}_{k+1}$以外は0\n\n$$\n\\alpha_{k+1} \\vec{p}_{k+1}^{T} A \\vec{p}_{k+1} =  \\vec{p}_{k+1}^{T} \\vec{r}_{k}\\\\\\\n\\alpha_{k+1} = \\frac{\\vec{p}_{k+1}^{T} \\vec{r}_{k}}{\\vec{p}_{k+1}^{T} A \\vec{p}_{k+1}}\n$$\n"
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "$$\n\\vec{x}_{k+1}\n= \\sum_{i=1}^{k+1} \\alpha_{i}\\vec{p}_{i}\n= \\sum_{i=1}^{k} \\alpha_{i}\\vec{p}_{i} + \\alpha_{k+1}\\vec{p}_{k+1}\n= \\vec{x}_{k} + \\alpha_{k+1}\\vec{p}_{k+1}\n$$"
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "$$\n\\vec{r}_{k+1} = A (\\sum_{i=1}^{n}\\alpha_{i}\\vec{p}_{i} - \\sum_{i=1}^{k+1}\\alpha_{i}\\vec{p}_{i}) \\\\\n = A \\sum_{i=(k+1)+1}^{n}\\alpha_{i}\\vec{p}_{i} \\\\\n = \\sum_{i=k+1}^{n}\\alpha_{i}A\\vec{p}_{i} - \\alpha_{k+1}A\\vec{p}_{k+1} \\\\\n = \\vec{r}_{k} - \\alpha_{k+1}A\\vec{p}_{k+1} \\\\\n$$"
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "",
      "execution_count": null,
      "outputs": []
    }
  ],
  "metadata": {
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3",
      "language": "python"
    },
    "language_info": {
      "name": "python",
      "version": "3.7.0",
      "mimetype": "text/x-python",
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "pygments_lexer": "ipython3",
      "nbconvert_exporter": "python",
      "file_extension": ".py"
    },
    "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": {
      "nav_menu": {},
      "number_sections": true,
      "sideBar": true,
      "skip_h1_title": false,
      "base_numbering": 1,
      "title_cell": "Table of Contents",
      "title_sidebar": "Contents",
      "toc_cell": true,
      "toc_position": {
        "height": "449px",
        "left": "0px",
        "right": "720px",
        "top": "161px",
        "width": "174px"
      },
      "toc_section_display": "block",
      "toc_window_display": true
    },
    "hide_input": false,
    "gist": {
      "id": "",
      "data": {
        "description": "Optimize 共役勾配CG法",
        "public": true
      }
    }
  },
  "nbformat": 4,
  "nbformat_minor": 2
 }
	{
	"cells": [
	{
	"metadata": {
	"toc": "true"
	},
	"cell_type": "markdown",
	"source": "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#共役勾配法-Conjugate-Gradient-Method\" data-toc-modified-id=\"共役勾配法-Conjugate-Gradient-Method-1\"><span class=\"toc-item-num\">1  </span>共役勾配法 Conjugate Gradient Method</a></span><ul class=\"toc-item\"><li><span><a href=\"#共役\" data-toc-modified-id=\"共役-1.1\"><span class=\"toc-item-num\">1.1  </span>共役</a></span><ul class=\"toc-item\"><li><span><a href=\"#コレスキー分解\" data-toc-modified-id=\"コレスキー分解-1.1.1\"><span class=\"toc-item-num\">1.1.1  </span>コレスキー分解</a></span></li></ul></li><li><span><a href=\"#wikipedia解説\" data-toc-modified-id=\"wikipedia解説-1.2\"><span class=\"toc-item-num\">1.2  </span>wikipedia解説</a></span><ul class=\"toc-item\"><li><span><a href=\"#直接法\" data-toc-modified-id=\"直接法-1.2.1\"><span class=\"toc-item-num\">1.2.1  </span>直接法</a></span></li><li><span><a href=\"#反復法\" data-toc-modified-id=\"反復法-1.2.2\"><span class=\"toc-item-num\">1.2.2  </span>反復法</a></span><ul class=\"toc-item\"><li><span><a href=\"#グラミシュミット法による共役ベクトルの作り方\" data-toc-modified-id=\"グラミシュミット法による共役ベクトルの作り方-1.2.2.1\"><span class=\"toc-item-num\">1.2.2.1  </span>グラミシュミット法による共役ベクトルの作り方</a></span></li><li><span><a href=\"#CGにおける共役ベクトルの作成\" data-toc-modified-id=\"CGにおける共役ベクトルの作成-1.2.2.2\"><span class=\"toc-item-num\">1.2.2.2  </span>CGにおける共役ベクトルの作成</a></span></li><li><span><a href=\"#各パラメータの求め方\" data-toc-modified-id=\"各パラメータの求め方-1.2.2.3\"><span class=\"toc-item-num\">1.2.2.3  </span>各パラメータの求め方</a></span></li></ul></li></ul></li></ul></li></ul></div>"
	},
	{
	"metadata": {},
	"cell_type": "markdown",
	"source": "# 共役勾配法 Conjugate Gradient Method"
	},
	{
	"metadata": {},
	"cell_type": "markdown",
	"source": "資料：\n\n- https://ja.wikipedia.org/wiki/%E5%85%B1%E5%BD%B9%E5%8B%BE%E9%85%8D%E6%B3%95\n\n- [システムの最適化－２．非線形計画法（NP : Nonlinear Programming）－, ２．５. 共役勾配法（Conjugate Gradient Method）,sist.ac.jp](http://www.sist.ac.jp/~suganuma/kougi/other_lecture/SE/opt/nonlinear/nonlinear.htm#2.5)\n- [共役勾配法, 矢部博, 日本オペレーションズ・リサーチ学会, 1987年6月号](http://www.orsj.or.jp/~archive/pdf/bul/Vol.32_06_363.pdf])\n- NumericalRecipe 3rd edition, 10.8 Conjugate Gradient Methods in Multidimensions, p.539"
	},
	{
	"metadata": {},
	"cell_type": "markdown",
	"source": "## 共役"
	},
	{
	"metadata": {},
	"cell_type": "markdown",
	"source": "- $A$: 正定値対称行列\n- 非零ベクトル$\\vec{u}、\\vec{v}$\n\n$$\n\\mathbf {u} ^{ \\mathrm {T} }\\mathbf {A} \\mathbf {v} \n= 0\n$$\n\nを満たすとき、u、vはAに関して共役であるという。\n\n対称行列$A$を単位行列にすれば内積0，つまり直交を表す．つまり，共役は直交関係を拡張したもの\n\nAは対称正定値より，左辺を展開すると，\n\n$$\n\\vec{u} ^{ \\mathrm {T} }\\mathbf {A} \\vec{v} \n= \\vec{u}^{T} \\cdot \\mathbf{A}\\vec{v}\n= \\mathbf{A}^{T} \\vec{u} \\cdot \\vec{v}\n= \\mathbf{A} \\vec{u} \\cdot \\vec{v}\n= 0\n$$\n\n片方を行列Aで変換したベクトルとの内積が0，つまり直交すれば共役となる．\n\n"
	},
	{
	"metadata": {},
	"cell_type": "markdown",
	"source": "### コレスキー分解"
	},
	{
	"metadata": {},
	"cell_type": "markdown",
	"source": "http://www.orsj.or.jp/~archive/pdf/bul/Vol.32_06_363.pdf\n\n\n対称行列$A$をコレスキー分解して\n\n$$\nA = P^{T}P\n$$\n\n$$\n\\vec{u}^{ \\mathrm {T} }\\mathbf {A} \\vec{v} \n= \\vec{u}^{\\mathrm {T}} P^{T} P \\vec{v}\n= (P\\vec{u})^{T}(P\\vec{v})\n= 0\n$$\n\n分解した行列で写像したベクトルが直交していることを表している．"
	},
	{
	"metadata": {
	"ExecuteTime": {
	"start_time": "2019-06-01T12:26:32.184466Z",
	"end_time": "2019-06-01T12:26:32.190415Z"
	},
	"trusted": true,
	"scrolled": true
	},
	"cell_type": "code",
	"source": "import scipy as sp\nfrom scipy import linalg\n\n# P = sp.matrix([[1, 0],\n# [2, 1]])\n\n# P.T.dot(P)\nA = sp.matrix([[2, 1],\n [1, 2]])\nP = linalg.cholesky(A)\nP",
	"execution_count": 9,
	"outputs": [
	{
	"output_type": "execute_result",
	"execution_count": 9,
	"data": {
	"text/plain": "array([[1.41421356, 0.70710678],\n [0. , 1.22474487]])"
	},
	"metadata": {}
	}
	]
	},
	{
	"metadata": {
	"ExecuteTime": {
	"start_time": "2019-06-01T12:26:42.320751Z",
	"end_time": "2019-06-01T12:26:42.324510Z"
	},
	"trusted": true
	},
	"cell_type": "code",
	"source": "P.T.dot(P)",
	"execution_count": 10,
	"outputs": [
	{
	"output_type": "execute_result",
	"execution_count": 10,
	"data": {
	"text/plain": "array([[2., 1.],\n [1., 2.]])"
	},
	"metadata": {}
	}
	]
	},
	{
	"metadata": {
	"collapsed": true,
	"trusted": false
	},
	"cell_type": "markdown",
	"source": "## wikipedia解説"
	},
	{
	"metadata": {},
	"cell_type": "markdown",
	"source": "### 直接法"
	},
	{
	"metadata": {},
	"cell_type": "markdown",
	"source": "$p_{k}$を$n$個の互いに共役なベクトル列とする。$p_{k}$ は基底Rnを構成する？(ここの理論がよくわからない)\n\n> 定数ベクトル$\\vec{b}$との直交基底を考えている？\n$A\\vec{x}^{}$ より，定数ベクトルと直交したベクトルは共役関係にあるので．\n\n共役より，\n\n$$\n\\vec{p}_{i}^{T} A \\vec{p}_{j} = 0, (i \\neq j)\n$$\n\nので、$A\\vec{x} = \\vec{b}$の解$x^{}$をこの基底で展開すると、線形結合で表されるので，\n\n$$\n\\vec{x}^{} = \\sum_{i=1}^{n} \\alpha_{i} \\vec{p} _{i}\n$$\n\nと書ける。\n\nこの解を代入すると，\n\n$$\nA\\vec{x}^{} = A \\sum_{i=1}^{n} \\alpha_{i} \\vec{p} _{i}\n= \\sum_{i=1}^{n} \\alpha_{i} A \\vec{p}_{i}\n= \\vec{b}\n$$\n\n左から任意の共役ベクトルをかけて\n\n$$\n\\vec{p}_{k}^{T} A\\vec{x}^{} \n= \\sum_{i=1}^{n} \\alpha_{i} \\vec{p}_{k}^{T} A \\vec{p} _{i}\n= \\vec{p}_{k}^{T} \\vec{b}\n$$\n\n\n共役により\n\n$$\n\\vec{p}_{k}^{T} A \\vec{p} _{i} = 0, (i \\neq k)\n$$\n\n残るのは二次形式の項のみ\n\n$$\n\\vec{p}_{k}^{T} A \\vec{p} _{k} > 0, (i==k)\n$$\n\n$$\n\\vec{p}_{k}^{T} A\\vec{x}^{} \n= \\alpha_{k} \\vec{p}_{k}^{T} A \\vec{p} _{k}\n= \\vec{p}_{k}^{T} \\vec{b}\n$$\n\n\n$$\n\\alpha_{k} = \\frac{\\vec{p}_{k}^{T} \\vec{b}}{ \\vec{p}_{k}^{T} A \\vec{p} _{k} }\n$$\n\n正定値より，分母$\\vec{p}_{k}^{T} A \\vec{p} _{k} > 0$\n\nよって，解は\n\n$$\n\\vec{x}^{*} = \\sum_{i=1}^{n} \\alpha_{i} \\vec{p} _{i}\n= \\sum_{i=1}^{n} \\frac{\\vec{p}_{i}^{T} \\vec{b}}{ \\vec{p}_{i}^{T} A \\vec{p} _{i} } \\vec{p} _{i}\n$$\n\n以上から、Ax = bを解くための方法が得られる。すなわち、まずn個の共役な方向を見つけ、それから係数αkを計算すればよい。"
	},
	{
	"metadata": {},
	"cell_type": "markdown",
	"source": "### 反復法"
	},
	{
	"metadata": {},
	"cell_type": "markdown",
	"source": "共役なベクトル列pkを注意深く選ぶことにより、一部のベクトルからxの良い近似を得られる可能性がある。そこで、共役勾配法を反復法として利用することを考える。こうすることで、nが非常に大きく、直接法では解くのに時間がかかりすぎるような問題にも適用することができる。\n\nxの初期値をx0 = 0 とする。x*が二次形式\n\n\n目的関数をいかにする．\n\n$$\nf(\\vec{x}) = \\frac{1}{2} \\vec{x}^{T} A \\vec{x} - \\vec{b}^{T}\\vec{x}, \n$$\n\n微分すると，\n\n$$\n\\nabla f \n= \\frac{d f(\\vec{x})}{d\\vec{x}} \n= A \\vec{x} - \\vec{b}\n$$\n\n勾配が元の連立方程式になっている．\n最小化問題なので，この勾配の負の方向に最適化していく．\nまた，最適解に近づくと，勾配自体も小さくなっていく．\n\n$$\n- \\nabla f \n= \\vec{b} - A \\vec{x} \n$$\n\n最初の基底ベクトルp1を$\\vec{x} = \\vec{x}_{0}$でのfの勾配$A \\vec{x}_{0}-\\vec{b}=−\\vec{b}$となるように取る。つまり，$\\vec{x}_{0} = -\\nabla f = \\vec{b}$.\nこのとき、基底の他のベクトルは勾配に共役である。そこで、この方法を共役勾配法と呼ぶ。\n\n$\\vec{r_{k}} =\n= - \\nabla f = \\vec{b} - A \\vec{x_{k}}$ を勾配かつ残差ベクトルとし，\n\n定数ベクトル$\\vec{b}$は，反復回数を重ねるごとにA_{x_k}で補完されていくので勾配・残差ベクトルは小さくなっていく．\n\n係数は，直接法の式の定数ベクトルを$\\vec{b} = $\n\n$$\n\\vec{r_{k}} = \\vec{b} - A \\vec{x_{k}} \\\\\nA \\vec{x_{k}} = \\vec{b} - \\vec{r_{k}} \\\\\nA \\vec{x_{k}} \n= \\sum_{i=1}^{k}\\alpha_i A \\vec{p_i}\n= \\vec{b} - \\vec{r_{k}} \\\\\n$$\n\n"
	},
	{
	"metadata": {},
	"cell_type": "markdown",
	"source": "#### グラミシュミット法による共役ベクトルの作り方"
	},
	{
	"metadata": {},
	"cell_type": "markdown",
	"source": "グラミシュミットの正規直交基底の共役ver 内積の式を共役に変えればいい．\n\nhttp://www.orsj.or.jp/~archive/pdf/bul/Vol.32_06_363.pdf\n\n\n1次独立な2つのベクトル $\\vec{v}_{0}$ と$\\vec{v}_{1}$が与えら\nれている場合を考える.このとき$\\vec{v}_{1}$を $\\vec{v}_{0}$ 共役の意味で正射影して得られるベクトルは $\\gamma_{10} \\vec{v}_{0}$ と表わされるから $\\vec{d}_{1} = \\vec{v}_{1} - \\gamma_{10}\\vec{v}_{0}$ が $\\vec{d_{1}, \\vec{v_0}}_{A} = \\vec{d_{1}}^{T}A\\vec{v}_{0} = 0$ を満たすように係\n$\\gamma_{10}$を決めれば\n\n$$\n\\vec{d_{1}} \n= \\vec{v}_{1} - (\\frac{(\\vec{v_0}, \\vec{v_1})_{A}}{(\\vec{v_0}, \\vec{v_0})_{A}})\\vec{v_{0}} \n= \n\\vec{v}_{1} - \\left(\\frac{\\vec{v_0}^{T} A \\vec{v_1}}{\\vec{v_0}^{T} A \\vec{v_0} }\\right) \\vec{v_{0}}\n$$"
	},
	{
	"metadata": {},
	"cell_type": "markdown",
	"source": "#### CGにおける共役ベクトルの作成"
	},
	{
	"metadata": {},
	"cell_type": "markdown",
	"source": "新しく作る$p_{k+1}$ベクトルは,今までの基底ベクトルとも共役関係になる．\n\n$$\n\\vec{p}_{k+1}^{T} A \\vec{p}_{k} = 0\n$$"
	},
	{
	"metadata": {},
	"cell_type": "markdown",
	"source": "勾配ベクトル$\\vec{r_{k}}$ に近いほうに共役である基底ベクトルを作る．\n$\\vec{r_{k}}$　から $p_{k}$との基底をつくる．\n\n$$\n\\vec{p_{k+1}} \n= \n\\vec{r}_{k} - \\left(\\frac{\\vec{p_{k}}^{T} A \\vec{r}_{k}}{\\vec{p_{k}}^{T} A \\vec{p_{k}} }\\right) \\vec{p_{k}}\n$$\n"
	},
	{
	"metadata": {},
	"cell_type": "markdown",
	"source": "#### 各パラメータの求め方"
	},
	{
	"metadata": {},
	"cell_type": "markdown",
	"source": "誤差・勾配ベクトルは\n\n$$\n\\vec{b} - A\\vec{x}_{k} = \\vec{r}_{k} \\\\\nA\\vec{x}^{} - A\\vec{x}_{k} = \\vec{r}_{k}\\\\\nA (\\vec{x}^{} - \\vec{x}_{k}) = \\vec{r}_{k}\\\\\nA (\\sum_{i=1}^{n}\\alpha_{i}\\vec{p}_{i} - \\sum_{i=1}^{k}\\alpha_{i}\\vec{p}_{i}) = \\vec{r}_{k}\\\\\nA \\sum_{i=k+1}^{n}\\alpha_{i}\\vec{p}_{i} = \\vec{r}_{k}\\\\\n\\sum_{i=k+1}^{n}\\alpha_{i} A\\vec{p}_{i} = \\vec{r}_{k}\\\\\\\n\\sum_{i=k+1}^{n}\\alpha_{i} \\vec{p}_{k+1}^{T} A \\vec{p}_{i} = \\vec{p}_{k+1}^{T} \\vec{r}_{k}\\\\\\\n$$\n\n共役により$\\vec{p}_{k+1}$以外は0\n\n$$\n\\alpha_{k+1} \\vec{p}_{k+1}^{T} A \\vec{p}_{k+1} = \\vec{p}_{k+1}^{T} \\vec{r}_{k}\\\\\\\n\\alpha_{k+1} = \\frac{\\vec{p}_{k+1}^{T} \\vec{r}_{k}}{\\vec{p}_{k+1}^{T} A \\vec{p}_{k+1}}\n$$\n"
	},
	{
	"metadata": {},
	"cell_type": "markdown",
	"source": "$$\n\\vec{x}_{k+1}\n= \\sum_{i=1}^{k+1} \\alpha_{i}\\vec{p}_{i}\n= \\sum_{i=1}^{k} \\alpha_{i}\\vec{p}_{i} + \\alpha_{k+1}\\vec{p}_{k+1}\n= \\vec{x}_{k} + \\alpha_{k+1}\\vec{p}_{k+1}\n$$"
	},
	{
	"metadata": {},
	"cell_type": "markdown",
	"source": "$$\n\\vec{r}_{k+1} = A (\\sum_{i=1}^{n}\\alpha_{i}\\vec{p}_{i} - \\sum_{i=1}^{k+1}\\alpha_{i}\\vec{p}_{i}) \\\\\n = A \\sum_{i=(k+1)+1}^{n}\\alpha_{i}\\vec{p}_{i} \\\\\n = \\sum_{i=k+1}^{n}\\alpha_{i}A\\vec{p}_{i} - \\alpha_{k+1}A\\vec{p}_{k+1} \\\\\n = \\vec{r}_{k} - \\alpha_{k+1}A\\vec{p}_{k+1} \\\\\n$$"
	},
	{
	"metadata": {
	"trusted": true
	},
	"cell_type": "code",
	"source": "",
	"execution_count": null,
	"outputs": []
	}
	],
	"metadata": {
	"kernelspec": {
	"name": "python3",
	"display_name": "Python 3",
	"language": "python"
	},
	"language_info": {
	"name": "python",
	"version": "3.7.0",
	"mimetype": "text/x-python",
	"codemirror_mode": {
	"name": "ipython",
	"version": 3
	},
	"pygments_lexer": "ipython3",
	"nbconvert_exporter": "python",
	"file_extension": ".py"
	},
	"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": {
	"nav_menu": {},
	"number_sections": true,
	"sideBar": true,
	"skip_h1_title": false,
	"base_numbering": 1,
	"title_cell": "Table of Contents",
	"title_sidebar": "Contents",
	"toc_cell": true,
	"toc_position": {
	"height": "449px",
	"left": "0px",
	"right": "720px",
	"top": "161px",
	"width": "174px"
	},
	"toc_section_display": "block",
	"toc_window_display": true
	},
	"hide_input": false,
	"gist": {
	"id": "",
	"data": {
	"description": "Optimize 共役勾配CG法",
	"public": true
	}
	}
	},
	"nbformat": 4,
	"nbformat_minor": 2
	}
No results found