2_equiangular_lines_2n_2.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "equiangular_lines_2n-2.ipynb",
      "provenance": [],
      "authorship_tag": "ABX9TyMu0EH6IxhS7bxFwzdQ7t1W",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/vinayprabhu/Grassmannian_tutorials/blob/master/2_equiangular_lines_2n_2.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tjJ9Uweq-e62",
        "colab_type": "text"
      },
      "source": [
        "# This notebook walks you through the procedure to generate $2n-2$ equiangular line packing for $n$-dimensional spaces and angle of separation $cos^{-1}(1/3)=70.529^{\\circ}$.\n",
        "\n",
        "For simplicty of visualization of the stages, we choose $n=5$"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "3EKlS7DDrmwG",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        },
        "outputId": "8a132186-9f15-4089-d88f-b7c538665f6c"
      },
      "source": [
        "import numpy as np\n",
        "from scipy.linalg import block_diag\n",
        "# Don't do linear algebra in Python without these two lines\n",
        "np.set_printoptions(suppress=True)\n",
        "%precision 3"
      ],
      "execution_count": 1,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "'%.3f'"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 1
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "hKXXd85y_CDZ",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Step-1: Generate the Grammian matrix"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "urJMCYLttcLn",
        "colab_type": "code",
        "outputId": "acc3e326-7941-4c64-84f5-c2f37b9f163f",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 153
        }
      },
      "source": [
        "M=np.array([[1,-1/3],\n",
        "            [-1/3,1]])\n",
        "G=block_diag(M,M,M,M)\n",
        "G[G==0]=1/3\n",
        "G"
      ],
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([[ 1.   , -0.333,  0.333,  0.333,  0.333,  0.333,  0.333,  0.333],\n",
              "       [-0.333,  1.   ,  0.333,  0.333,  0.333,  0.333,  0.333,  0.333],\n",
              "       [ 0.333,  0.333,  1.   , -0.333,  0.333,  0.333,  0.333,  0.333],\n",
              "       [ 0.333,  0.333, -0.333,  1.   ,  0.333,  0.333,  0.333,  0.333],\n",
              "       [ 0.333,  0.333,  0.333,  0.333,  1.   , -0.333,  0.333,  0.333],\n",
              "       [ 0.333,  0.333,  0.333,  0.333, -0.333,  1.   ,  0.333,  0.333],\n",
              "       [ 0.333,  0.333,  0.333,  0.333,  0.333,  0.333,  1.   , -0.333],\n",
              "       [ 0.333,  0.333,  0.333,  0.333,  0.333,  0.333, -0.333,  1.   ]])"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 3
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QRF8c3BR_Qyx",
        "colab_type": "text"
      },
      "source": [
        "Verify that this is a rank deficient matrix whose rank is $2n-2=8$"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_DUxcCdnt3z3",
        "colab_type": "code",
        "outputId": "6710bbff-b404-4792-f2e7-2465c1b0529b",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        }
      },
      "source": [
        "np.linalg.matrix_rank(G)"
      ],
      "execution_count": 4,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "5"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 4
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QvBGUO3w_atv",
        "colab_type": "text"
      },
      "source": [
        "# Step-2: Perform the Jordan decomposition of the grammian and extract the codebook (F). \n",
        "The grammian, by definition is $G=FF^{T}$.\n",
        "So, the easiest way to obtain the codebook is perform Jordan decomposition ($G=PDP^{-1}$) that yields the codebook to be: $F=PD^{1/2}$"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "QMyzGfPskFDV",
        "colab_type": "code",
        "outputId": "c71ba8c5-b53d-498b-da91-03d25edfdfc7",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 306
        }
      },
      "source": [
        "D,P=np.linalg.eigh(G) # Don't use np.linalg.eig ... aaargh. I hate Python!\n",
        "P,D, P.T-np.linalg.inv(P)"
      ],
      "execution_count": 5,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(array([[ 0.088,  0.606,  0.   ,  0.   ,  0.697,  0.   , -0.117,  0.354],\n",
              "        [ 0.088,  0.606,  0.   ,  0.   , -0.697,  0.   ,  0.117,  0.354],\n",
              "        [ 0.542, -0.285, -0.   ,  0.   , -0.117,  0.   , -0.697,  0.354],\n",
              "        [ 0.542, -0.285, -0.   ,  0.   ,  0.117,  0.   ,  0.697,  0.354],\n",
              "        [-0.315, -0.16 , -0.5  ,  0.707, -0.   ,  0.002, -0.   ,  0.354],\n",
              "        [-0.315, -0.16 , -0.5  , -0.707, -0.   , -0.002,  0.   ,  0.354],\n",
              "        [-0.315, -0.16 ,  0.5  ,  0.002, -0.   , -0.707,  0.   ,  0.354],\n",
              "        [-0.315, -0.16 ,  0.5  , -0.002, -0.   ,  0.707,  0.   ,  0.354]]),\n",
              " array([-0.   , -0.   ,  0.   ,  1.333,  1.333,  1.333,  1.333,  2.667]),\n",
              " array([[ 0.,  0., -0., -0.,  0.,  0.,  0.,  0.],\n",
              "        [ 0.,  0.,  0.,  0., -0., -0.,  0.,  0.],\n",
              "        [-0.,  0.,  0., -0.,  0.,  0.,  0.,  0.],\n",
              "        [ 0.,  0.,  0.,  0., -0.,  0.,  0.,  0.],\n",
              "        [ 0.,  0., -0., -0., -0., -0., -0., -0.],\n",
              "        [-0.,  0.,  0.,  0.,  0.,  0., -0., -0.],\n",
              "        [ 0., -0.,  0., -0.,  0.,  0.,  0.,  0.],\n",
              "        [ 0.,  0.,  0.,  0., -0., -0.,  0.,  0.]]))"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 5
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "M8HPdA3L1ZrS",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 153
        },
        "outputId": "e22880be-78ff-4738-bf09-c42604a625ed"
      },
      "source": [
        "D[np.isclose(D,0)]=0 # For numerical stability \n",
        "F=P @ np.sqrt(np.diag(D))\n",
        "F"
      ],
      "execution_count": 6,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([[ 0.   ,  0.   ,  0.   ,  0.   ,  0.805,  0.   , -0.135,  0.577],\n",
              "       [ 0.   ,  0.   ,  0.   ,  0.   , -0.805,  0.   ,  0.135,  0.577],\n",
              "       [ 0.   ,  0.   ,  0.   ,  0.   , -0.135,  0.   , -0.805,  0.577],\n",
              "       [ 0.   ,  0.   ,  0.   ,  0.   ,  0.135,  0.   ,  0.805,  0.577],\n",
              "       [ 0.   ,  0.   ,  0.   ,  0.816, -0.   ,  0.003, -0.   ,  0.577],\n",
              "       [ 0.   ,  0.   ,  0.   , -0.816, -0.   , -0.003,  0.   ,  0.577],\n",
              "       [ 0.   ,  0.   ,  0.   ,  0.003, -0.   , -0.816,  0.   ,  0.577],\n",
              "       [ 0.   ,  0.   ,  0.   , -0.003, -0.   ,  0.816,  0.   ,  0.577]])"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 6
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "5FniPyHOxlDQ",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 374
        },
        "outputId": "dcf567c1-3a55-40e0-be45-3ef5c3ee3de9"
      },
      "source": [
        "P_5=np.zeros((8,5))\n",
        "D_5=np.zeros((5,5))\n",
        "j=0\n",
        "for i in range(8):\n",
        "  if (~ np.isclose(D[i],0)):\n",
        "    P_5[:,j]=P[:,i]\n",
        "    D_5[j,j]=D[i]\n",
        "    j+=1\n",
        "F_5=P_5 @ np.sqrt(D_5)\n",
        "P_5,D_5,F_5"
      ],
      "execution_count": 7,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(array([[ 0.   ,  0.697,  0.   , -0.117,  0.354],\n",
              "        [ 0.   , -0.697,  0.   ,  0.117,  0.354],\n",
              "        [ 0.   , -0.117,  0.   , -0.697,  0.354],\n",
              "        [ 0.   ,  0.117,  0.   ,  0.697,  0.354],\n",
              "        [ 0.707, -0.   ,  0.002, -0.   ,  0.354],\n",
              "        [-0.707, -0.   , -0.002,  0.   ,  0.354],\n",
              "        [ 0.002, -0.   , -0.707,  0.   ,  0.354],\n",
              "        [-0.002, -0.   ,  0.707,  0.   ,  0.354]]),\n",
              " array([[1.333, 0.   , 0.   , 0.   , 0.   ],\n",
              "        [0.   , 1.333, 0.   , 0.   , 0.   ],\n",
              "        [0.   , 0.   , 1.333, 0.   , 0.   ],\n",
              "        [0.   , 0.   , 0.   , 1.333, 0.   ],\n",
              "        [0.   , 0.   , 0.   , 0.   , 2.667]]),\n",
              " array([[ 0.   ,  0.805,  0.   , -0.135,  0.577],\n",
              "        [ 0.   , -0.805,  0.   ,  0.135,  0.577],\n",
              "        [ 0.   , -0.135,  0.   , -0.805,  0.577],\n",
              "        [ 0.   ,  0.135,  0.   ,  0.805,  0.577],\n",
              "        [ 0.816, -0.   ,  0.003, -0.   ,  0.577],\n",
              "        [-0.816, -0.   , -0.003,  0.   ,  0.577],\n",
              "        [ 0.003, -0.   , -0.816,  0.   ,  0.577],\n",
              "        [-0.003, -0.   ,  0.816,  0.   ,  0.577]]))"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 7
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "cAbMXBtcAl8V",
        "colab_type": "text"
      },
      "source": [
        "# Step-3: Verify that the dot products are all of magnitude 1/3 "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "i1_Q3WsZ0X-K",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 153
        },
        "outputId": "264c4e45-d59e-48d0-8bf4-b506468a714a"
      },
      "source": [
        "F_5 @ F_5.T"
      ],
      "execution_count": 8,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([[ 1.   , -0.333,  0.333,  0.333,  0.333,  0.333,  0.333,  0.333],\n",
              "       [-0.333,  1.   ,  0.333,  0.333,  0.333,  0.333,  0.333,  0.333],\n",
              "       [ 0.333,  0.333,  1.   , -0.333,  0.333,  0.333,  0.333,  0.333],\n",
              "       [ 0.333,  0.333, -0.333,  1.   ,  0.333,  0.333,  0.333,  0.333],\n",
              "       [ 0.333,  0.333,  0.333,  0.333,  1.   , -0.333,  0.333,  0.333],\n",
              "       [ 0.333,  0.333,  0.333,  0.333, -0.333,  1.   ,  0.333,  0.333],\n",
              "       [ 0.333,  0.333,  0.333,  0.333,  0.333,  0.333,  1.   , -0.333],\n",
              "       [ 0.333,  0.333,  0.333,  0.333,  0.333,  0.333, -0.333,  1.   ]])"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 8
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "To30q-SxBOQ0",
        "colab_type": "text"
      },
      "source": [
        "# Now, save the codebooks"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "74fOvluEBNH9",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "np.save('emb_8_5.npy',F_5)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "KWzpGMtcC6ip",
        "colab_type": "text"
      },
      "source": [
        "# Laste step: Sew this numerical recipe up as a function!\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "fhzk-MeSGQmq",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def codebook_gen(n):\n",
        "  print(f'About to generate equi-angular line packings in {n}-D euclidean space.\\n  This will yield {2*n-2} equiangular lines with angle 1/3')\n",
        "  # 1: Generate the Grammian!\n",
        "  M=np.array([[1,-1/3],\n",
        "            [-1/3,1]])\n",
        "  n_rep=int(n-1)\n",
        "  G = block_diag(*[M]*n_rep)\n",
        "  G[G==0]=1/3\n",
        "  # 2: Perform rank-check!\n",
        "  if(np.linalg.matrix_rank(G)!=n):\n",
        "    print('Rank of the Grammian not matching!')\n",
        "    F_code=np.empty(0)\n",
        "  else:\n",
        "    # 3: Jordan decomposition and extracting codebooks\n",
        "    D,P=np.linalg.eigh(G)\n",
        "    code_columns=np.where(~ np.isclose(D,0))[0]\n",
        "    P_code=P[:,code_columns]\n",
        "    F_code=P_code @ np.sqrt(np.diag(D[code_columns]))\n",
        "    # 4: Pairwise distance comparisons\n",
        "    dist_mat=F_code @ F_code.T\n",
        "    pairwise_dist=np.abs(dist_mat[np.triu_indices(2*n-2,k=1)])\n",
        "    if(np.allclose(pairwise_dist,1/3)):\n",
        "      file_name=f'embedding_{n}_{2*n-2}.npy'\n",
        "      print(f'All the lines are 1/3 away!\\n   Saving the codebook as {file_name}')\n",
        "      np.save(file_name, F_code)\n",
        "    else:\n",
        "      print('Lines are not 1/3 away :( \\n Something went wrong. No codebooks for you!)')\n",
        "      F_code=np.empty(0)\n",
        "  return F_code\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "G4gY1s0MM41l",
        "colab_type": "text"
      },
      "source": [
        "Now, run it for $n=501$ which yields 1000 lines that can potentially be used for imagenet"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "yqMoHCLKHvXr",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 221
        },
        "outputId": "a2f19bd5-76b6-44ec-9ec7-63757c45d45d"
      },
      "source": [
        "# First, let's verify for the n=5 case\n",
        "x=codebook_gen(5)\n",
        "print(x@x.T)"
      ],
      "execution_count": 11,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "About to generate equi-angular line packings in 5-D euclidean space.\n",
            "  This will yield 8 equiangular lines with angle 1/3\n",
            "All the lines are 1/3 away!\n",
            "   Saving the codebook as embedding_5_8.npy\n",
            "[[ 1.    -0.333  0.333  0.333  0.333  0.333  0.333  0.333]\n",
            " [-0.333  1.     0.333  0.333  0.333  0.333  0.333  0.333]\n",
            " [ 0.333  0.333  1.    -0.333  0.333  0.333  0.333  0.333]\n",
            " [ 0.333  0.333 -0.333  1.     0.333  0.333  0.333  0.333]\n",
            " [ 0.333  0.333  0.333  0.333  1.    -0.333  0.333  0.333]\n",
            " [ 0.333  0.333  0.333  0.333 -0.333  1.     0.333  0.333]\n",
            " [ 0.333  0.333  0.333  0.333  0.333  0.333  1.    -0.333]\n",
            " [ 0.333  0.333  0.333  0.333  0.333  0.333 -0.333  1.   ]]\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "IJ6BEPNUI_oq",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 119
        },
        "outputId": "8e7f651b-bec0-4bd3-ff54-202b21a1d052"
      },
      "source": [
        "%%time\n",
        "x=codebook_gen(501)\n",
        "x@x.T"
      ],
      "execution_count": 12,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "About to generate equi-angular line packings in 501-D euclidean space.\n",
            "  This will yield 1000 equiangular lines with angle 1/3\n",
            "All the lines are 1/3 away!\n",
            "   Saving the codebook as embedding_501_1000.npy\n",
            "CPU times: user 1.26 s, sys: 320 ms, total: 1.58 s\n",
            "Wall time: 850 ms\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "wE9Oor-DN6js",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 119
        },
        "outputId": "d6be9747-95a0-450a-846a-e9bc7c4a5b80"
      },
      "source": [
        "%%time\n",
        "x=codebook_gen(512)\n",
        "x@x.T"
      ],
      "execution_count": 13,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "About to generate equi-angular line packings in 512-D euclidean space.\n",
            "  This will yield 1022 equiangular lines with angle 1/3\n",
            "All the lines are 1/3 away!\n",
            "   Saving the codebook as embedding_512_1022.npy\n",
            "CPU times: user 1.34 s, sys: 246 ms, total: 1.59 s\n",
            "Wall time: 810 ms\n"
          ],
          "name": "stdout"
        }
      ]
    }
  ]
}