_programs/python/logistic/logistic_regression_02.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "import struct\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MNISTのファイル (あらかじめダウンロードしておく)\n",
    "train_image_file = 'mnist/train-images-idx3-ubyte'\n",
    "train_label_file = 'mnist/train-labels-idx1-ubyte'\n",
    "test_image_file = 'mnist/t10k-images-idx3-ubyte'\n",
    "test_label_file = 'mnist/t10k-labels-idx1-ubyte'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_images(filename):\n",
    "    \"\"\" MNISTの画像データを読み込む \"\"\"\n",
    "\n",
    "    fp = open(filename, 'rb')\n",
    "    \n",
    "    # マジックナンバー\n",
    "    magic = struct.unpack('>i', fp.read(4))[0]\n",
    "    if magic != 2051:\n",
    "        raise Exception('Invalid MNIST file!')\n",
    "        \n",
    "    # 各種サイズ\n",
    "    n_images, height, width = struct.unpack('>iii', fp.read(4 * 3))\n",
    "    \n",
    "    # 画像の読み込み\n",
    "    total_pixels = n_images * height * width\n",
    "    images = struct.unpack('>' + 'B' * total_pixels, fp.read(total_pixels))\n",
    "    \n",
    "    images = np.asarray(images, dtype='uint8')\n",
    "    images = images.reshape((n_images, height, width, 1))\n",
    "    \n",
    "    # 値の範囲を[0, 1]に変更する\n",
    "    images = images.astype('float32') / 255.0\n",
    "    \n",
    "    fp.close()\n",
    "    \n",
    "    return images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_labels(filename):\n",
    "    \"\"\" MNISTのラベルデータを読み込む \"\"\"\n",
    "\n",
    "    fp = open(filename, 'rb')\n",
    "    \n",
    "    # マジックナンバー\n",
    "    magic = struct.unpack('>i', fp.read(4))[0]\n",
    "    if magic != 2049:\n",
    "        raise Exception('Invalid MNIST file!')\n",
    "        \n",
    "    # 各種サイズ\n",
    "    n_labels = struct.unpack('>i', fp.read(4))[0]\n",
    "    \n",
    "    # ラベルの読み込み\n",
    "    labels = struct.unpack('>' + 'B' * n_labels, fp.read(n_labels))\n",
    "    labels = np.asarray(labels, dtype='int32')\n",
    "    \n",
    "    fp.close()\n",
    "    \n",
    "    return labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def to_onehot(labels):\n",
    "    \"\"\" one-hot形式への変換 \"\"\"\n",
    "    return np.identity(10)[labels]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## モデルの訓練"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "images = load_images(train_image_file)\n",
    "labels = load_labels(train_label_file)\n",
    "onehot = to_onehot(labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_data = len(images)\n",
    "X = images.reshape((num_data, -1))\n",
    "y = onehot.reshape((num_data, -1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def softmax(x, axis=-1):\n",
    "    \"\"\" softmax関数 \"\"\"\n",
    "    x = x - np.max(x, axis=axis, keepdims=True)\n",
    "    ex = np.exp(x)\n",
    "    return ex / np.sum(ex, axis=axis, keepdims=True)\n",
    "\n",
    "def logsumexp(x, axis=-1):\n",
    "    \"\"\" logsumexp関数 \"\"\"\n",
    "    x_max = np.max(x, axis=axis, keepdims=True)\n",
    "    x = x - x_max\n",
    "    return x_max + np.log(np.sum(np.exp(x), axis=axis, keepdims=True))\n",
    "\n",
    "def log_softmax(x, axis=-1):\n",
    "    \"\"\" log-softmax関数 \"\"\"\n",
    "    ex = np.exp(x)\n",
    "    return x - logsumexp(x, axis=axis)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch:1, step:60000, loss=0.615525, acc=0.843750\n",
      "epoch:2, step:60000, loss=0.513534, acc=0.843750\n",
      "epoch:3, step:60000, loss=0.460950, acc=0.906250\n",
      "epoch:4, step:60000, loss=0.269180, acc=0.968750\n",
      "epoch:5, step:60000, loss=0.650967, acc=0.875000\n",
      "epoch:6, step:60000, loss=0.218641, acc=0.906250\n"
     ]
    }
   ],
   "source": [
    "# ミニバッチのサイズ\n",
    "epochs = 6\n",
    "batch_size = 32\n",
    "\n",
    "# モメンタムSGDのパラメータ\n",
    "alpha = 0.01\n",
    "momentum = 0.9\n",
    "\n",
    "# パラメータの初期化\n",
    "in_features = X.shape[-1]\n",
    "out_features = y.shape[-1]\n",
    "AA = np.random.normal(size=(out_features, in_features)) / np.sqrt(0.5 * (out_features + in_features))\n",
    "bb = np.zeros((out_features))\n",
    "\n",
    "# エポック\n",
    "grad_AA = np.zeros_like(AA)\n",
    "grad_bb = np.zeros_like(bb)\n",
    "for epoch in range(epochs):\n",
    "    # データの順番は偏りをなくすためにランダムシャッフルする\n",
    "    indices = np.random.permutation(np.arange(num_data))\n",
    "    for b in range(0, num_data, batch_size):\n",
    "        bs = min(batch_size, num_data - b)        \n",
    "        X_batch = X[indices[b:b+bs], :]\n",
    "        y_batch = y[indices[b:b+bs], :]\n",
    "\n",
    "        loss = 0.0\n",
    "        acc = 0.0\n",
    "        grad_AA_new = np.zeros_like(AA)\n",
    "        grad_bb_new = np.zeros_like(bb)\n",
    "        \n",
    "        # バッチ内の各データに対してロス、精度、勾配を求める\n",
    "        t = np.einsum('ij,bj->bi', AA, X_batch) + bb\n",
    "        log_y_pred = log_softmax(t)\n",
    "        losses = np.sum(-y_batch * log_y_pred, axis=-1, keepdims=True)\n",
    "\n",
    "        delta = np.identity(AA.shape[0])\n",
    "        y_pred = np.exp(log_y_pred)\n",
    "        dLdy = -y_batch / y_pred\n",
    "        dydt = np.einsum('bi,ij->bij', y_pred, delta) - np.einsum('bi,bj->bij', y_pred, y_pred)\n",
    "        dLdt = np.einsum('bi,bij->bj', dLdy, dydt)\n",
    "\n",
    "        dtdA = np.einsum('bk,ij->bijk', X_batch, delta)\n",
    "        dtdb = np.ones((bb.shape[-1], bb.shape[-1]))\n",
    "        dLdA = np.einsum('bi,bijk->bjk', dLdt, dtdA)\n",
    "        dLdb = np.einsum('bi,ij->bj', dLdt, dtdb)\n",
    "            \n",
    "        y_pred_id = np.argmax(y_pred, axis=-1)\n",
    "        y_real_id = np.argmax(y_batch, axis=-1)\n",
    "        acc = np.mean(y_pred_id == y_real_id)\n",
    "        loss = np.mean(losses)\n",
    "\n",
    "        grad_AA_new = np.mean(dLdA, axis=0)\n",
    "        grad_bb_new = np.mean(dLdb, axis=0)\n",
    "        \n",
    "        # 勾配の更新\n",
    "        grad_AA = grad_AA * momentum + grad_AA_new * (1.0 - momentum)\n",
    "        grad_bb = grad_bb * momentum + grad_bb_new * (1.0 - momentum)\n",
    "\n",
    "        # 最急降下法による値の更新\n",
    "        AA -= alpha * grad_AA\n",
    "        bb -= alpha * grad_bb\n",
    "        \n",
    "        # printの代わりにsys.stdout.writeを使うとcarrige returnが使える\n",
    "        sys.stdout.write('\\repoch:{}, step:{}, loss={:.6f}, acc={:.6f}'.format(epoch + 1, b + bs, loss, acc))\n",
    "        sys.stdout.flush()\n",
    "\n",
    "    sys.stdout.write('\\n')\n",
    "    sys.stdout.flush()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## モデルのテスト"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_images = load_images(test_image_file)\n",
    "test_labels = load_labels(test_label_file)\n",
    "test_onehot = to_onehot(test_labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_data = len(images)\n",
    "X = images.reshape((num_data, -1))\n",
    "y = onehot.reshape((num_data, -1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "loss = 0.0\n",
    "acc = 0.0\n",
    "for b in range(0, num_data, batch_size):\n",
    "    bs = min(batch_size, num_data - b)        \n",
    "    X_batch = X[indices[b:b+bs], :]\n",
    "    y_batch = y[indices[b:b+bs], :]\n",
    "    \n",
    "    t = np.einsum('ij,bj->bi', AA, X_batch) + bb\n",
    "    log_y_pred = log_softmax(t)\n",
    "    losses = np.sum(-y_batch * log_y_pred, axis=-1, keepdims=True)\n",
    "\n",
    "    y_pred = np.exp(log_y_pred)\n",
    "    y_pred_id = np.argmax(y_pred, axis=-1)\n",
    "    y_real_id = np.argmax(y_batch, axis=-1)\n",
    "\n",
    "    acc += np.sum(y_pred_id == y_real_id)\n",
    "    loss += np.sum(losses)\n",
    "\n",
    "loss /= num_data\n",
    "acc /= num_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loss: 0.3474\n",
      "Accuracy: 0.9038\n"
     ]
    }
   ],
   "source": [
    "print('Loss: {:.4f}'.format(loss))\n",
    "print('Accuracy: {:.4f}'.format(acc))"
   ]
  }
 ],
 "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.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}