Random_Transform.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.linalg import orth\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Fast SVD approximation using the Moore-Penrose pseudoinverse !\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 244,
   "metadata": {},
   "outputs": [],
   "source": [
    "# nxm matrix\n",
    "n=500\n",
    "m=800\n",
    "k=10 # maximum number of singular values we use for estimation\n",
    "A = np.random.rand(m,n)\n",
    "U,S,V = np.linalg.svd(A,full_matrices=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 245,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(800, 500)\n",
      "(500,)\n",
      "(500, 500)\n"
     ]
    }
   ],
   "source": [
    "print(U.shape)\n",
    "print(S.shape)\n",
    "print(V.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 246,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(500, 500)\n",
      "(800, 500)\n"
     ]
    }
   ],
   "source": [
    "S[k+1:] = S[k+1:]*(np.arange(k+1,min(n,m),dtype=float).conjugate().T)[:]**-1 #  is the numpy complex conjugate operator\n",
    "S = np.diag(S,)\n",
    "print(S.shape)\n",
    "A = U@S@V.conjugate().T\n",
    "print(A.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 247,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(500,)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(S[k+1].shape)\n",
    "print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 248,
   "metadata": {},
   "outputs": [],
   "source": [
    "def random_sketch_approx(A,k_over):\n",
    "    n = min(A.shape)\n",
    "    if A.shape[0] > A.shape[1]:\n",
    "        min_ax = 0\n",
    "    else:\n",
    "        min_ax = 1\n",
    "    k_over = 2*k # The fix? Take k+10\n",
    "    ds = np.exp(1j*2*np.pi*np.random.rand(n,1)) # generate random gausian vectors\n",
    "    sm = np.random.randint(0,high=n,size=k_over) \n",
    "    \n",
    "    a_sketch = A.T\n",
    "    a_sketch = a_sketch*ds\n",
    "    a_sketch = np.fft.fft(a_sketch,n,axis=min_ax)/np.sqrt(n) # n point fft over the longer axis\n",
    "    a_sketch = a_sketch[sm,:]\n",
    "    a_sketch = a_sketch.T\n",
    "\n",
    "    q = orth(a_sketch)\n",
    "    A_small = q.conjugate().T@A\n",
    "    [uu,ss,vv] = np.linalg.svd(A_small,full_matrices=False)\n",
    "\n",
    "    return [uu,ss,vv], q\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 249,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Absolute Error (norm):  4.117761511988246\n",
      "Relative error: [8.66681750e-06 1.74732326e-03 2.87798434e-03 1.65844722e-03\n",
      " 2.26399367e-03 2.07205488e-03 2.51777592e-03 2.49454076e-03\n",
      " 3.76928998e-03 4.05329820e-03]\n"
     ]
    }
   ],
   "source": [
    "[uu,ss,vv], q = random_sketch_approx(A,15)\n",
    "\n",
    "A_error = np.real(q@q.conjugate().T@A)\n",
    "norm = np.linalg.norm(A-A_error)\n",
    "print(\"Absolute Error (norm): \",norm)\n",
    "sigs_exact = np.diag(S)\n",
    "print(\"Relative error:\",(sigs_exact[:k]-ss[:k])/sigs_exact[:k])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Time Trials \n",
    "## Numpy svd (np.linalg.svd(A))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 250,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "50.3 ms ± 1.94 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "U,S,V = np.linalg.svd(A,full_matrices=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numpy svd (np.linalg.svd(A))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 251,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "17 ms ± 2.9 ms per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "[uu,ss,vv], q = random_sketch_approx(A,10)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}