randomized_numerical_linear_algebra.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook demonstrates the idea of approximating matrix multiplication using sampling. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "metadata": {},
   "outputs": [],
   "source": [
    "def norm_sq_distribution(A, B):\n",
    "    if A.shape[1] != B.shape[0]:\n",
    "        raise ValueError(\"dimensions don't match\")\n",
    "    Pr_A = np.linalg.norm(A, axis = 0) # norm of cols of A\n",
    "    Pr_B = np.linalg.norm(B, axis = 1) # norm of rows of B\n",
    "    Pr = np.multiply(Pr_A, Pr_B)  # element wise multiplication\n",
    "    return Pr/(np.sum(Pr))  \n",
    "\n",
    "def approx_multiply(A, B, s):\n",
    "    m, n = A.shape\n",
    "    n_, p = B.shape\n",
    "    if n != n_:\n",
    "        raise ValueError(\"dimensions don't match\")\n",
    "    Pr = norm_sq_distribution(A, B)\n",
    "    C = np.zeros((m, s))\n",
    "    R = np.zeros((s, p))\n",
    "    for i in range(s):\n",
    "        j = np.random.choice(n, p = Pr)\n",
    "        C[:, i] = A[:, j]/((s*Pr[j])**0.5)\n",
    "        R[i, :] = B[j, :]/((s*Pr[j])**0.5)\n",
    "    return C, R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7328.2854237891197"
      ]
     },
     "execution_count": 125,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.random.randn(1000,50)\n",
    "np.linalg.norm(A @ A.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 150,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10282.836637106488"
      ]
     },
     "execution_count": 150,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C, R = approx_multiply(A, A.T, 70)\n",
    "np.linalg.norm(C @ R)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Approximaion bound:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [],
   "source": [
    "def col_norm_sq_distribution(A):\n",
    "    n = A.shape[1] # number of cols\n",
    "    Pr = np.linalg.norm(A, axis = 0)**2 \n",
    "    return Pr/(np.sum(Pr))  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}