tutorial/106_linear_regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "accredited-passion",
   "metadata": {},
   "source": [
    "# linear regression\n",
    "\n",
    "Here we describe linear regression using annealing.\n",
    "\n",
    "Reference: https://arxiv.org/abs/2008.02355"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "latest-smooth",
   "metadata": {},
   "source": [
    "### Loss function for linear regression"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "vocal-taxation",
   "metadata": {},
   "source": [
    "Let $\\boldsymbol{x}$ be the data (vector) with $n$ variables, and $y$ be the target value you want to predict for it.  \n",
    "\n",
    "Assuming that $y$ is predictable by the weighted sum of each variable, we can consider the following weight vector $\\boldsymbol{w}$.\n",
    "\n",
    "$$\n",
    "\\boldsymbol{x}^T \\cdot \\boldsymbol{w} = y\n",
    "$$\n",
    "\n",
    "Next, suppose we have a training data set $X$ ($m \\times n$ matrix) consisting of $m$ data and target data $Y$ ($m$-dimensional vector).  \n",
    "The goal of linear regression is to find a single weight vector $\\boldsymbol{w}$ that predicts the corresponding $y$ for all the data $\\boldsymbol{x}$ in $X$.\n",
    "\n",
    "This becomes a function minimization problem as follows.\n",
    "\n",
    "$$\n",
    "\\min_{\\boldsymbol{w}} E(\\boldsymbol{w}) =  || X\\boldsymbol{w}  - Y ||^2\n",
    "$$\n",
    "\n",
    "By transforming $E(\\boldsymbol{w})$\n",
    "\n",
    "$$\n",
    "\\min_{\\boldsymbol{w}} E(\\boldsymbol{w}) = \\boldsymbol{w}^T X^T X \\boldsymbol{w} - 2\\boldsymbol{w}^T X^T Y + Y^T Y\n",
    "$$\n",
    "\n",
    "Let $\\boldsymbol{w} = P\\hat{\\boldsymbol{w}}$ to encode the weights $\\boldsymbol{w}$ into qubit measurement results.  \n",
    "$\\hat{\\boldsymbol{w}}$ is a vector consisting of qubit measurement results $\\{0, 1\\}$.  \n",
    "We also omit $Y^T Y$, which does not affect minimization.\n",
    "\n",
    "$$\n",
    "\\min_{\\boldsymbol{w}} E(\\boldsymbol{w}) = 　\\hat{\\boldsymbol{w}}^T P^T  X^T X P\\hat{\\boldsymbol{w}} - 2\\hat{\\boldsymbol{w}}^T P^T X^T Y\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "radio-webcam",
   "metadata": {},
   "source": [
    "Then it can be reduced to a general QUBO problem, which can be solved by annealing.\n",
    "\n",
    "$$\n",
    "\\min_{\\boldsymbol{w}} E(\\boldsymbol{w}) = 　\\hat{\\boldsymbol{w}}^T A \\hat{\\boldsymbol{w}} + \\hat{\\boldsymbol{w}}^T \\boldsymbol{b}\n",
    "$$\n",
    "\n",
    "where \n",
    "\n",
    "$$\n",
    "A = P^T  X^T X P\n",
    "$$\n",
    "\n",
    "and \n",
    "\n",
    "$$\n",
    "\\boldsymbol{b} = -2P^T X^T Y\n",
    "$$\n",
    "\n",
    "It is worth noting that the size of the QUBO matrix does not depend on the number of data contained in the data set.  \n",
    "The size is  (the number of variables in the data) $\\times$ (the number of qubits used to encode the weight values).  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "effective-october",
   "metadata": {},
   "source": [
    "Let's implement this using blueqat."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "olive-nursery",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "worth-packaging",
   "metadata": {},
   "source": [
    "The first step is to create a data set.  \n",
    "This time, we will set the weights we estimate in advance.  \n",
    "Generate the data set $X$ randomly, and calculate the weighted and noise-added value to be the target data $Y$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "manual-breath",
   "metadata": {},
   "outputs": [],
   "source": [
    "w = np.array([0.25, 0.75, 0.5])\n",
    "X = np.random.rand(100, 3)\n",
    "\n",
    "y = X @ w + np.random.normal(scale = 0.05, size = X.shape[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "perceived-production",
   "metadata": {},
   "source": [
    "Let's do a classical linear regression using scikit-learn.  \n",
    "Although there is a slight deviation due to the addition of noise to the target data $Y$, the predetermined weights can be estimated with good accuracy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "sunset-enemy",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Predicted weight: [0.22016659 0.73691686 0.53090383]\n",
      "True weight: [0.25 0.75 0.5 ]\n"
     ]
    }
   ],
   "source": [
    "from sklearn import linear_model\n",
    "\n",
    "skmodel = linear_model.LinearRegression()\n",
    "skmodel.fit(X, y)\n",
    "w_sk = skmodel.coef_\n",
    "print(\"Predicted weight:\", w_sk)\n",
    "print(\"True weight:\", w)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "light-miami",
   "metadata": {},
   "source": [
    "QAOA encodes each value of the weights $\\boldsymbol{w}$ into a number of qubits.  \n",
    "What value to encode can be set arbitrarily with the transformation matrix $P$ in the above formula.  \n",
    "\n",
    "Here we simply encode one weight parameter into two qubits and predict it from four values of $[0, 0.25, 0.5, 0.75]$.  \n",
    "The number of qubits to use is $(\\text{number of qubits used for encoding}) \\times (\\text{number of variables}) = 2 \\times 3 = 6$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "guilty-samba",
   "metadata": {},
   "outputs": [],
   "source": [
    "K = 2 # bit number for weight\n",
    "d = 3 # Number of features\n",
    "\n",
    "p = [2 ** (-i-1) for i in range(K)]\n",
    "I = np.eye(d)\n",
    "P = np.kron(I, p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "median-vintage",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.dot(np.dot(P.T, X.T), np.dot(X, P))\n",
    "b = -2 * np.dot(np.dot(P.T, X.T), y)\n",
    "QUBO = A + np.diag(b)\n",
    "QUBO = np.triu(QUBO + QUBO.T) - np.eye(QUBO.shape[0]) * QUBO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "539a5e47-4d25-42c6-9561-dacab14f3efe",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "011110\n"
     ]
    }
   ],
   "source": [
    "from blueqat.utils import qaoa\n",
    "from blueqat.pauli import from_qubo\n",
    "\n",
    "step = 2\n",
    "hamiltonian = from_qubo(QUBO)\n",
    "\n",
    "result = qaoa(hamiltonian, step)\n",
    "b = result.circuit.run(shots=10)\n",
    "sample = b.most_common(1)[0][0]\n",
    "print(sample)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ruled-explanation",
   "metadata": {},
   "source": [
    "Convert the result of annealing into the value of the weight parameter and compare it with the value that was originally set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "adapted-companion",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Predicted weight: [0.25 0.75 0.5 ]\n",
      "True weight: [0.25 0.75 0.5 ]\n"
     ]
    }
   ],
   "source": [
    "w_qa = P @ [int(i) for i in list(sample)]\n",
    "\n",
    "print(\"Predicted weight:\", w_qa)\n",
    "print(\"True weight:\", w)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "excited-ghana",
   "metadata": {},
   "source": [
    "It was correctly predicted.  \n",
    "Here, the encoding is determined according to the value of the originally set weight $\\boldsymbol{w}$, so it can be estimated exactly.  \n",
    "\n",
    "In reality, there is an error due to the encoding.  \n",
    "It is necessary to prepare more qubits to reduce the encoding error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "a6630e75-6623-46fa-9367-bad707e61b65",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Counter({'110110': 1,\n",
       "         '111110': 1,\n",
       "         '011110': 3,\n",
       "         '111111': 1,\n",
       "         '011111': 3,\n",
       "         '101110': 1})"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}