03_kernel_random_fourier_empty.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4ad26283",
   "metadata": {},
   "source": [
    "# Random Fourier Features for Scalable Kernels\n",
    "\n",
    "\n",
    "[Random Features for Large-Scale Kernel Machines](https://proceedings.neurips.cc/paper/2007/file/013a006f03dbc5392effeb8f18fda755-Paper.pdf) by Ali Rahimi and Benjamin Recht."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64e7d6e7",
   "metadata": {},
   "source": [
    "[`sklearn.metrics.pairwise.rbf_kernel`](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.pairwise.rbf_kernel.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42f09263",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import ticker\n",
    "import numpy as np\n",
    "from sklearn.metrics.pairwise import rbf_kernel\n",
    "from sklearn.datasets import make_s_curve\n",
    "\n",
    "\n",
    "N    = 1000\n",
    "D    = 3\n",
    "X, t = make_s_curve(N, noise=0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f328995",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_3d(points, points_color, title):\n",
    "    x, y, z = points.T\n",
    "\n",
    "    fig, ax = plt.subplots(\n",
    "        figsize=(6, 6),\n",
    "        facecolor=\"white\",\n",
    "        tight_layout=True,\n",
    "        subplot_kw={\"projection\": \"3d\"},\n",
    "    )\n",
    "    fig.suptitle(title, size=16)\n",
    "    col = ax.scatter(x, y, z, c=points_color, s=50, alpha=0.8)\n",
    "    ax.view_init(azim=-60, elev=9)\n",
    "    ax.xaxis.set_major_locator(ticker.MultipleLocator(1))\n",
    "    ax.yaxis.set_major_locator(ticker.MultipleLocator(1))\n",
    "    ax.zaxis.set_major_locator(ticker.MultipleLocator(1))\n",
    "\n",
    "    fig.colorbar(col, ax=ax, orientation=\"horizontal\", shrink=0.6, aspect=60, pad=0.01)\n",
    "    plt.show()\n",
    "\n",
    "plot_3d(X, t, \"Original S-curve samples\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45310c47",
   "metadata": {},
   "outputs": [],
   "source": [
    "X = X[t.argsort()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fcd72450",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 5)\n",
    "fig.set_size_inches(15, 4)\n",
    "font = {'fontname': 'serif', 'fontsize': 18}\n",
    "\n",
    "gamma = 1/2.\n",
    "# The RBF kernel is the Gaussian kernel.\n",
    "K = rbf_kernel(X, gamma=gamma)\n",
    "\n",
    "axes[0].imshow(K, cmap=plt.cm.Blues)\n",
    "axes[0].set_title('Exact RBF kernel', **font)\n",
    "axes[0].set_xticks([])\n",
    "axes[0].set_yticks([])\n",
    "\n",
    "sigma = 1.\n",
    "\n",
    "for R, ax in zip([1, 10, 100, 1000], axes[1:]):\n",
    "    W    = np.random.normal(loc=0, scale=1, size=(R, D))\n",
    "    B    = np.random.uniform(0, 2*np.pi, size=(R,1))\n",
    "    norm = 1./ np.sqrt(R)\n",
    "    Z    = norm * np.sqrt(2) * np.cos(sigma*W @ X.T + B)\n",
    "    ZZ   = Z.T@Z\n",
    "\n",
    "    ax.imshow(ZZ, cmap=plt.cm.Blues)\n",
    "    ax.set_title(r'$\\mathbf{Z} \\mathbf{Z}^{\\top}$, $R=%s$' % R, **font)\n",
    "    ax.set_xticks([])\n",
    "    ax.set_yticks([])\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8d241cf",
   "metadata": {},
   "source": [
    "Radial basis function kernel [`RBF`](https://scikit-learn.org/stable/modules/generated/sklearn.gaussian_process.kernels.RBF.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7d024c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "from   sklearn.gaussian_process.kernels import RBF\n",
    "\n",
    "N     = 100\n",
    "X     = np.linspace(-10, 10, 2*N)[:, None]\n",
    "\n",
    "mean  = np.zeros(X.shape[0])\n",
    "cov   = RBF()(X)\n",
    "y_true = np.random.multivariate_normal(mean, cov)\n",
    "\n",
    "noise = np.random.normal(0, 0.5, 2*N)\n",
    "y    = y_true + noise\n",
    "\n",
    "X_train =  X[::2]#[:, None]\n",
    "y_train = y[::2]#[:, None]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1f8be0c",
   "metadata": {},
   "outputs": [],
   "source": [
    "cmap      = plt.cm.get_cmap('Blues')\n",
    "plt.plot(X, y_true, c=cmap(1.))\n",
    "plt.scatter(X_train, y_train, s=30, c=[cmap(0.3)])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0a3ae26",
   "metadata": {},
   "source": [
    "[`KernelRidge`](https://scikit-learn.org/stable/modules/generated/sklearn.kernel_ridge.KernelRidge.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4ef7d5d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "from   sklearn.kernel_ridge import KernelRidge\n",
    "\n",
    "# Fit kernel ridge regression using an RBF kernel.\n",
    "\n",
    "y_pred = ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21aa4ea5",
   "metadata": {},
   "outputs": [],
   "source": [
    "cmap      = plt.cm.get_cmap('Blues')\n",
    "plt.plot(X, y_true, c=cmap(1.))\n",
    "plt.scatter(X, y, s=30, c=[cmap(0.3)])\n",
    "plt.plot(X, y_pred, c='red')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "692a3b88",
   "metadata": {},
   "source": [
    "[`Ridge`](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.Ridge.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e1475cf7",
   "metadata": {},
   "outputs": [],
   "source": [
    "from   sklearn.linear_model import Ridge\n",
    "from   sklearn.exceptions import NotFittedError\n",
    "\n",
    "class RFFRidgeRegression:\n",
    "\n",
    "    def __init__(self, rff_dim=1, alpha=1.0, sigma=1.0):\n",
    "        \"\"\"Kernel ridge regression using random Fourier features.\n",
    "        rff_dim : Dimension of random feature.\n",
    "        alpha :   Regularization strength. Should be a positive float.\n",
    "        \"\"\"\n",
    "        self.fitted  = False\n",
    "        self.rff_dim = rff_dim\n",
    "        self.sigma   = sigma\n",
    "        self.lm      = Ridge(alpha=alpha)\n",
    "        self.B_      = None\n",
    "        self.W_      = None\n",
    "\n",
    "    def fit(self, X, y):\n",
    "        \"\"\"Fit model with training data X and target y.\n",
    "        \"\"\"\n",
    "        Z = self._get_rffs(X)\n",
    "        self.lm.fit(Z.T, y)\n",
    "        self.fitted = True\n",
    "        pass\n",
    "\n",
    "    def predict(self, X):\n",
    "        \"\"\"Predict using fitted model and testing data X.\n",
    "        \"\"\"\n",
    "        if not self.fitted:\n",
    "            msg = \"Call 'fit' with appropriate arguments first.\"\n",
    "            raise NotFittedError(msg)\n",
    "        Z = self._get_rffs(X)\n",
    "        return self.lm.predict(Z.T)\n",
    "\n",
    "    def _get_rffs(self, X):\n",
    "        \"\"\"Return random Fourier features based on data X, as well as random\n",
    "        variables W and b.\n",
    "        \"\"\"\n",
    "        # your code here\n",
    "        return Z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a3d2e67",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up figure and plot data.\n",
    "fig, axes = plt.subplots(2, 1)\n",
    "fig.set_size_inches(10, 5)\n",
    "ax1, ax2  = axes\n",
    "cmap      = plt.cm.get_cmap('Blues')\n",
    "\n",
    "ax1.scatter(X, y, s=30, c=[cmap(0.3)])\n",
    "ax2.scatter(X, y, s=30, c=[cmap(0.3)])\n",
    "ax1.plot(X, y_true, c=cmap(1.))\n",
    "ax2.plot(X, y_true, c=cmap(1.))\n",
    "\n",
    "# Fit kernel ridge regression using an RBF kernel.\n",
    "clf    = KernelRidge(kernel=RBF())\n",
    "clf    = clf.fit(X_train, y_train)\n",
    "y_pred = clf.predict(X)\n",
    "ax1.plot(X, y_pred, c='red')\n",
    "\n",
    "# Fit kernel ridge regression using random Fourier features.\n",
    "rff_dim = 20\n",
    "clf     = RFFRidgeRegression(rff_dim=rff_dim)\n",
    "clf.fit(X_train, y_train)\n",
    "y_pred  = clf.predict(X)\n",
    "ax2.plot(X, y_pred, c='red')\n",
    "\n",
    "# Labels, etc.\n",
    "ax1.margins(0, 0.1)\n",
    "ax1.set_title('RBF kernel regression')\n",
    "ax1.set_ylabel(r'$y$', fontsize=14)\n",
    "ax1.set_xticks([])\n",
    "ax2.margins(0, 0.1)\n",
    "ax2.set_title(rf'RFF ridge regression, $R = {rff_dim}$')\n",
    "ax2.set_ylabel(r'$y$', fontsize=14)\n",
    "ax2.set_xlabel(r'$x$', fontsize=14)\n",
    "ax2.set_xticks(np.arange(-10, 10.1, 1))\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e48f46c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "import torchvision\n",
    "from torchvision import datasets\n",
    "import os\n",
    "from pathlib import Path\n",
    "\n",
    "ROOT_DIR = Path.home()\n",
    "data_path = os.path.join(ROOT_DIR,'data/MNIST/')\n",
    "#root_dir = './data/MNIST/'\n",
    "#torchvision.datasets.MNIST(root=data_path,download=True)\n",
    "\n",
    "train_set = torchvision.datasets.MNIST(root=data_path, train=True, download=True)\n",
    "test_set = torchvision.datasets.MNIST(root=data_path, train=False, download=True)\n",
    "\n",
    "def MNIST_numpy(dataset):\n",
    "    images = dataset.data.numpy().astype(np.float32)/255-0.5\n",
    "    labels = dataset.targets.numpy()\n",
    "    return images.reshape((-1,28*28)), labels\n",
    "\n",
    "MNIST_train, y_train = MNIST_numpy(train_set)\n",
    "MNIST_test, y_test = MNIST_numpy(test_set)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1d496ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(MNIST_train.shape,y_train.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "efcf4e53",
   "metadata": {},
   "outputs": [],
   "source": [
    "def mat_transform_pca(Xs, dim=128):\n",
    "    _, V = np.linalg.eig(Xs.T.dot(Xs))\n",
    "    return np.real(V[:, :dim])\n",
    "\n",
    "T = mat_transform_pca(MNIST_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ebf2ac8",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train = MNIST_train.dot(T)\n",
    "X_test = MNIST_test.dot(T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cfc69cf3",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(X_test.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a991f94",
   "metadata": {},
   "source": [
    "[Support Vector Machines](https://scikit-learn.org/stable/modules/svm.html)\n",
    "\n",
    "[`SVC`](https://scikit-learn.org/stable/modules/generated/sklearn.svm.SVC.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7efb93a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.svm import SVC\n",
    "\n",
    "svc = SVC(kernel = \"rbf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be848df8",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "svc.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5119616",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "score = 100 * svc.score(X_test, y_test)\n",
    "score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fdf74d63",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.svm import LinearSVC\n",
    "\n",
    "class RFFClassif:\n",
    "\n",
    "    def __init__(self, rff_dim=1, sigma=1.0):\n",
    "        \"\"\"Classification using random Fourier features.\n",
    "        rff_dim : Dimension of random feature.\n",
    "        \"\"\"\n",
    "        self.fitted  = False\n",
    "        self.rff_dim = rff_dim\n",
    "        self.sigma   = sigma\n",
    "        self.lm      = LinearSVC()\n",
    "        self.B_      = None\n",
    "        self.W_      = None\n",
    "\n",
    "    # your code here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "396af65f",
   "metadata": {},
   "outputs": [],
   "source": [
    "rffsvc = RFFClassif(rff_dim=640, sigma=0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f2e1d21",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "rffsvc.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d8ee397",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "rffsvc.score(X_test, y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ebbb60dd",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ml-icfp",
   "language": "python",
   "name": "ml-icfp"
  },
  "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}