数値計算の誤差

numerical_errors.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import (\n",
    "    print_function,\n",
    "    division,\n",
    "    absolute_import,\n",
    "    unicode_literals)\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 有限で数表現することによる誤差\n",
    "\n",
    "単に格納するだけで誤差を生じる．内部が10進数ではなく2進数なので，10進数でキリがよくても表現できない"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.2]\n",
      "[ 0.200000000000000011102230246252]\n"
     ]
    }
   ],
   "source": [
    "a = np.empty(1, dtype=np.float64)\n",
    "a[:] = 0.2\n",
    "np.set_printoptions(precision=8) # default\n",
    "print(a)\n",
    "np.set_printoptions(precision=30)\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 桁落ちの例\n",
    "\n",
    "本来は負にならない（二乗平均 - 平均の二乗）が負になる"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-1.16415321827e-10\n"
     ]
    }
   ],
   "source": [
    "n = 1000\n",
    "x = 1000 + np.random.randn(n) * 10**-5\n",
    "print(np.mean(x**2) - np.mean(x)**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "大きな絶対値を引いて，誤差に対する有効桁数を相対的に大きくしていると問題を生じにくい"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9.50575571567e-11\n"
     ]
    }
   ],
   "source": [
    "y = x - 1000\n",
    "print(np.mean(y**2) - np.mean(y)**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## シグモイド関数の桁あふれ\n",
    "\n",
    "機械学習では非常によく使われるシグモイド関数は容易に桁あふれを起こす\n",
    "\\\\[ f(x) = \\frac{1}{1 + \\exp(-x)} \\\\]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sig(x):\n",
    "    return 1 / (1 + np.exp(-x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "数学上は 0 や 1 の値をとることはないシグモイド関数だが，指数関数を含むため容易に桁あふれを起こす"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 0.5\n",
      "5 0.993307149076\n",
      "10 0.999954602131\n",
      "15 0.999999694098\n",
      "20 0.999999997939\n",
      "25 0.999999999986\n",
      "30 1.0\n",
      "35 1.0\n",
      "40 1.0\n",
      "45 1.0\n",
      "True\n"
     ]
    }
   ],
   "source": [
    "for x in range(0, 50, 5):\n",
    "    print(x, sig(x))\n",
    "print(sig(50) == 1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 大きな数と小さな数の和\n",
    "\n",
    "小さな数の総和：最後の方は大きな数と小さな数の和は，大きな数にとって誤差となる"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.91624494050302e-12\n"
     ]
    }
   ],
   "source": [
    "k = 5\n",
    "a = 0\n",
    "for i in range(10**k):\n",
    "    a += 10**-k\n",
    "print(1 - a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "分けて計算すれば影響は小さい"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-6.661338147750939e-16\n"
     ]
    }
   ],
   "source": [
    "a = 0.\n",
    "for i1 in range(10):\n",
    "    b = 0.\n",
    "    for i2 in range(10):\n",
    "        c = 0.\n",
    "        for i2 in range(10):\n",
    "            d = 0.\n",
    "            for i3 in range(10):\n",
    "                e = 0.\n",
    "                for i3 in range(10):\n",
    "                    e += 10**-k\n",
    "                d +=e\n",
    "            c += d\n",
    "        b += c\n",
    "    a += b\n",
    "print(1 - a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 桁落ちの例\n",
    "\n",
    "※ p.8，伊理正夫，藤野和建「数値計算の常識」共立出版\n",
    "\n",
    "\\\\[1 - \\frac{1}{\\sqrt{1 + x}} = \\frac{\\sqrt{1 + x} - 1}{\\sqrt{1 + x}}\\\\]\n",
    "\n",
    "$x$が小さいと $\\sqrt{1+x}$ と $1$ の差が小さくなり桁落ちを生じる"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.57079682594e-10\n",
      "1.5707968257e-10\n"
     ]
    }
   ],
   "source": [
    "x = np.pi * 10 ** -10\n",
    "\n",
    "a = 1 - 1 / np.sqrt(1 + x)\n",
    "b = (np.sqrt(1 + x) - 1)/ np.sqrt(1 + x)\n",
    "\n",
    "print(a)\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\sqrt{1-x}$ を分母，分子に掛けて，$x / (1 + x + \\sqrt{1 + x})$ のように差をうまく消すと，桁落ちを回避しやすい"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.57079632642e-10\n",
      "4.99520067125e-17 4.99273326903e-17\n"
     ]
    }
   ],
   "source": [
    "c = x / (1 + x + np.sqrt(1 + x))\n",
    "print(c)\n",
    "print(np.abs(a - c), np.abs(b - c))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:py3k]",
   "language": "python",
   "name": "conda-env-py3k-py"
  },
  "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}