Conjugate_Gradient.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose we want to solve the system of linear equations\n",
    "\n",
    "$${\\displaystyle \\mathbf {A} \\mathbf {x} =\\mathbf {b} }$$\n",
    "for the vector $x$, where the known $n × n$ matrix $A$ is symmetric (i.e., $A^T = A$), positive-definite (i.e. $xTAx > 0$ for all non-zero vectors $x$ in $R^n$), and real, and $b$ is known as well. We denote the unique solution of this system by $x$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.sparse.linalg import cg\n",
    "\n",
    "from sklearn.datasets import make_sparse_spd_matrix as mspdm\n",
    "#from scipy.sparse.linalg import spsolve\n",
    "from numpy.linalg import solve, norm\n",
    "#import scipy.stats as stats\n",
    "#import scipy.sparse as sparse\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "#from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "np.random.seed(210)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here i use __make_sparse_spd_matrix__ function from __sklearn.datasets__ in order to generate  $n × n$ matrix $A$ which must be _sparse symmetric positive-definite and real_. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# def sprandsym(n, density):\n",
    "#     rvs = stats.norm().rvs\n",
    "#     X = sparse.random(n, n, density=density, data_rvs=rvs)\n",
    "#     upper_X = sparse.triu(X) \n",
    "#     result = upper_X + upper_X.T - sparse.diags(X.diagonal())\n",
    "#     return result\n",
    "\n",
    "# A = sprandsym(size, 0.01)\n",
    "# #print(repr(A))\n",
    "# # <5000x5000 sparse matrix of type '<class 'numpy.float64'>'\n",
    "# #   with 249909 stored elements in Compressed Sparse Row format>\n",
    "\n",
    "# A = A.toarray()\n",
    "# #A\n",
    "\n",
    "def gen_matrix(size):\n",
    "    A = mspdm(size, 0.99)\n",
    "    return(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#checking if matrix A satisfy all requirements\n",
    "def check_matrix(A):\n",
    "    #Check A for being Positive Definite\n",
    "    np.linalg.cholesky(A)\n",
    "\n",
    "    # # check that the matrix A is symmetric. The difference should have no non-zero elements\n",
    "    #A = A.toarray()\n",
    "    #assert (A - A.T).nnz == 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here i generate vectors $x$ and $b$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#generating x and b vectors\n",
    "def gen_xb(size):\n",
    "    b = np.random.rand(size,1)\n",
    "    #np.dtype(b) #why error?\n",
    "    x = (np.zeros(size)+1.0).reshape(-1, 1)\n",
    "    return x, b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here i generate a Jacobi preconditioner $n × n$ matrix $M$: $M =$ diag$(A)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#generating Jacobi preconditioner\n",
    "def gen_precond(A, size):\n",
    "    M = (np.zeros(size*size)).reshape(-1, size)\n",
    "    i = 0\n",
    "    j = 0\n",
    "    while i < size:\n",
    "        while j < size:\n",
    "            if(j==i):\n",
    "                M[i][j]=A[i][j]\n",
    "            j=j+1\n",
    "        i=i+1\n",
    "        j=0\n",
    "    return(M)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simple examples:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From Wiki article about Conjugate Gradient method https://en.wikipedia.org/wiki/Conjugate_gradient_method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#From Wiki\n",
    "A_e1 = np.array([[4,1],[1,3]])\n",
    "b_e1 = np.array([1,2]).reshape(-1, 1)\n",
    "x_e1 = np.array([2,1]).reshape(-1, 1)\n",
    "\n",
    "#%timeit cg(A_e1,b_e1)\n",
    "#print(cg(A_e1,b_e1)[0]) #correct answer\n",
    "#[0.09090909 0.63636364]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From Амосов А.А., Дубинский Ю.А., Копченова Н.В. Вычислительные методы для инженеров. http://pnu.edu.ru/media/filer_public/2013/04/10/6-1_amosov-dubinsky-kopchenova_1994.pdf p284."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#From Amosov\n",
    "A_e2 = np.array([[2,0],[0,4]])\n",
    "b_e2 = np.array([4,4]).reshape(-1, 1)\n",
    "x_e2 = np.array([1,1]).reshape(-1, 1)\n",
    "\n",
    "#%timeit cg(A_e2,b_e2)\n",
    "#print(cg(A_e2,b_e2)[0]) #correct answer\n",
    "#[2. 1.]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "My realization of the method. Reworked code Wiki article about Conjugate Gradient method https://en.wikipedia.org/wiki/Conjugate_gradient_method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "k - step\n",
    "r - residual of the k step\n",
    "p - conjugate vector\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#conjugate gradients method realisation\n",
    "def my_cg(A,b,x):\n",
    "    r = b - A @ x # the first residual\n",
    "    p = r # the first basis vectop p\n",
    "    rsold = r.T @ r\n",
    "    k = 0\n",
    "    while k < len(A)*10: # should converge in n steps (but it does not :( )\n",
    "        Ap = A @ p\n",
    "        alpha = rsold / (p.T @ Ap)\n",
    "        x = x + alpha * p\n",
    "        r = r - alpha * Ap\n",
    "        rsnew = r.T @ r\n",
    "        if norm(r) < 1e-10:\n",
    "            print(\"1e-10 Broken\")\n",
    "            break\n",
    "        p = r + (rsnew / rsold) * p\n",
    "        rsold = rsnew\n",
    "\n",
    "        k = k + 1\n",
    "    print(\"Counted in:\", k, \"steps\")\n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "My realization of the method. Reworked code Wiki article about Conjugate Gradient method https://en.wikipedia.org/wiki/Conjugate_gradient_method.<br>\n",
    "Also used this: An Introduction to the Conjugate Gradient Method Without the Agonizing Pain http://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf p39 to learn about preconditioning."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#conjugate gradients method realisation with Jacobi preconditioner\n",
    "def my_cg_prec(A,b,x,M):\n",
    "    M_inv = np.linalg.inv(M)\n",
    "    r = b - A @ x\n",
    "    z = M_inv @ r\n",
    "    p = z\n",
    "    i = 0\n",
    "    while i < len(A)*10: # should converge in n steps (but it does not :( )\n",
    "        Ap = A @ p\n",
    "        alpha = (r.T @ z) / (p.T @ Ap)\n",
    "        x = x + alpha * p\n",
    "        ztr_old = z.T @ r\n",
    "        r = r - alpha * Ap\n",
    "        if norm(r) < 1e-10:\n",
    "            print(\"1e-10 Broken\")\n",
    "            break\n",
    "        z = M_inv @ r\n",
    "        betta = (z.T @ r)/(ztr_old)\n",
    "        p = z + betta * p\n",
    "\n",
    "        i = i+1\n",
    "    print(\"Counted in:\", i, \"steps\")\n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Testing my code on the known problems(from Wiki and Amosov)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "#my_cg(A_e1,b_e1,x_e1) #correct"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "#my_cg(A_e2,b_e2,x_e2) #correct"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1e-10 Broken\n",
      "Counted in: 35 steps\n",
      "1e-10 Broken\n",
      "Counted in: 33 steps\n",
      "**********NORMAL**********\n",
      "1 LOOP COMPLETED size = 100\n",
      "mean error = [2.8134203e-06]\n",
      "**********PRECONDITIONED**********\n",
      "1 LOOP COMPLETED size = 100\n",
      "mean error = [2.81342086e-06]\n",
      "\n",
      "1e-10 Broken\n",
      "Counted in: 64 steps\n",
      "1e-10 Broken\n",
      "Counted in: 61 steps\n",
      "**********NORMAL**********\n",
      "2 LOOP COMPLETED size = 200\n",
      "mean error = [7.37651602e-06]\n",
      "**********PRECONDITIONED**********\n",
      "2 LOOP COMPLETED size = 200\n",
      "mean error = [7.37651509e-06]\n",
      "\n",
      "1e-10 Broken\n",
      "Counted in: 101 steps\n",
      "1e-10 Broken\n",
      "Counted in: 94 steps\n",
      "**********NORMAL**********\n",
      "3 LOOP COMPLETED size = 300\n",
      "mean error = [1.10326547e-05]\n",
      "**********PRECONDITIONED**********\n",
      "3 LOOP COMPLETED size = 300\n",
      "mean error = [1.10326539e-05]\n",
      "\n",
      "1e-10 Broken\n",
      "Counted in: 157 steps\n",
      "1e-10 Broken\n",
      "Counted in: 142 steps\n",
      "**********NORMAL**********\n",
      "4 LOOP COMPLETED size = 400\n",
      "mean error = [1.0635308e-05]\n",
      "**********PRECONDITIONED**********\n",
      "4 LOOP COMPLETED size = 400\n",
      "mean error = [1.06353041e-05]\n",
      "\n",
      "1e-10 Broken\n",
      "Counted in: 206 steps\n",
      "1e-10 Broken\n",
      "Counted in: 177 steps\n",
      "**********NORMAL**********\n",
      "5 LOOP COMPLETED size = 500\n",
      "mean error = [2.86935522e-05]\n",
      "**********PRECONDITIONED**********\n",
      "5 LOOP COMPLETED size = 500\n",
      "mean error = [2.86935502e-05]\n",
      "\n",
      "1e-10 Broken\n",
      "Counted in: 279 steps\n",
      "1e-10 Broken\n",
      "Counted in: 237 steps\n",
      "**********NORMAL**********\n",
      "6 LOOP COMPLETED size = 600\n",
      "mean error = [3.17296767e-05]\n",
      "**********PRECONDITIONED**********\n",
      "6 LOOP COMPLETED size = 600\n",
      "mean error = [3.17296779e-05]\n",
      "\n",
      "1e-10 Broken\n",
      "Counted in: 375 steps\n",
      "1e-10 Broken\n",
      "Counted in: 306 steps\n",
      "**********NORMAL**********\n",
      "7 LOOP COMPLETED size = 700\n",
      "mean error = [3.8268308e-05]\n",
      "**********PRECONDITIONED**********\n",
      "7 LOOP COMPLETED size = 700\n",
      "mean error = [3.82683058e-05]\n",
      "\n",
      "1e-10 Broken\n",
      "Counted in: 442 steps\n",
      "1e-10 Broken\n",
      "Counted in: 359 steps\n",
      "**********NORMAL**********\n",
      "8 LOOP COMPLETED size = 800\n",
      "mean error = [3.61678364e-05]\n",
      "**********PRECONDITIONED**********\n",
      "8 LOOP COMPLETED size = 800\n",
      "mean error = [3.61678356e-05]\n",
      "\n",
      "1e-10 Broken\n",
      "Counted in: 503 steps\n",
      "1e-10 Broken\n",
      "Counted in: 411 steps\n",
      "**********NORMAL**********\n",
      "9 LOOP COMPLETED size = 900\n",
      "mean error = [4.45296137e-05]\n",
      "**********PRECONDITIONED**********\n",
      "9 LOOP COMPLETED size = 900\n",
      "mean error = [4.45296198e-05]\n",
      "\n",
      "1e-10 Broken\n",
      "Counted in: 709 steps\n",
      "1e-10 Broken\n",
      "Counted in: 546 steps\n",
      "**********NORMAL**********\n",
      "10 LOOP COMPLETED size = 1000\n",
      "mean error = [5.22395686e-05]\n",
      "**********PRECONDITIONED**********\n",
      "10 LOOP COMPLETED size = 1000\n",
      "mean error = [5.22395648e-05]\n",
      "\n",
      "1e-10 Broken\n",
      "Counted in: 905 steps\n",
      "1e-10 Broken\n",
      "Counted in: 694 steps\n",
      "**********NORMAL**********\n",
      "11 LOOP COMPLETED size = 1100\n",
      "mean error = [8.13513698e-05]\n",
      "**********PRECONDITIONED**********\n",
      "11 LOOP COMPLETED size = 1100\n",
      "mean error = [8.1351388e-05]\n",
      "\n",
      "1e-10 Broken\n",
      "Counted in: 1210 steps\n",
      "1e-10 Broken\n",
      "Counted in: 889 steps\n",
      "**********NORMAL**********\n",
      "12 LOOP COMPLETED size = 1200\n",
      "mean error = [8.45014517e-05]\n",
      "**********PRECONDITIONED**********\n",
      "12 LOOP COMPLETED size = 1200\n",
      "mean error = [8.450128e-05]\n",
      "\n",
      "1e-10 Broken\n",
      "Counted in: 1544 steps\n",
      "1e-10 Broken\n",
      "Counted in: 1127 steps\n",
      "**********NORMAL**********\n",
      "13 LOOP COMPLETED size = 1300\n",
      "mean error = [0.0001437]\n",
      "**********PRECONDITIONED**********\n",
      "13 LOOP COMPLETED size = 1300\n",
      "mean error = [0.0001437]\n",
      "\n",
      "1e-10 Broken\n",
      "Counted in: 1851 steps\n",
      "1e-10 Broken\n",
      "Counted in: 1357 steps\n",
      "**********NORMAL**********\n",
      "14 LOOP COMPLETED size = 1400\n",
      "mean error = [0.0001701]\n",
      "**********PRECONDITIONED**********\n",
      "14 LOOP COMPLETED size = 1400\n",
      "mean error = [0.0001701]\n",
      "\n",
      "1e-10 Broken\n",
      "Counted in: 2386 steps\n",
      "1e-10 Broken\n",
      "Counted in: 1706 steps\n",
      "**********NORMAL**********\n",
      "15 LOOP COMPLETED size = 1500\n",
      "mean error = [0.00022801]\n",
      "**********PRECONDITIONED**********\n",
      "15 LOOP COMPLETED size = 1500\n",
      "mean error = [0.00022811]\n",
      "\n",
      "1e-10 Broken\n",
      "Counted in: 3224 steps\n",
      "1e-10 Broken\n",
      "Counted in: 2252 steps\n",
      "**********NORMAL**********\n",
      "16 LOOP COMPLETED size = 1600\n",
      "mean error = [0.00027489]\n",
      "**********PRECONDITIONED**********\n",
      "16 LOOP COMPLETED size = 1600\n",
      "mean error = [0.00027487]\n",
      "\n",
      "1e-10 Broken\n",
      "Counted in: 3451 steps\n",
      "1e-10 Broken\n",
      "Counted in: 2943 steps\n",
      "**********NORMAL**********\n",
      "17 LOOP COMPLETED size = 1700\n",
      "mean error = [0.00027149]\n",
      "**********PRECONDITIONED**********\n",
      "17 LOOP COMPLETED size = 1700\n",
      "mean error = [0.00027218]\n",
      "\n",
      "1e-10 Broken\n",
      "Counted in: 4692 steps\n",
      "1e-10 Broken\n",
      "Counted in: 3258 steps\n",
      "**********NORMAL**********\n",
      "18 LOOP COMPLETED size = 1800\n",
      "mean error = [0.00035487]\n",
      "**********PRECONDITIONED**********\n",
      "18 LOOP COMPLETED size = 1800\n",
      "mean error = [0.00035581]\n",
      "\n",
      "1e-10 Broken\n",
      "Counted in: 6208 steps\n",
      "1e-10 Broken\n",
      "Counted in: 4319 steps\n",
      "**********NORMAL**********\n",
      "19 LOOP COMPLETED size = 1900\n",
      "mean error = [0.00051525]\n",
      "**********PRECONDITIONED**********\n",
      "19 LOOP COMPLETED size = 1900\n",
      "mean error = [0.00152907]\n",
      "\n",
      "1e-10 Broken\n",
      "Counted in: 7645 steps\n",
      "1e-10 Broken\n",
      "Counted in: 5200 steps\n",
      "**********NORMAL**********\n",
      "20 LOOP COMPLETED size = 2000\n",
      "mean error = [0.00074879]\n",
      "**********PRECONDITIONED**********\n",
      "20 LOOP COMPLETED size = 2000\n",
      "mean error = [0.00325894]\n",
      "\n",
      "1e-10 Broken\n",
      "Counted in: 9528 steps\n",
      "1e-10 Broken\n",
      "Counted in: 6538 steps\n",
      "**********NORMAL**********\n",
      "21 LOOP COMPLETED size = 2100\n",
      "mean error = [0.01396835]\n",
      "**********PRECONDITIONED**********\n",
      "21 LOOP COMPLETED size = 2100\n",
      "mean error = [0.04499398]\n",
      "\n",
      "1e-10 Broken\n",
      "Counted in: 11440 steps\n",
      "1e-10 Broken\n",
      "Counted in: 7478 steps\n",
      "**********NORMAL**********\n",
      "22 LOOP COMPLETED size = 2200\n",
      "mean error = [0.00937217]\n",
      "**********PRECONDITIONED**********\n",
      "22 LOOP COMPLETED size = 2200\n",
      "mean error = [0.22849979]\n",
      "\n",
      "1e-10 Broken\n",
      "Counted in: 14626 steps\n",
      "1e-10 Broken\n",
      "Counted in: 9656 steps\n",
      "**********NORMAL**********\n",
      "23 LOOP COMPLETED size = 2300\n",
      "mean error = [0.08038733]\n",
      "**********PRECONDITIONED**********\n",
      "23 LOOP COMPLETED size = 2300\n",
      "mean error = [1.04281559]\n",
      "\n",
      "1e-10 Broken\n",
      "Counted in: 16873 steps\n",
      "1e-10 Broken\n",
      "Counted in: 10803 steps\n",
      "**********NORMAL**********\n",
      "24 LOOP COMPLETED size = 2400\n",
      "mean error = [1.52057264]\n",
      "**********PRECONDITIONED**********\n",
      "24 LOOP COMPLETED size = 2400\n",
      "mean error = [2.97912052]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "count = 25\n",
    "y = np.zeros(count-1)\n",
    "z = np.zeros(count-1)\n",
    "\n",
    "i = 1\n",
    "while i < count:\n",
    "    size = 100*i\n",
    "    A = gen_matrix(size)\n",
    "#    check_matrix(A)\n",
    "    x,b = gen_xb(size)\n",
    "    M = gen_precond(A, size)\n",
    "#    check_matrix(M)\n",
    "    CG = cg(A,b,x)[0]\n",
    "    My = my_cg(A,b,x)\n",
    "    My_p = my_cg_prec(A,b,x,M)\n",
    "    \n",
    "    print('**********NORMAL**********')\n",
    "    #normal\n",
    "    j = 0\n",
    "    summa = 0\n",
    "    while j < size:\n",
    "        #print(CG[j]-My[j])\n",
    "        summa += abs(CG[j]-My[j])\n",
    "        j = j + 1\n",
    "    print(i,\"LOOP COMPLETED size =\",size)\n",
    "    mer = summa/size\n",
    "    y[i-1] = mer\n",
    "    print('mean error =', mer)\n",
    "    \n",
    "    print('**********PRECONDITIONED**********')\n",
    "    #preconditioned\n",
    "    j = 0\n",
    "    summa = 0\n",
    "    while j < size:\n",
    "        #print(CG[j]-My_p[j])\n",
    "        summa += abs(CG[j]-My_p[j])\n",
    "        j = j + 1\n",
    "    print(i,\"LOOP COMPLETED size =\",size)\n",
    "    mer = summa/size\n",
    "    z[i-1] = mer\n",
    "    print('mean error =', mer)\n",
    "    \n",
    "#    M_inv @ A is better conditioned than A\n",
    "#    print('A condition number =', np.linalg.cond(A))\n",
    "#    M_inv = np.linalg.inv(M)\n",
    "#    print('M_inv @ A condition number =', np.linalg.cond(M_inv @ A))\n",
    "\n",
    "    print()\n",
    "    i = i + 1\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets plot mean error for each loop:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEMCAYAAAAIx/uNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAGM1JREFUeJzt3X+QZWV95/H3xxHILFriD8zqMITRGgfH6AbsRSMxRaUSGIwRJFaWQXddnRSLFdbkj5AM626JbtXCxiRb60KtjhFRY2aWNYSd0tkdDcYFVxKnEQKM1Og4atGDFcYf+JMo4Hf/uLeZpu3uuX37nHu6b79fVV1979P3nvu9p273p8/zPOc5qSokSWrSk7ouQJI0fgwXSVLjDBdJUuMMF0lS4wwXSVLjDBdJUuMMF0lS4wwXSVLjntx1AceS5JXA6+nVurmqXtFxSZKkY+jkyCXJ9UkeTHLvrPYtSQ4kOZhkO0BV3VZVlwEfAz7YRb2SpMXpqlvsBmDLzIYka4DrgPOBzcDWJJtnPOQSYOeoCpQkDa+TbrGqujXJabOazwIOVtUhgCS7gAuALyQ5FfhOVX13vm0muRS4FODEE0986emnn95G6ZI0tu64445vVNXJTWxrOY25rAPun3F/CnhZ//Y24AMLPbmqdgA7ACYmJmpycrKNGiVpbCX5WlPbWk7hkjnaCqCq3j7iWiRJS7CcpiJPAetn3D8FeKCjWiRJS7CcwmUfsDHJhiTHAxcDuzuuSZI0hK6mIu8Ebgc2JZlKsq2qHgUuB/YC9wE3VtX+LuqTJC1NV7PFts7TvgfYM+JyJEkNW07dYpKkMWG4SJIaZ7hIkhpnuEiSGme4SJIaZ7hIkhpnuEiSGme4SJIaZ7hIkhpnuEiSGme4SJIaZ7hIkhpnuEiSGme4SJIaZ7hIkhpnuEiSGme4SJIaZ7hIkhpnuEiSGme4SJIaZ7hIkhq37MMlyTlJbkvyniTndF2PJOnYOgmXJNcneTDJvbPatyQ5kORgku395gK+D/wMMDXqWiVJi9fVkcsNwJaZDUnWANcB5wObga1JNgO3VdX5wB8C7xhxnZKkIXQSLlV1K/CtWc1nAQer6lBV/RjYBVxQVT/p//zbwAnzbTPJpUkmk0weOXKklbolSYNZTmMu64D7Z9yfAtYluSjJe4EPA9fO9+Sq2lFVE1U1cfLJJ7dcqiRpIU/uuoAZMkdbVdVNwE2jLkaSNLzldOQyBayfcf8U4IGOapEkLcFyCpd9wMYkG5IcD1wM7O64JknSELqairwTuB3YlGQqybaqehS4HNgL3AfcWFX7u6hPkrQ0nYy5VNXWedr3AHtGXI4kqWHLqVtMkjQmDBdJUuMMF0lS4wwXSVLjDBdJUuMMF0lS4wwXSVLjDBdJUuMMF0lS4wwXSVLjDBdJUuMMF0lS4wwXSVLjDBdJUuMMF0lS4wwXSVLjDBdJUuMMF0lS4wwXSVLjDBdJUuMMF0lS45Z9uCR5YZL3JPlokrd0XY8k6dg6CZck1yd5MMm9s9q3JDmQ5GCS7QBVdV9VXQb8FjDRRb2SpMXp6sjlBmDLzIYka4DrgPOBzcDWJJv7P3sN8BngltGWKUkaRifhUlW3At+a1XwWcLCqDlXVj4FdwAX9x++uqlcAr59vm0kuTTKZZPLIkSNtlS5JGsCTuy5ghnXA/TPuTwEvS3IOcBFwArBnvidX1Q5gB8DExES1V6Yk6ViWU7hkjraqqk8Dnx5tKZKkpVhOs8WmgPUz7p8CPNBRLZKkJVhO4bIP2JhkQ5LjgYuB3R3XJEkaQldTkXcCtwObkkwl2VZVjwKXA3uB+4Abq2p/F/VJkpamkzGXqto6T/seFhi0lyStDMupW0ySNCYMF0lS4wwXSVLjDBdJUuMMF0lS4wwXSVLjDBdJUuMMF0lS4wwXSVLjDBdJUuMMF0lS4wwXSVLjDBdJUuMMF0lS4xYMlyRrkvz1qIqRJI2HBcOlqh4DfpjkaSOqR5I0Bga5WNg/Avck+STwg+nGqnpra1VJkla0QcLl4/0vSZIGcsxwqaoPJjkeeEG/6UBVPdJuWZKkleyY4ZLkHOCDwFeBAOuTvLGqbm23NEnSSjVIt9ifAOdW1QGAJC8AdgIvbbMwSdLKNch5LsdNBwtAVX0ROK69kp4oyfOSvD/JR0f1mpKkpRkkXCb7f9zP6X+9D7hjKS+a5PokDya5d1b7liQHkhxMsh2gqg5V1balvJ4kabQGCZe3APuBtwK/C3wBuGyJr3sDsGVmQ5I1wHXA+cBmYGuSzUt8HUlSBxYcc+n/wX9/Vb0B+NOmXrSqbk1y2qzms4CDVXWo/9q7gAvohZkkaQUZ5Az9k/tTkdu2Drh/xv0pYF2SZyZ5D3BGkivne3KSS5NMJpk8cuRI27VKkhYwyGyxrwL/L8lunniGfmNHMn2Zo62q6psM0A1XVTuAHQATExPVcG2SpEUYJFwe6H89CXhqi7VMAetn3D+l/7qSpBVmkDGXp1TVFSOoZR+wMckG4DBwMXDJCF5XktSwQcZczmz6RZPsBG4HNiWZSrKtqh4FLgf2AvcBN1bV/qZfW5LUvkG6xe7qj7f8T5445nLTsC9aVVvnad8D7Bl2u5Kk5WGQcHkG8E3gV2a0FTB0uEiSxtsgqyK/aRSFSJLGxzHP0E/ygiS3TC/VkuQlSf59+6VJklaqQZZ/eR9wJfAIQFXdTW8mlyRpGbn5zsOcfc2n2LD945x9zae4+c7DndUyyJjLP6mqzyVPOMfx0ZbqkSQN4eY7D3PlTffw8COPAXD4oYe58qZ7ALjwjHUjr2eQI5dvJHk+vUF8krwO+HqrVUmSFuVdew88HizTHn7kMd6198A8z2jXIEcuv0NvWZXTkxwGvgK8vtWqJEmL8sBDDy+qvW2DzBY7BPxqkhOBJ1XV99ovS5K0GM89aS2H5wiS5560toNqBusWA6CqfmCwSNLydMV5m1h73JontK09bg1XnLepk3oG6RaTJC1z04P279p7gAceepjnnrSWK87b1MlgPhgukjQ2LjxjXWdhMttA4ZLkFcBpMx9fVR9qqSZJ0gp3zHBJ8mHg+cBdwPQ8twIMF0nSnAY5cpkANleVV3eUJA1kkNli9wL/tO1CJEnjY5Ajl2cBX0jyOeBH041V9ZrWqpIkrWiDhMtVbRchSRovg5yh/39HUYgkaXwMcj2XlyfZl+T7SX6c5LEk3x1FcZKklWmQAf1rga3Al4C1wG/32yRJmtNAJ1FW1cEka6rqMeADST7bcl2SpBVskHD5YZLjgbuS/BG9a7mc2G5ZRyV5HvA24GlV9bpRva4kaXiDdIv9y/7jLgd+AKwHfnOQjSe5PsmDSe6d1b4lyYEkB5NsX2gbVXWoqrYN8nqSpOVhkNliX0uyFnhOVb1jkdu/gd74zONLxSRZA1wH/BowBexLshtYA1w96/lvrqoHF/makrTi3Xzn4WWzwvEwBllb7DeAPwaOBzYk+QXgnYOcRFlVtyY5bVbzWcDB/kXISLILuKCqrgZevbjyJWn83HznYa686Z7HL1t8+KGHufKmewBWTMAM0i12Fb1AeAigqu6it0LysNYB98+4P9Vvm1OSZyZ5D3BGkisXeNylSSaTTB45cmQJ5UlSt96198DjwTLt4Uce4117D3RU0eINMqD/aFV9J0lTrznXhuZdFLOqvglcdqyNVtUOYAfAxMSEi2xKWrHmu+79fO3L0UALVya5BFiTZGOS/wYsZSryFL1JAdNOAR5YwvYkaazMd937+dqXo0HC5d8CL6K3aOVO4LvA7y3hNfcBG5Ns6E9xvhjYvYTtSdJYueK8Taw9bs0T2tYet4YrztvUUUWLN8hssR/SO8/kbYvdeJKdwDnAs5JMAW+vqvcnuRzYS2+G2PVVtX+x25akcTU9aL+SZ4vlWNcASzIB/Dt++jLHL2m1siWYmJioycnJrsuQpBUlyR1VNdHEtgYZ0P8IcAVwD/CTJl5UkjTeBgmXI1XlmIgkaWCDhMvbk/wZcAtPvBLlTa1VJUla0QYJlzcBpwPHcbRbrADDRZI0p0HC5Z9V1Ytbr0SSlrmVvt7XKA1ynsvfJtnceiWStIxNr/d1+KGHKY6u93XznYe7Lm1ZGiRcfonetVwOJLk7yT1J7m67MElaTsZhva9RGqRbbEvrVUjSMjfsel+rtSttoOu5jKIQSVrOnnvSWg7PESQLrfc1DkvnD2uQbjFJWvWGWe9rNXelDdItJkmr3jDrfY3D0vnDMlwkaUAXnrFuUd1Zw3SljQu7xSSpJeOwdP6wPHKRpJaMw9L5wzJcJKlFi+1KGxd2i0mSGme4SJIaZ7hIkhpnuEiSGme4SJIaZ7hIkhq37MMlyQuTvCfJR5O8pet6JEnH1mq4JLk+yYNJ7p3VvqV/fZiDSbYvtI2quq+qLgN+C5hos15JUjPaPnK5gVnXg0myBrgOOB/YDGxNsjnJi5N8bNbXs/vPeQ3wGeCWluuVJDWg1TP0q+rWJKfNaj4LOFhVhwCS7AIuqKqrgVfPs53dwO4kHwf+or2KJa0Wq/UiXqPSxfIv64D7Z9yfAl4234OTnANcBJwA7FngcZcClwKceuqpTdQpaUyt5ot4jUoX4ZI52mq+B1fVp4FPH2ujVbUD2AEwMTEx7/YkjZdhjkAWuoiX4dKMLsJlClg/4/4pwAMd1CFphRv2CGQ1X8RrVLqYirwP2JhkQ5LjgYuB3R3UIWmFG/YywvNdrGs1XMRrVNqeirwTuB3YlGQqybaqehS4HNgL3AfcWFX726xD0spw852HOfuaT7Fh+8c5+5pPcfOdhxd8/LBHIKv5Il6j0vZssa3ztO9hgcF5ScvHMGMawz5nsV1cw15GeDVfxGtUvFiYpHkN8wd/2HGQYQbZrzhv0xNeCwY/AlmtF/EalWW//Iuk5iy222mYMY1hx0GG6eK68Ix1XH3Ri1l30loCrDtpLVdf9GJDYxnwyEVaJYY5ohjmD/6w4yBL6eIyTJYfj1ykFWixRyAw3BHFMLOqhp2J5SD7eDFcpI4tNiimj0AOP/QwxdEjkDZmVg3zB3/YkLCLa7zYLSZ1aJiuqmHPLh+m22mYWVVLmYllF9f4MFykDg0TFEs5t2OYmVXD/ME3JGS4SA1a7PkdwwSF53ZoJTBcpDks55MAPbdDK4ED+tIsww6YDzMba5jBbwe+tRJ45CLNMuyA+bAnAU6/5mKOkjwC0XJnuGjsjWIcBDwJUJrJbjGtGMOcODhMF5cnAUpL55GLlmwUq+auhMUQnY0lHWW46HGjmiE1yhMHRzkOMv1cw0QyXMbWcj4yGOWJg46DSN1wzGUMDTPOMMpl0oc9cXAx7dMcB5G6YbiMoWGCYilHBotpH/Y5LoYorSx2i42hUS4pMszg9zDPcRxEWlkMlzE0yiVFRrlqriEhrRypqq5raNzExERNTk52XUZnZg/OQy8ojtUdNMxsMUnjI8kdVTXRxLaW/ZFLknOA/wjsB3ZV1ac7LWgJRnE+CHhkIKl7rYZLkuuBVwMPVtXPz2jfAvxXYA3wZ1V1zQKbKeD7wM8AUy2W26pRnQ8yzaCQ1KW2Z4vdAGyZ2ZBkDXAdcD6wGdiaZHOSFyf52KyvZwO3VdX5wB8C72i53tYMM4Nr2OnBktS1Vo9cqurWJKfNaj4LOFhVhwCS7AIuqKqr6R3lzOfbwAlt1DmMUSyGOOz0YEnqWhfnuawD7p9xf6rfNqckFyV5L/Bh4NoFHndpkskkk0eOHGms2LmMajHEYU8clKSudREumaNt3ilrVXVTVf2bqvoXCw3mV9WOqpqoqomTTz65iTrnNaqLQnl2uaSVqovZYlPA+hn3TwEe6KCOoY1qMURX2ZW0UnURLvuAjUk2AIeBi4FLOqhjaKNcDNFZX5JWola7xZLsBG4HNiWZSrKtqh4FLgf2AvcBN1bV/jbraJrdVZK0sLZni22dp30PsKfN1x7UKE9SlKTVYtmfod8mT1KUpHas6iX3PUlRktqxqsPFkxQlqR2rOlw8SVGS2rGqw8VZX5LUjlU9oO+sL0lqx6oOF3DWlyS1YVV3i0mS2mG4SJIaZ7hIkhpnuEiSGme4SJIaZ7hIkhpnuEiSGme4SJIaZ7hIkhpnuEiSGme4SJIaZ7hIkhpnuEiSGme4SJIat+yX3E/ySuD19GrdXFWv6LgkSdIxtHrkkuT6JA8muXdW+5YkB5IcTLJ9oW1U1W1VdRnwMeCDbdYrSWpG20cuNwDXAh+abkiyBrgO+DVgCtiXZDewBrh61vPfXFUP9m9fAvx2y/VKkhrQarhU1a1JTpvVfBZwsKoOASTZBVxQVVcDr55rO0lOBb5TVd9tsVxJUkO6GHNZB9w/4/4U8LJjPGcb8IGFHpDkUuDS/t0fze6KW8WeBXyj6yKWCffFUe6Lo9wXR21qakNdhEvmaKuFnlBVbz/WRqtqB7ADIMlkVU0MV954cV8c5b44yn1xlPviqCSTTW2ri6nIU8D6GfdPAR7ooA5JUku6CJd9wMYkG5IcD1wM7O6gDklSS9qeirwTuB3YlGQqybaqehS4HNgL3AfcWFX7G37pHQ1vbyVzXxzlvjjKfXGU++KoxvZFqhYc7pAkadFc/kWS1DjDRZLUuLEKl8UsKzMuknw1yT1J7pqeRpjkGUk+meRL/e9P77cnybv7++fuJGd2W/3SzLW80DDvPckb+4//UpI3dvFelmqefXFVksP9z8ZdSV4142dX9vfFgSTnzWhf8b9DSdYn+Zsk9yXZn+R3++2r7rOxwL5o/7NRVWPxRW/5mC8DzwOOB/6e3kKXndfW8vv+KvCsWW1/BGzv394O/Of+7VcB/5veuUYvB/6u6/qX+N5/GTgTuHfY9w48AzjU//70/u2nd/3eGtoXVwG/P8djN/d/P04ANvR/b9aMy+8Q8BzgzP7tpwJf7L/nVffZWGBftP7ZGKcjl8eXlamqHwO7gAs6rqkrF3B0kc8PAhfOaP9Q9fwtcFKS53RRYBOq6lbgW7OaF/vezwM+WVXfqqpvA58EtrRffbPm2RfzuQDYVVU/qqqvAAfp/f6Mxe9QVX29qj7fv/09erNS17EKPxsL7Iv5NPbZGKdwmWtZmYV24rgo4BNJ7ugvgQPws1X1deh9uIBn99tXwz5a7Hsf931yeb+r5/rpbiBW0b7or214BvB3rPLPxqx9AS1/NsYpXBa9rMyYOLuqzgTOB34nyS8v8NjVuo9g/vc+zvvkvwPPB34B+DrwJ/32VbEvkjwF+Evg92rhRW/Hfn/MsS9a/2yMU7isymVlquqB/vcHgb+id/j6D9PdXf3v05ctWA37aLHvfWz3SVX9Q1U9VlU/Ad5H77MBq2BfJDmO3h/Tj1TVTf3mVfnZmGtfjOKzMU7hsuqWlUlyYpKnTt8GzgXupfe+p2e2vBH4X/3bu4F/1Z8d83J6lzH4+ojLbtti3/te4NwkT+93DZzbb1vxZo2nvZbeZwN6++LiJCck2QBsBD7HmPwOJQnwfuC+qvrTGT9adZ+N+fbFSD4bXc9maHhmxKvozYb4MvC2rusZwft9Hr1ZG38P7J9+z8AzgVuAL/W/P6PfHnoXavsycA8w0fV7WOL730nvkP4Rev9ZbRvmvQNvpjdweRB4U9fvq8F98eH+e727/4fgOTMe/7b+vjgAnD+jfcX/DgG/RK/L5m7grv7Xq1bjZ2OBfdH6Z8PlXyRJjRunbjFJ0jJhuEiSGme4SJIaZ7hIkhpnuEiSGme4aNVK8tokleT0AR9/VZLf799+Z5JfbbfCOWt4zUpdrViri1ORtWoluZHeqrG3VNVVAzz+KuD7VfXHLZcmrXgeuWhV6q+1dDa9kw0vXuBxb+tfw+KvgU0z2m9I8rr+7a8m+U9Jbk8ymeTMJHuTfDnJZTOec0WSff3FAt/Rbzutf62N9/Wvt/GJJGv7P3trki/0H7+r3/avk1zbv/1zSW7p//yWJKfOqO3dST6b5NB0ndIoGS5arS4E/k9VfRH4Vua4cFqSl9ILnjOAi4B/vsD27q+qXwRuA24AXkfv2iDv7G/rXHpLaZxFb7HAl85YZHQjcF1VvQh4CPjNfvt24IyqegnweEjNcC29peJfAnwEePeMnz2H3tnZrwauWaBuqRWGi1arrfSuSUH/+9Y5HvNK4K+q6ofVW0l2obWUpn92D72LTX2vqo4A/5jkJHrrUp0L3Al8HjidXqgAfKWq7urfvgM4rX/7buAjSd4APDrHa/4i8Bf92x+mFybTbq6qn1TVF4CfXaBuqRVP7roAadSSPBP4FeDnkxS9q+xVkj+onx6EHHRQ8kf97z+ZcXv6/pPprV91dVW9d1Ytp816/GPA2v7tX6d3hcnXAP8hyYuOUcPMWmduc67l0qVWeeSi1eh19LqTfq6qTquq9cBXeOJ//gC3Aq9Nsra/+vRvLOE19wJv7o/1kGRdkmfP9+AkTwLWV9XfAH8AnAQ8ZdbDPsvR8aLXA59ZQn1Sozxy0Wq0lZ8eh/hL4BJ6YyYAVNXnk/wPeivJfm3mzxarqj6R5IXA7b1V0Pk+8AZ6RypzWQP8eZKn0Tvy+C9V9VD/udPeClyf5ArgCPCmYeuTmuZUZElS4+wWkyQ1znCRJDXOcJEkNc5wkSQ1znCRJDXOcJEkNc5wkSQ17v8DvmOXJVGClxQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#for normal\n",
    "plt.xlim([0, 2500])\n",
    "plt.ylim([1e-7, 1e+7])\n",
    "\n",
    "x_axis = np.arange(1,i)*100\n",
    " \n",
    "# Plot\n",
    "plt.scatter(x_axis, y)\n",
    "plt.xlabel('A dimension')\n",
    "plt.ylabel('mean error')\n",
    "plt.yscale('log')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEMCAYAAAAIx/uNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAGL1JREFUeJzt3X+wXGd93/H3B2E7qmEifpgUZDkWjJERMY3NjSEQMp5MYsuUYOMwqWXTUlDGdSYupDNxYpd2MHSmpiE/pow9ARGMgVKpLnFcDVEriBMqU5wgGTtYwiMQMoyvzEQCY34asM23f+xe6/py79Xevefs7t19v2buaPe5Z8/57pm9+9E5z3Oek6pCkqQmPWXYBUiSxo/hIklqnOEiSWqc4SJJapzhIklqnOEiSWqc4SJJapzhIklq3FOHXcDxJHkVcDmdWjdW1SuGXJIk6TiGcuSS5KYkR5Lsm9O+KcmBJAeTXANQVXdU1ZXAx4EPDaNeSdLSDOu02M3AptkNSVYBNwIXAhuBzUk2zlrkMmDboAqUJPVvKKfFqmp3ktPnNJ8LHKyqQwBJtgMXAV9Ichrwrar69kLrTHIFcAXAySef/NIzzzyzjdIlaWzdddddX6+qU5pY1yj1uawFHpj1fBp4WffxFuCDi724qrYCWwGmpqZq7969bdQoSWMryVebWtcohUvmaSuAqnr7gGuRJC3DKA1FngbWzXp+KvDgkGqRJC3DKIXLHuCMJOuTnAhcCuwYck2SpD4MayjyNuBOYEOS6SRbquox4CpgF3AfcEtV7R9GfZKk5RnWaLHNC7TvBHYOuBxJUsNG6bSYJGlMGC6SpMYZLpKkxhkukqTGGS6SpMYZLpKkxhkukqTGGS6SpMYZLpKkxhkukqTGGS6SpMYZLpKkxhkukqTGGS6SpMYZLpKkxhkukqTGGS6SpMYZLpKkxhkukqTGGS6SpMYZLpKkxo18uCQ5L8kdSd6b5Lxh1yNJOr6hhEuSm5IcSbJvTvumJAeSHExyTbe5gO8CPwVMD7pWSdLSDevI5WZg0+yGJKuAG4ELgY3A5iQbgTuq6kLgD4B3DLhOSVIfhhIuVbUbeGhO87nAwao6VFU/ArYDF1XVj7u//yZw0kLrTHJFkr1J9h49erSVuiVJvRmlPpe1wAOznk8Da5NckuR9wEeAGxZ6cVVtraqpqpo65ZRTWi5VkrSYpw67gFkyT1tV1a3ArYMuRpLUv1E6cpkG1s16firw4JBqkSQtwyiFyx7gjCTrk5wIXArsGHJNkqQ+DGso8jbgTmBDkukkW6rqMeAqYBdwH3BLVe0fRn2SpOUZSp9LVW1eoH0nsHPA5UiSGjZKp8UkSWPCcJEkNc5wkSQ1znCRJDXOcJEkNc5wkSQ1znCRJDXOcJEkNc5wkSQ1znCRJDXOcJEkNc5wkSQ1znCRJDXOcJEkNc5wkSQ1znCRJDXOcJEkNc5wkSQ1znCRJDXOcJEkNc5wkSQ1buTDJcmLkrw3yceS/Paw65EkHd9QwiXJTUmOJNk3p31TkgNJDia5BqCq7quqK4HfBKaGUa8kaWmGdeRyM7BpdkOSVcCNwIXARmBzko3d370W+DRw+2DLlCT1YyjhUlW7gYfmNJ8LHKyqQ1X1I2A7cFF3+R1V9Qrg8oXWmeSKJHuT7D169GhbpUuSevDUYRcwy1rggVnPp4GXJTkPuAQ4Cdi50IuraiuwFWBqaqraK1OSdDyjFC6Zp62q6lPApwZbiiRpOUZptNg0sG7W81OBB4dUiyRpGUYpXPYAZyRZn+RE4FJgx5BrkiT1YVhDkbcBdwIbkkwn2VJVjwFXAbuA+4Bbqmr/MOqTJC3PUPpcqmrzAu07WaTTXpK0MozSaTFJ0pgwXCRJjTNcJEmNM1wkSY0zXCRJjTNcJEmNM1wkSY0zXCRJjTNcJEmNM1wkSY0zXCRJjTNcJEmNM1wkSY0zXCRJjVs0XJKsSvLXgypGkjQeFg2Xqnoc+H6Snx5QPZKkMdDLzcJ+ANyb5JPA92Yaq+otrVUlSVrRegmXv+r+SJLUk+OGS1V9KMmJwAu7TQeq6tF2y5IkrWTHDZck5wEfAr4CBFiX5I1Vtbvd0iRJK1Uvp8X+GDi/qg4AJHkhsA14aZuFSZJWrl6uczlhJlgAquqLwAntlfRkSZ6f5ANJPjaobUqSlqeXcNnb/XI/r/vzfuCu5Ww0yU1JjiTZN6d9U5IDSQ4muQagqg5V1ZblbE+SNFi9hMtvA/uBtwBvBb4AXLnM7d4MbJrdkGQVcCNwIbAR2Jxk4zK3I0kagkX7XLpf+B+oqjcAf9LURqtqd5LT5zSfCxysqkPdbW8HLqITZpKkFaSXK/RP6Q5Fbtta4IFZz6eBtUmeleS9wNlJrl3oxUmuSLI3yd6jR4+2XaskaRG9jBb7CvD/kuzgyVfoN3Yk05V52qqqvkEPp+GqaiuwFWBqaqoark2StAS9hMuD3Z+nAE9vsZZpYN2s56d2tytJWmF66XN5WlVdPYBa9gBnJFkPHAYuBS4bwHYlSQ3rpc/lnKY3mmQbcCewIcl0ki1V9RhwFbALuA+4par2N71tSVL7ejktdk+3v+V/8uQ+l1v73WhVbV6gfSews9/1StIku+3uw7x71wEefPgRnrdmNVdfsIGLz147lFp6CZdnAt8AfmVWWwF9h4skqVm33X2Ya2+9l0cefRyAww8/wrW33gswlIDpZVbkNw2iEElS/96968ATwTLjkUcf5927DgwlXI57hX6SFya5fWaqliQvSfIf2i9NktSrBx9+ZEntbetl+pf3A9cCjwJU1efpjOSSJI2I561ZvaT2tvUSLv+kqj47p+2xNoqRJPXn6gs2sPqEVU9qW33CKq6+YMNQ6umlQ//rSV5ApxOfJK8HvtZqVZKkJZnpV1lJo8V+h860KmcmOQzcD1zealWSNOH6GVZ88dlrhxYmc/UyWuwQ8KtJTgaeUlXfab8sSZpcozasuB+99LkAUFXfM1gkqX2LDSteKXoOF0nSYIzasOJ+GC6SNGJGbVhxP3rp0CfJK4DTZy9fVR9uqSZJmmhXX7DhSX0uMNxhxf04brgk+QjwAuAeYOadFmC4SFILRm1YcT96OXKZAjZWlXd3lKQBGaVhxf3opc9lH/BP2y5EkjQ+ejlyeTbwhSSfBX4401hVr22tKknSitZLuFzXdhGSNK5G6QZeg9TLFfr/dxCFSNK4GYcr7fvVy/1cXp5kT5LvJvlRkseTfHsQxUnSSjYOV9r3q5cO/RuAzcCXgNXAb3XbJEmLGIcr7fvV0xX6VXUQWFVVj1fVB4HzWq1KksbAOFxp369ewuX7SU4E7knyh0n+HXByy3U9Icnzk3wgyccGtU1JasKo3cBrkHoJl3/ZXe4q4HvAOuA3ell5kpuSHEmyb077piQHkhxMcs1i66iqQ1W1pZftSdIoufjstVx/yVmsXbOaAGvXrOb6S84a+8586G202FeTrAaeW1XvWOL6b6bTP/PEVDFJVgE3Ar8GTAN7kuwAVgHXz3n9m6vqyBK3KUkjY6Vfad+vXuYW+3Xgj4ATgfVJfh54Zy8XUVbV7iSnz2k+FzjYvQkZSbYDF1XV9cBrlla+JGkU9XJa7Do6gfAwQFXdQ2eG5H6tBR6Y9Xy62zavJM9K8l7g7CTXLrLcFUn2Jtl79OjRZZQnSVquXq7Qf6yqvpWkqW3Ot6IFJ8Wsqm8AVx5vpVW1FdgKMDU15SSbkjREvYTLviSXAauSnAG8BfjMMrY5TWdQwIxTgQeXsT5JGohJncqlH72cFvu3wIvpTFq5Dfg28LvL2OYe4Iwk67tDnC8FdixjfZLUupmpXA4//AjFsalcbrv78LBLG0nHDZeq+n5Vva2qfqGqprqPf9DLypNsA+4ENiSZTrKlqh6jM6x5F3AfcEtV7V/Om5Cktk3yVC796GW02BTw7/nJ2xy/5HivrarNC7TvBHb2XKUkDdkkT+XSj176XD4KXA3cC/y43XIkaTQ9b81qDs8TJJMwlUs/eulzOVpVO6rq/qr66sxP65VJ0giZ5Klc+tHLkcvbk/w5cDtPvhPlra1VJUkjZmZUmKPFetNLuLwJOBM4gWOnxQowXCRNlEmdyqUfvYTLP6uqs1qvRJI0Nnrpc/m7JBtbr0SSNDZ6OXL5JeCNSe6n0+cSoHoZiixJmky9hMum1quQJI2Vnu7nMohCJEnjo5c+F0mSlsRwkSQ1znCRJDXOcJEkNc5wkSQ1rpehyJI0dryrZLsMF0kTZ+aukjM3/5q5qyRgwDTE02KSJo53lWyf4SJp4nhXyfYZLpImzkJ3j/Suks0xXCStaLfdfZhXvutvWH/NX/HKd/0Nt919+Liv8a6S7bNDX9KK1W/HvHeVbN/Ih0uSFwFvBZ4N3F5VfzbkkiSNiMU65o8XFN5Vsl2tnhZLclOSI0n2zWnflORAkoNJrllsHVV1X1VdCfwmMNVmvZJWFjvmR1fbfS43M+d+MElWATcCFwIbgc1JNiY5K8nH5/w8p/ua1wKfBm5vuV5JK4gd86Or1XCpqt3AQ3OazwUOVtWhqvoRsB24qKrurarXzPk50l3Pjqp6BXB5m/VKWlnsmB9dw+hzWQs8MOv5NPCyhRZOch5wCXASsHOR5a4ArgA47bTTmqhT0oizY350DSNcMk9bLbRwVX0K+NTxVlpVW4GtAFNTUwuuT9J4sWN+NA3jOpdpYN2s56cCDw6hDklSS4YRLnuAM5KsT3IicCmwYwh1SJJa0uppsSTbgPOAZyeZBt5eVR9IchWwC1gF3FRV+9usQ9LK4DT446PVcKmqzQu072SRznlJo6OfL/x+X+M0+OPDucUkLWjmC//ww49QHPvCX2z+rn5eA06DP24MF2mCLHWSx36+8PsNCa+2Hy+GizQh+jmi6OcLv9+Q8Gr78WK4SCtQP9PM93NE0c8Xfr8h4dX248VwkYZsqUHRb59GP0cU/Xzh9xsSF5+9lusvOYu1a1YTYO2a1Vx/yVl25q9QIz/lvjTO+hkh1e80889bs5rD8wTJYkcU/UyvspwpWbzafnwYLtIQ9RMU/fZpXH3BhicFGfR+RLHUL3xDQoaL1KClXt/RT1D0cwQCTvKowTJcpHkM6iLAfoKi3yOQmToMEw2CHfrSHIO8CLCfzm87vrUSeOQizdFvh3k/p7j6PVXlEYhGneGisTeIfhBYXl+IQaFx42kxrRj9XDjYzykuLwKUls8jFy3bIGbN7XfG3H5OcS1nyO7MNh2NpUlnuOgJgxohNcgLBwfZDzLzWsNEMlzG1igfGQzywkH7QaThsM9lDPXTzzDIadL7vXBwKe0z7AeRhsNwGUP9BMUgp0nv5zVOhiitLJ4WG0ODnFKkn87vfl5jP4i0shguY2iQU4oMctZcQ0JaOVJVw66hcVNTU7V3795hlzE0czvnoRMUxzsd1M9oMUnjI8ldVTXVxLpG/sglyXnAfwL2A9ur6lNDLWgZBnE9CHhkIGn4Wg2XJDcBrwGOVNXPzWrfBPxXYBXw51X1rkVWU8B3gZ8Cplsst1WDuh5khkEhaZjaHi12M7BpdkOSVcCNwIXARmBzko1Jzkry8Tk/zwHuqKoLgT8A3tFyva3pZwRXv8ODJWnYWj1yqardSU6f03wucLCqDgEk2Q5cVFXX0znKWcg3gZPaqLMfg5gMsd/hwZI0bMO4zmUt8MCs59PdtnkluSTJ+4CPADcsstwVSfYm2Xv06NHGip3PoCZD7PfCQUkatmGES+ZpW3DIWlXdWlX/pqr+xWKd+VW1taqmqmrqlFNOaaLOBQ3qplBeXS5ppRrGaLFpYN2s56cCDw6hjr4NajJEZ9mVtFINI1z2AGckWQ8cBi4FLhtCHX0b5GSIjvqStBK1elosyTbgTmBDkukkW6rqMeAqYBdwH3BLVe1vs46mebpKkhbX9mixzQu07wR2trntXg3yIkVJmhQjf4V+m7xIUZLaMdFT7nuRoiS1Y6LDxYsUJakdEx0uXqQoSe2Y6HBx1JcktWOiO/Qd9SVJ7ZjocAFHfUlSGyb6tJgkqR2GiySpcYaLJKlxhoskqXGGiySpcYaLJKlxhoskqXGGiySpcYaLJKlxhoskqXGGiySpcYaLJKlxhoskqXGGiySpcSM/5X6SVwGX06l1Y1W9YsglSZKOo9UjlyQ3JTmSZN+c9k1JDiQ5mOSaxdZRVXdU1ZXAx4EPtVmvJKkZbR+53AzcAHx4piHJKuBG4NeAaWBPkh3AKuD6Oa9/c1Ud6T6+DPitluuVJDWg1XCpqt1JTp/TfC5wsKoOASTZDlxUVdcDr5lvPUlOA75VVd9usVxJUkOG0eeyFnhg1vNp4GXHec0W4IOLLZDkCuCK7tMfzj0VN8GeDXx92EWMCPfFMe6LY9wXx2xoakXDCJfM01aLvaCq3n68lVbVVmArQJK9VTXVX3njxX1xjPviGPfFMe6LY5LsbWpdwxiKPA2sm/X8VODBIdQhSWrJMMJlD3BGkvVJTgQuBXYMoQ5JUkvaHoq8DbgT2JBkOsmWqnoMuArYBdwH3FJV+xve9NaG17eSuS+OcV8c4744xn1xTGP7IlWLdndIkrRkTv8iSWqc4SJJatxYhctSppUZF0m+kuTeJPfMDCNM8swkn0zype6/z+i2J8l7uvvn80nOGW71yzPf9EL9vPckb+wu/6UkbxzGe1muBfbFdUkOdz8b9yR59azfXdvdFweSXDCrfcX/DSVZl+Rvk9yXZH+St3bbJ+6zsci+aP+zUVVj8UNn+pgvA88HTgT+gc5El0OvreX3/RXg2XPa/hC4pvv4GuC/dB+/GvjfdK41ejnw98Ouf5nv/ZeBc4B9/b534JnAoe6/z+g+fsaw31tD++I64PfmWXZj9+/jJGB99+9m1bj8DQHPBc7pPn468MXue564z8Yi+6L1z8Y4Hbk8Ma1MVf0I2A5cNOSahuUijk3y+SHg4lntH66OvwPWJHnuMApsQlXtBh6a07zU934B8Mmqeqiqvgl8EtjUfvXNWmBfLOQiYHtV/bCq7gcO0vn7GYu/oar6WlV9rvv4O3RGpa5lAj8bi+yLhTT22RincJlvWpnFduK4KOATSe7qToED8DNV9TXofLiA53TbJ2EfLfW9j/s+uap7quemmdNATNC+6M5teDbw90z4Z2POvoCWPxvjFC5LnlZmTLyyqs4BLgR+J8kvL7LspO4jWPi9j/M++TPgBcDPA18D/rjbPhH7IsnTgL8AfrcWn/R27PfHPPui9c/GOIXLRE4rU1UPdv89AvwlncPXf5w53dX9d+a2BZOwj5b63sd2n1TVP1bV41X1Y+D9dD4bMAH7IskJdL5MP1pVt3abJ/KzMd++GMRnY5zCZeKmlUlycpKnzzwGzgf20XnfMyNb3gj8r+7jHcC/6o6OeTmd2xh8bcBlt22p730XcH6SZ3RPDZzfbVvx5vSnvY7OZwM6++LSJCclWQ+cAXyWMfkbShLgA8B9VfUns341cZ+NhfbFQD4bwx7N0PDIiFfTGQ3xZeBtw65nAO/3+XRGbfwDsH/mPQPPAm4HvtT995nd9tC5UduXgXuBqWG/h2W+/210DukfpfM/qy39vHfgzXQ6Lg8Cbxr2+2pwX3yk+14/3/0ieO6s5d/W3RcHgAtnta/4vyHgl+icsvk8cE/359WT+NlYZF+0/tlw+hdJUuPG6bSYJGlEGC6SpMYZLpKkxhkukqTGGS6SpMYZLppYSV6XpJKc2ePy1yX5ve7jdyb51XYrnLeG167U2Yo1WRyKrImV5BY6s8beXlXX9bD8dcB3q+qPWi5NWvE8ctFE6s619Eo6Fxteushyb+vew+KvgQ2z2m9O8vru468k+c9J7kyyN8k5SXYl+XKSK2e95uoke7qTBb6j23Z6914b7+/eb+MTSVZ3f/eWJF/oLr+92/avk9zQffyzSW7v/v72JKfNqu09ST6T5NBMndIgGS6aVBcD/6eqvgg8lHlunJbkpXSC52zgEuAXFlnfA1X1i8AdwM3A6+ncG+Sd3XWdT2cqjXPpTBb40lmTjJ4B3FhVLwYeBn6j234NcHZVvQR4IqRmuYHOVPEvAT4KvGfW755L5+rs1wDvWqRuqRWGiybVZjr3pKD77+Z5lnkV8JdV9f3qzCS72FxKM7+7l87Npr5TVUeBHyRZQ2deqvOBu4HPAWfSCRWA+6vqnu7ju4DTu48/D3w0yRuAx+bZ5i8C/737+CN0wmTGbVX146r6AvAzi9QtteKpwy5AGrQkzwJ+Bfi5JEXnLnuV5PfrJzshe+2U/GH33x/Pejzz/Kl05q+6vqreN6eW0+cs/ziwuvv4n9O5w+Rrgf+Y5MXHqWF2rbPXOd906VKrPHLRJHo9ndNJP1tVp1fVOuB+nvw/f4DdwOuSrO7OPv3ry9jmLuDN3b4ekqxN8pyFFk7yFGBdVf0t8PvAGuBpcxb7DMf6iy4HPr2M+qRGeeSiSbSZn+yH+AvgMjp9JgBU1eeS/A86M8l+dfbvlqqqPpHkRcCdnVnQ+S7wBjpHKvNZBfy3JD9N58jjT6vq4e5rZ7wFuCnJ1cBR4E391ic1zaHIkqTGeVpMktQ4w0WS1DjDRZLUOMNFktQ4w0WS1DjDRZLUOMNFktS4/w8Zhs1q0h5kDwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#for preconditioned\n",
    "plt.xlim([0, 2500])\n",
    "plt.ylim([1e-7, 1e+7])\n",
    "\n",
    "x_axis = np.arange(1,i)*100\n",
    " \n",
    "# Plot\n",
    "plt.scatter(x_axis, z)\n",
    "plt.xlabel('A dimension')\n",
    "plt.ylabel('mean error')\n",
    "plt.yscale('log')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we have unlimited steps preconditioned algorith is cheaper, but gives worse answer\n",
    "With steps limited to n the preconditioned algorithm is cheaper and more precise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Literature:\n",
    "1. Амосов А.А., Дубинский Ю.А., Копченова Н.В. Вычислительные методы для инженеров. http://pnu.edu.ru/media/filer_public/2013/04/10/6-1_amosov-dubinsky-kopchenova_1994.pdf p284.\n",
    "2. Wiki article about Conjugate Gradient method \n",
    "https://en.wikipedia.org/wiki/Conjugate_gradient_method\n",
    "3. Preconditioning: An Introduction to the Conjugate Gradient Method Without the Agonizing Pain http://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf p39\n",
    "3. How to randomly generate sparse matrices in Python: https://docs.scipy.org/doc/scipy/reference/sparse.html <br>Or BETTER: https://stackoverflow.com/questions/26895011/python-how-to-use-python-to-generate-a-random-sparse-symmetric-matrix <br> Or even BETTER (used in current version): https://scikit-learn.org/stable/modules/generated/sklearn.datasets.make_sparse_spd_matrix.html\n",
    "\n",
    "Honorable mentions: \n",
    "1. The Concept of Conjugate Gradient Descent in Python. http://ikuz.eu/2015/04/15/the-concept-of-conjugate-gradient-descent-in-python/\n",
    "2. Iterative methods lab. Jacobi preconditioned conjugate gradients method. https://www.kth.se/social/files/5885d039f2765429974418ce/Lab1.pdf\n"
   ]
  }
 ],
 "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}