psets/pset2/PS2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem set 2 (145 pts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 1 (LU decomposition) 25 pts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. LU for band matrices (7 pts)\n",
    "\n",
    "The complexity to find an LU decomposition of a dense $n\\times n$ matrix is $\\mathcal{O}(n^3)$.\n",
    "Significant reduction in complexity can be achieved if the matrix has a certain structure, e.g. it is sparse. \n",
    "In the following task we consider an important example of $LU$ for a special type of sparse matrices –– band matrices with the bandwidth $m$ equal to 3 or 5 which called tridiagonal and pentadiagonal respectively.\n",
    "\n",
    "- (5 pts) Write a function ```band_lu(diag_broadcast, n)``` which computes LU decomposition for tridiagonal or pentadiagonal matrix with given diagonal values. \n",
    "For example, input parametres ```(diag_broadcast = [4,-2,1], n = 4)``` mean that we need to find LU decomposition for the triangular matrix of the form:\n",
    "\n",
    "$$A = \\begin{pmatrix}\n",
    "-2 & 1 & 0 & 0\\\\\n",
    "4 & -2 & 1 & 0 \\\\\n",
    "0 & 4 & -2 & 1 \\\\\n",
    "0 & 0 & 4 & -2 \\\\\n",
    "\\end{pmatrix}.$$\n",
    "\n",
    "As an output it is considered to make ```L``` and ```U``` - 2D arrays representing diagonals in factors $L$ (```L[0]``` keeps first lower diagonal, ```L[1]``` keeps second lower, ...), and $U$ (```U[:,0]``` keeps main diagonal, ```U[:,1]``` keeps first upper, ...).\n",
    "- (2 pts) Compare execution time of the band LU decomposition using standard function from ```scipy```, i.e. which takes the whole matrix and does not know about its special structure, and band decomposition of yours implementation. Comment on the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.sparse import diags # can be used with broadcasting of scalars if desired dimensions are large\n",
    "import numpy as np\n",
    "\n",
    "# INPUT : diag_broadcast - list of diagonals value to broadcast,length equal to 3 or 5; n - integer, band matrix shape.\n",
    "# OUTPUT : L - 2D np.ndarray, L.shape[0] depends on bandwidth, L.shape[1] = n-1, do not store main diagonal, where all ones;                  add zeros to the right side of rows to handle with changing length of diagonals.\n",
    "#          U - 2D np.ndarray, U.shape[0] = n, U.shape[1] depends on bandwidth;\n",
    "#              add zeros to the bottom of columns to handle with changing length of diagonals.\n",
    "def band_lu(diag_broadcast, n):\n",
    "    # enter your code here\n",
    "    raise NotImplementedError() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your solution is here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Stability of LU (8 pts)\n",
    "\n",
    "Let\n",
    "$A = \\begin{pmatrix}\n",
    "\\varepsilon & 1 & 0\\\\\n",
    "1 & 1 & 1 \\\\\n",
    "0 & 1 & 1\n",
    "\\end{pmatrix}.$ \n",
    "* (5 pts) Find analytically LU decomposition with and without pivoting for the matrix $A$.\n",
    "* (3 pts) Explain, why can the LU decomposition fail to approximate factors $L$ and $U$ for $|\\varepsilon|\\ll 1$ in computer arithmetic?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your solution is here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Block LU (10 pts)\n",
    "\n",
    "Let $A = \\begin{bmatrix} A_{11} & A_{12} \\\\ A_{21} & A_{22} \\end{bmatrix}$ be a block matrix. The goal is to solve the linear system\n",
    "\n",
    "$$\n",
    "     \\begin{bmatrix} A_{11} & A_{12} \\\\ A_{21} & A_{22} \\end{bmatrix} \\begin{bmatrix} u_1 \\\\ u_2 \\end{bmatrix} = \\begin{bmatrix} f_1 \\\\ f_2 \\end{bmatrix}.\n",
    "$$\n",
    "\n",
    "* (2 pts) Using block elimination find matrix $S$ and right-hand side $\\hat{f_2}$ so that $u_2$ can be found from $S u_2 = \\hat{f_2}$. Note that the matrix $S$ is called <font color='red'> Schur complement </font> of the block $A_{11}$.\n",
    "* (4 pts) Using Schur complement properties prove that \n",
    "\n",
    "$$\\det(X+AB) = \\det(X)\\det(I+BX^{-1}A), $$\n",
    "\n",
    "\n",
    "where $X$ - nonsingular square matrix.\n",
    "* (4 pts) Let matrix $F \\in \\mathbb{R}^{m \\times n}$ and $G \\in \\mathbb{R}^{n \\times m}$. Prove that \n",
    "\n",
    "$$\\det(I_m - FG) = \\det(I_n - GF).$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your solution is here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "## Problem 2 (QR decomposition) 30 pts\n",
    "\n",
    "### 1. Standard Gram-Schmidt algorithm (18 pts)\n",
    "Our goal is to orthogonalize a system of linearly independent vectors $v_1,\\dots,v_n$.\n",
    "The standard algorithm for this task is the Gram-Schmidt process:\n",
    "\n",
    "$$\n",
    "\\begin{split}\n",
    "u_1 &= v_1, \\\\\n",
    "u_2 &= v_2 - \\frac{(v_2, u_1)}{(u_1, u_1)} u_1, \\\\\n",
    "\\dots \\\\\n",
    "u_n &= v_n - \\frac{(v_n, u_1)}{(u_1, u_1)} u_1 - \\frac{(v_n, u_2)}{(u_2, u_2)} u_2 - \\dots - \\frac{(v_n, u_{n-1})}{(u_{n-1}, u_{n-1})} u_{n-1}.\n",
    "\\end{split}\n",
    "$$\n",
    "\n",
    "Obtained $u_1, \\dots, u_n$ are orthogonal vectors in exact arithmetics. Then to make the system orthonormal you should divide each of the vectors by its norm: $u_i := u_i/\\|u_i\\|$.\n",
    "The Gram-Schmidt process can be considered as a QR decomposition. Let us show that."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* (2 pts) Write out what is matrices $Q$ and $R$ obtained in the process above. \n",
    "\n",
    "* (5 pts) Implement the described Gram-Schmidt algorithm as a function ```gram_schmidt_qr(A)``` that takes a rectangular matrix ```A``` and outputs ```Q,R```.\n",
    "\n",
    "* (3 pts) Create a rectangular matrix $K = [x, Ax, A^2x, \\ldots, A^kx]$ where $A \\in \\mathbb{R}^{100 \\times 100}$ is a random matrix and $x$ is a random vector of the appropriate dimension, $k = 20$.  \n",
    "The loss of orthogonality can be described by the following error: $\\|Q^{\\top}Q-I\\|_2$, where $Q^{\\top}Q$ is called a Gram matrix. Compute QR decomposition of the created matrix $K$ with function that you have implemented and calculate error $\\|Q^{\\top}Q-I\\|_2$. Comment on the result.\n",
    "\n",
    "* (5 pts) The observed loss of orthogonality is a problem of this particular algorithm. Luckily, there is [a simple improvement to the algorithm above](https://en.wikipedia.org/wiki/Gram–Schmidt_process#Numerical_stability) that reduces the loss of orthogonality. Implement this modification as a function ```modified_gram_schmidt_qr(A)``` such that input and output are similar to ```gram_schmidt_qr(A)```. \n",
    "* (3 pts) Compute QR decomposition of the matrix $K$ from the previous task with the function ```modified_gram_schmidt_qr(A)```.\n",
    "Compute error $\\|Q^{\\top}Q-I\\|_2$. Compare this error to the error obtained with a \"pure\" Gram-Schmidt and comment on the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# INPUT : rectangular matrix A\n",
    "# OUTPUT: matrices Q - orthogonal and R - upper triangular such that A = QR\n",
    "def gram_schmidt_qr(A): # 5 pts\n",
    "    # your code is here\n",
    "    return Q, R\n",
    "\n",
    "# INPUT : rectangular matrix A\n",
    "# OUTPUT: matrices Q - orthogonal and R - upper triangular such that A = QR\n",
    "def modified_gram_schmidt_qr(A): # 5 pts\n",
    "    # your code is here\n",
    "    return Q, R"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Updating QR decompositions (12 pts)\n",
    "\n",
    "Suppose we are dealing with overdetermined system:\n",
    "\n",
    "$$ Ax = y, \\;\\; A \\in \\mathbb{R}^{m \\times n}, \\;\\;  m \\geq n. $$\n",
    "\n",
    "If you were attentive during the course, you should be aware of \n",
    "the canonical way to solve this system. \n",
    "\n",
    "- No, implementing **Householder $QR$** will not be your task :)\n",
    "\n",
    "- We assume that $QR$ decomposition of $A$ **is given**:\n",
    "\n",
    "$$ QR = A = \\begin{bmatrix}a_1 & a_2 & \\ldots & a_n \\end{bmatrix}, \\;\\; a_i \\in \\mathbb{R}^m,$$\n",
    "\n",
    "where $Q \\in \\mathbb{R}^{m \\times m}$ is unitary matrix and $R \\in \\mathbb{R}^{m \\times n}$ is upper triangular.\n",
    "\n",
    "- Let $k$ be an arbitrary column index, $1 \\leq k \\leq n$, and $z \\in \\mathbb{R}^{m}$ be an arbitrary vector.\n",
    "\n",
    "We are interested in matrices:\n",
    "\n",
    "$$ \\tilde{A} = \\begin{bmatrix}a_1 & \\ldots & a_{k-1} & a_{k+1} & \\ldots a_n \\end{bmatrix} \\in \\mathbb{R}^{m \\times (n-1)},$$\n",
    "\n",
    "and:\n",
    "\n",
    "$$ \\hat{A} = \\begin{bmatrix}a_1 & \\ldots & a_{k} & z & a_{k+1} & \\ldots a_n \\end{bmatrix} \\in \\mathbb{R}^{m \\times (n+1)},$$\n",
    "\n",
    "which are obtained by deleting / inserting one column from / to the matrix $A$.\n",
    "\n",
    "For example, if $A$ comes from a least-squares problem, such matrices appear if one wants to examine the significance of the $k^{th}$ parameter in the underlying model (e.g. think in terms of linear regression).\n",
    "\n",
    "- We want to compute decompostions of these updated matrices:\n",
    "\n",
    "$$ \\tilde{A} = \\tilde{Q}\\tilde{R}, \\;\\; \\hat{A} = \\hat{Q}\\hat{R}. $$\n",
    "\n",
    "Performing full orthogonalization process will be redundant. Instead, the known $QR$ factors of $A$ should be **updated** in a proper way.\n",
    "\n",
    "- Different methods for performing $QR$ decomposition were studied during the course. We expect the factors of updated matrices can be recovered with a small number of selective updates, hence the **Givens rotations** seem to be the most appropriate.\n",
    "\n",
    "- Recall the Givens rotation matrix :\n",
    "\n",
    "$$G_{ij}(\\theta) = \\begin{bmatrix} \n",
    "1      & \\ldots & 0      & \\ldots & 0      & \\ldots & 0      \\\\\n",
    "\\vdots & \\ddots & \\vdots &        & \\vdots &        & \\vdots \\\\\n",
    "0      & \\ldots & c      & \\ldots & s      & \\ldots & 0      \\\\\n",
    "\\vdots &        & \\vdots & \\ddots & \\vdots &        & \\vdots \\\\\n",
    "0      & \\ldots & -s     & \\ldots & c      & \\ldots & 0      \\\\\n",
    "\\vdots &        & \\vdots &        & \\vdots & \\ddots & \\vdots \\\\\n",
    "0      & \\ldots & 0      & \\ldots & 0      & \\ldots & 1      \n",
    "\\end{bmatrix}  \\;\\; \\in \\mathbb{R}^{m \\times m}, \\;\\; c = \\cos \\theta, \\;\\; s = \\sin \\theta,$$ \n",
    "\n",
    "which differs from the identity $\\mathrm{I}^{m \\times m}$ by a $2 \\times 2$ sub-matrix:\n",
    "\n",
    "$$\n",
    "{G}(\\theta) = \\begin{bmatrix} c & s \\\\ -s & c \\end{bmatrix} \\in \\mathbb{R}^{2 \\times 2},\n",
    "$$\n",
    "placed on the rows and columns with indices $i$ and $j$. \n",
    "\n",
    "1. (1 pts)\n",
    "Matrix ${G}$ has clear geometric interpretation in 2 dimensions. Let  $a = \\begin{pmatrix} a_1 \\\\ a_2 \\end{pmatrix} \\in \\mathbb{R}^2$ be an arbitrary non-zero vector. According to the template provided below, implement the function to construct matrix $G$, such that:\n",
    "\n",
    "$$ {G}^T \\begin{pmatrix} a_1 \\\\ a_2 \\end{pmatrix} = \\begin{pmatrix} \\alpha \\\\ 0 \\end{pmatrix}.$$\n",
    "\n",
    "What are geometrical meanings of $\\alpha$ and $\\theta$? \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def givens(a):\n",
    "    '''\n",
    "    Input: \n",
    "        a: np.array of size (2)\n",
    "    Output: \n",
    "        G: np.array of size (2, 2)\n",
    "            Rotation matrix, G^T eliminates the second component of a\n",
    "    '''\n",
    "\n",
    "    # Your code is here\n",
    "    \n",
    "    raise NotImplementedError()\n",
    "    \n",
    "# Check yourself \n",
    "# a = np.random.rand(2)\n",
    "# givens(a).T.dot(a) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Similarly, matrix $G_{ij}$ can be constructed in order to eliminate the $j^{th}$ component of an arbitrary vector of size $m$. This leads to the well-known **Givens $QR$ algorithm** (see lecture notes!).\n",
    "\n",
    "- Givens rotations typically appear in matrix form:\n",
    "\n",
    "$$ G_{ij}^T X, \\;\\; X^{'}G_{ij},$$\n",
    "\n",
    "where $X$ and $X^{'}$ are some matrices with compatible sizes, e.g. intermediate factors in the **Givens $QR$ algorithm**. \n",
    "\n",
    "How to perform these multiplications efficiently? \n",
    "Obviously, $G_{ij}$ does not need to be stored in explicit form. It affects only two rows / columns, hence only appropriate submatrices should be modified.\n",
    "\n",
    "2. (5 pts) Your second task is to propose an efficient procedure for computing $\\tilde{Q}, \\tilde{R}$ factors of $\\tilde{A}$ by a proper update of $Q, R$ factors of $A$. \n",
    "\n",
    "Hint: consider the following block partition of $R$ (assuming $k^{th}$ column deleted):\n",
    "$$ R = \n",
    "\\begin{bmatrix} \n",
    "R_{11} & v & R_{13} \\\\\n",
    "0 & r_{kk} & \\omega^T \\\\\n",
    "0 & 0 & R^{33}\n",
    "\\end{bmatrix},\n",
    "$$\n",
    "and note, that:\n",
    "$$Q^T\\tilde{A} = \n",
    "\\begin{bmatrix}\n",
    "R_{11} & R_{13} \\\\\n",
    "0 & \\omega^T \\\\\n",
    "0 & R_{33}\n",
    "\\end{bmatrix} = H.$$\n",
    "\n",
    "- (2 pts) What structure has matrix $H$ and how does it help in understanding the proper way to compute updated QR decomposition?\n",
    "\n",
    "- (3 pts) Implement the function according to the template provided below.\n",
    "Use Givens matrices $G_{ij}$ with $j = i+1$ for the convinience.\n",
    "\n",
    "- Only one loop over non-zero subdiagonal elements of $H$ is allowed. You are supposed to perform computations in a proper vectorized form. For example, the *in-place* multiplication $G_{ij}^T X$ with a random matrix $X \\in \\mathbb{R}^{m \\times n}$ can be implemented as follows: ```X[[i,j],:] = G.T.dot(X[[i,j],:])```. \n",
    "\n",
    "- Do not modify elements which are known to be zero.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def delete_column_QR(A, k):\n",
    "    '''\n",
    "    Input: \n",
    "        A: np.array of size (m, n) \n",
    "            \n",
    "        k: positive integer\n",
    "            Index of deleted column, 0 \\leq k \\leq n - 1\n",
    "        \n",
    "    Output: \n",
    "        Q: np.array of size (m, m)\n",
    "        \n",
    "        R: np.array of size (m, n - 1)\n",
    "\n",
    "    '''\n",
    "    Q, R = np.linalg.qr(A, mode='complete') # factors of $A$ are known\n",
    "    R = np.delete(R, k, axis = 1) # delete column\n",
    "        \n",
    "    # Your code is here\n",
    "        \n",
    "    raise NotImplementedError()\n",
    "\n",
    "# Check yourself\n",
    "m, n = (200, 100)\n",
    "A = np.random.rand(m, n)\n",
    "k = np.random.randint(0, n)\n",
    "A_tilde = np.delete(A, k, axis=1)\n",
    "Q_tilde, R_tilde = delete_column_QR(A, k)\n",
    "print(np.linalg.norm(A_tilde - Q_tilde.dot(R_tilde)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. (3 pts) When a column is inserted, the factors $\\hat{Q}, \\hat{R}$ of $\\hat{A}$ again can be computed efficiently by a proper update of $Q, R$ factors of $A$.\n",
    "\n",
    "\n",
    "- Proceed similarly to the previous task. Implement the function according to the template provided below. \n",
    "\n",
    "- Choose optimal order of elimination to minimize the number of rotations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def insert_column_QR(A, k, z):\n",
    "    '''\n",
    "    Input: \n",
    "        A: np.array of size (m, n) \n",
    "            \n",
    "        k: positive integer\n",
    "            Index of inserted column, 0 \\leq k leq n - 1\n",
    "            \n",
    "        z: np.array of size (m)\n",
    "            Inserted column\n",
    "        \n",
    "    Output: \n",
    "        Q: np.array of size (m, m)\n",
    "        \n",
    "        R: np.array of size (m, n + 1)\n",
    "    '''\n",
    "    Q, R = np.linalg.qr(A, mode='complete') # factors of $A$ are known\n",
    "\n",
    "    # Your code is here\n",
    "    \n",
    "    raise NotImplementedError()\n",
    "\n",
    "\n",
    "# Check yourself\n",
    "m, n = (200, 100)\n",
    "A = np.random.rand(m, n)\n",
    "k = np.random.randint(0, n)\n",
    "z = np.random.rand(m)\n",
    "A_hat = np.insert(A, k, z, axis=1)\n",
    "Q_hat, R_hat = insert_column_QR(A, k, z)\n",
    "print(np.linalg.norm(A_hat - Q_hat.dot(R_hat)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. (3 pts) In the context of least-squares problems, appending rows to a matrix is just as important as inserting columns (e.g. new samples in linear regression).\n",
    "\n",
    "Consider matrix:\n",
    "\n",
    "$$ \\overline{A} = \\begin{bmatrix} \\omega^T \\\\ A \\end{bmatrix} \\in \\mathbb{R}^{(m+1) \\times n},$$\n",
    "\n",
    "which is obtained by appending a row $\\omega^T \\in \\mathbb{R}^{n}$ to the matrix $A$.\n",
    "\n",
    "Similarly to the previous tasks, you need to compute $\\overline{Q}, \\overline{R}$ factors of $\\overline{A}$ by a proper update of $Q,R$ factors of $A$.\n",
    "\n",
    "- Implement the function according to the template provided below. Your code is still expected to be efficient!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def insert_row_QR(A, w):\n",
    "    '''\n",
    "    Input: \n",
    "        A: np.array of size (m, n)\n",
    "\n",
    "        w: np.array of size (n)\n",
    "            Inserted row\n",
    "        \n",
    "    Output:\n",
    "        Q: np.array of size (m + 1, m + 1)\n",
    "        \n",
    "        R: np.array of size (m + 1, n)\n",
    "    '''\n",
    "    Q, R = np.linalg.qr(A, mode='complete')\n",
    "\n",
    "    # Your code is here\n",
    "    \n",
    "    raise NotImplementedError()\n",
    "\n",
    "G\n",
    "# Check yourself\n",
    "m, n = (200, 100)\n",
    "A = np.random.rand(m, n)\n",
    "w = np.random.rand(n)\n",
    "A_overline = np.insert(A, 0, w, axis=0)\n",
    "Q_overline, R_overline = insert_row_QR(A, w)\n",
    "print(np.linalg.norm(A_overline - Q_overline.dot(R_overline)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 3 (Word2Vec as Matrix Factorization) 45 pts\n",
    "\n",
    "In this assignment you are supposed to apply SVD to training your own [word embedding model](https://en.wikipedia.org/wiki/Word_embedding) which maps English words to vectors of real numbers.\n",
    "\n",
    "Skip-Gram Negative Sampling (SGNS) word embedding model, commonly known as **word2vec** ([Mikolov et al., 2013](http://papers.nips.cc/paper/5021-distributed-representations-of-words-and-phrases-and-their-compositionality.pdf)), is usually optimized by stochastic gradient descent. However, the optimization of SGNS objective can be viewed as implicit matrix factorization objective as was shown in ([Levy and Goldberg, 2015](http://papers.nips.cc/paper/5477-neural-word-embedding-as-implicit-matrix-factorization.pdf))."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### 1. Notation\n",
    "Assume we have a text corpus given as a sequence of words $\\{w_1,w_2,\\dots,w_n\\}$ where $n$ may be larger than $10^{12}$ and $w_i \\in \\mathcal{V}$ belongs to a vocabulary of words $\\mathcal{V}$. A word $c \\in \\mathcal{V}$ is called *a context* of word $w_i$ if they are found together in the text. More formally, given some measure $L$ of closeness between two words (typical choice is $L=2$), a word $c \\in \\mathcal{V}$ is called a context if $c \\in \\{w_{i-L}, \\dots, w_{i-1}, w_{i+1}, \\dots, w_{i+L} \\}$ Let $\\mathbf{w},\\mathbf{c}\\in\\mathbb{R}^d$ be the *word embeddings* of word $w$ and context $c$, respectively. Assume they are specified by the mapping  $\\Phi:\\mathcal{V}\\rightarrow\\mathbb{R}^d$, so $\\mathbf{w}=\\Phi(w)$. The ultimate goal of SGNS word embedding model is to fit a good mapping $\\Phi$.\n",
    "\n",
    "Let $\\mathcal{D}$ be a multiset of all word-contexts pairs observed in the corpus. In the SGNS model, the probability that word-context pair $(w,c)$ is observed in the corpus is modeled as the following distribution:\n",
    "\n",
    "$$\n",
    "P(\\#(w,c)\\neq 0|w,c) = \\sigma(\\mathbf{w}^\\top \\mathbf{c}) = \\frac{1}{1 + \\exp(-\\mathbf{w}^\\top \\mathbf{c})},\n",
    "$$\n",
    "\n",
    "where $\\#(w,c)$ is the number of times the pair $(w,c)$ appears in $\\mathcal{D}$ and $\\mathbf{w}^\\top\\mathbf{c}$ is the scalar product of vectors $\\mathbf{w}$ and $\\mathbf{c}$. Two important quantities which we will also use further are the number of times the word $w$ and the context $c$ appear in $\\mathcal{D}$, which can be computed as\n",
    "\n",
    "$$\n",
    "\\#(w) = \\sum_{c\\in\\mathcal{V}} \\#(w,c), \\quad \\#(c) = \\sum_{w\\in\\mathcal{V}} \\#(w,c).\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### 2. Optimization objective"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vanilla word embedding models are trained by maximizing log-likelihood of observed word-context pairs, namely\n",
    "\n",
    "$$\n",
    "\\mathcal{L} = \\sum_{w \\in \\mathcal{V}} \\sum_{c \\in \\mathcal{V}} \\#(w,c) \\log \\sigma(\\mathbf{w}^\\top\\mathbf{c}) \\rightarrow \\max_{\\mathbf{w},\\mathbf{c} \\in \\mathbb{R}^d}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Skip-Gram Negative Sampling approach modifies the objective by additionally minimizing the log-likelihood of random word-context pairs, so called *negative samples*. This idea incorporates some useful linguistic information that some number ($k$, usually $k=5$) of word-context pairs *are not* found together in the corpus which usually results in word embeddings of higher quality. The resulting optimization problem is\n",
    "\n",
    "$$\n",
    "\\mathcal{L} = \\sum_{w \\in \\mathcal{V}} \\sum_{c \\in \\mathcal{V}} \\left( \\#(w,c) \\log \\sigma(\\mathbf{w}^\\top\\mathbf{c}) + k \\cdot \\mathbb{E}_{c'\\sim P_\\mathcal{D}} \\log \\sigma (-\\mathbf{w}^\\top\\mathbf{c}) \\right) \\rightarrow \\max_{\\mathbf{w},\\mathbf{c} \\in \\mathbb{R}^d},\n",
    "$$\n",
    "\n",
    "where $P_\\mathcal{D}(c)=\\frac{\\#(c)}{|\\mathcal{D}|}$ is a probability distribution over word contexts from which negative samples are drawn."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Levy and Goldberg, 2015](http://papers.nips.cc/paper/5477-neural-word-embedding-as-implicit-matrix-factorization.pdf) showed that this objective can be equivalently written as\n",
    "\n",
    "$$\n",
    "\\mathcal{L} = \\sum_{w \\in \\mathcal{V}} \\sum_{c \\in \\mathcal{V}} f(w,c) = \\sum_{w \\in \\mathcal{V}} \\sum_{c \\in \\mathcal{V}} \\left( \\#(w,c) \\log \\sigma(\\mathbf{w}^\\top\\mathbf{c}) + \\frac{k\\cdot\\#(w)\\cdot\\#(c)}{|\\mathcal{D}|} \\log \\sigma (-\\mathbf{w}^\\top\\mathbf{c}) \\right) \\rightarrow \\max_{\\mathbf{w},\\mathbf{c} \\in \\mathbb{R}^d},\n",
    "$$\n",
    "\n",
    "A crucial observation is that this loss function depends only on the scalar product $\\mathbf{w}^\\top\\mathbf{c}$ but not on embedding $\\mathbf{w}$ and $\\mathbf{c}$ separately."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### 3. Matrix factorization problem statement"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let $|\\mathcal{V}|=m$, $W \\in \\mathbb{R}^{m\\times d}$ and $C \\in \\mathbb{R}^{m\\times d}$ be matrices, where each row $\\mathbf{w}\\in\\mathbb{R}^d$ of matrix $W$ is the word embedding of the corresponding word $w$ and each row $\\mathbf{c}\\in\\mathbb{R}^d$ of matrix $C$ is the context embedding of the corresponding context $c$. SGNS embeds both words and their contexts into a low-dimensional space $\\mathbb{R}^d$, resulting in the word and context matrices $W$ and $C$. The rows of matrix $W$ are typically used in NLP tasks (such as computing word similarities) while $C$ is ignored. It is nonetheless instructive to consider the product $W^\\top C = M$. Viewed this way, SGNS can be described as factorizing an implicit matrix $M$ of dimensions $m \\times m$ into two smaller matrices.\n",
    "\n",
    "Which matrix is being factorized? A matrix entry $M_{wc}$ corresponds to the dot product $\\mathbf{w}^\\top\\mathbf{c}$ . Thus, SGNS is factorizing a matrix in which each row corresponds to a word $w \\in \\mathcal{V}$ , each column corresponds to a context $c \\in \\mathcal{V}$, and each cell contains a quantity $f(w,c)$ reflecting the strength of association between that particular word-context pair. Such word-context association matrices are very common in the NLP and word-similarity literature. That said, the objective of SGNS does not explicitly state what this association metric is. What can we say about the association function $f(w,c)$? In other words, which matrix is SGNS factorizing? Below you will find the answers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 1 (theoretical) 9 pts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve SGNS optimization problem with respect to the $\\mathbf{w}^\\top\\mathbf{c}$ and show that the matrix being factorized is\n",
    "\n",
    "$$\n",
    "M_{wc} = \\mathbf{w}^\\top\\mathbf{c} = \\log \\left( \\frac{\\#(w,c) \\cdot |\\mathcal{D}|}{k\\cdot\\#(w)\\cdot\\#(c)} \\right)\n",
    "$$\n",
    "\n",
    "**Hint:** Denote $x=\\mathbf{w}^\\top\\mathbf{c}$, rewrite SGNG optimization problem in terms of $x$ and solve it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note:** This matrix is called Shifted Pointwise Mutual Information (SPMI) matrix, as its elements can be written as\n",
    "\n",
    "$$\n",
    "\\text{SPMI}(w,c) = M_{wc} = \\mathbf{w}^\\top\\mathbf{c} = \\text{PMI}(w,c) - \\log k\n",
    "$$\n",
    "\n",
    "and $\\text{PMI}(w,c) = \\log \\left( \\frac{\\#(w,c) \\cdot |\\mathcal{D}|}{\\#(w)\\cdot\\#(c)} \\right)$ is the well-known [pointwise mutual information](https://en.wikipedia.org/wiki/Pointwise_mutual_information) of $(w,c)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your proof is here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 2 (practical) 36 pts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "from sklearn.metrics.pairwise import cosine_similarity\n",
    "from scipy.sparse.linalg import svds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Download dataset [enwik8](http://mattmahoney.net/dc/enwik8.zip) of compressed Wikipedia articles and preprocess raw data with Perl script **main_.pl**. This script will clean all unnecessary symbols, make all words to lowercase, and produce only sentences with words. If you are Windows user, you can install [linux subsystem in windows](https://www.microsoft.com/en-us/p/ubuntu-1804-lts/9n9tngvndl3q?activetab=pivot:overviewtab) and run all these bash commands from ubuntu console (it was checked on the Windows and it works). Also, you can just download file from this [link](https://nla.skoltech.ru/homeworks/files/enwik8.txt).\n",
    "```\n",
    "wget http://mattmahoney.net/dc/enwik8.zip\n",
    "unzip enwik8.zip\n",
    "mkdir data\n",
    "perl main_.pl enwik8 > data/enwik8.txt\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load enwik 8\n",
    "\n",
    "import re\n",
    "file = open(\"data/enwik8.txt\", \"r\")\n",
    "doclist = [line for line in file]\n",
    "docstr = ''.join(doclist)\n",
    "sentences = re.split(r'[.!?]', docstr)\n",
    "sentences = [sentence.split() for sentence in sentences if len(sentence) > 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['achilles', 'wrath', 'is', 'terrible', 'and', 'he', 'slays', 'many', 'trojan', 'warriors', 'and', 'allies', 'including', 'priam', 's', 'son', 'lycaon', 'whom', 'achilles', 'had', 'previously', 'captured', 'and', 'sold', 'into', 'slavery', 'but', 'who', 'had', 'been', 'returned', 'to', 'troy']\n"
     ]
    }
   ],
   "source": [
    "print (sentences[1249])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Construct the word vocabulary from the obtained sentences which enumerates words which occur more than $r=200$ times in the corpus."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_vocabulary(sentences, r=200):\n",
    "    vocabulary = {}\n",
    "    # Your code is here\n",
    "    return vocabulary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "vocab = create_vocabulary(sentences)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. Scan the text corpus with sliding window of size $5$ and step $1$ (which corresponds to $L$=2) and construct co-occurrence word-context matrix $D$ with elements $D_{wc}=\\#(w,c)$. Please, ignore words which occur less than $r=200$ times, but include them into the sliding window. Please, see the graphical illustration of the procedure described."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Sliding window](sliding_window.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_corpus_matrix(sentences, vocabulary):\n",
    "    # Your code is here\n",
    "    return corpus_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "D = create_corpus_matrix(sentences, vocab)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. To find good word embeddings, [Levy and Goldberg, 2015](http://papers.nips.cc/paper/5477-neural-word-embedding-as-implicit-matrix-factorization.pdf) proposed to find rank-$d$ SVD of Shifted Positive Pointwise Mutual Information (SPPMI) matrix\n",
    "\n",
    "$$ U \\Sigma V^\\top \\approx \\text{SPPMI}, $$\n",
    "\n",
    "where $\\text{SPPMI}(w, c) = \\max\\left(\\text{SPMI}(w, c), 0 \\right)$ and $\\text{SPMI}(w, c)$ is the element of the matrix $\\text{SPPMI}$ at position $(w, c)$.\n",
    "Then use $W=U\\sqrt{\\Sigma}$ as word embedding matrix. Your task is to reproduce their results. Write function constructs $\\text{SPPMI}$ matrix, computes its SVD and produces word-vectors matrix $W$. Pay attention that $\\text{SPPMI}$ matrix is **sparse**!  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_embeddings(D, k, d=200):\n",
    "    # Your code is here\n",
    "    return embedding_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 5 # negative sampling parameter\n",
    "W = compute_embeddings(D, k)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "5. Write class **WordVectors** using provided template."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "class WordVectors:\n",
    "    \n",
    "    def __init__(self, vocabulary, embedding_matrix):\n",
    "        self.vocab = vocabulary\n",
    "        self.W = embedding_matrix\n",
    "        self.inv_vocab = {v: k for k, v in self.vocab.items()}\n",
    "        \n",
    "    def word_vector(self, word):\n",
    "        \"\"\" \n",
    "        Takes word and returns its word vector.\n",
    "        \"\"\"\n",
    "#         Your code is here\n",
    "        return word_vector\n",
    "    \n",
    "    def nearest_words(self, word, top_n=10):\n",
    "        \"\"\" \n",
    "        Takes word from the vocabulary and returns its top_n\n",
    "        nearest neighbors in terms of cosine similarity.\n",
    "        \"\"\"\n",
    "#         Your code is here\n",
    "        return neighbors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = WordVectors(vocab, W)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('communism', 0.792),\n",
       " ('anarcho', 0.787),\n",
       " ('capitalism', 0.784),\n",
       " ('socialism', 0.752),\n",
       " ('liberalism', 0.727),\n",
       " ('criticisms', 0.705),\n",
       " ('capitalist', 0.662),\n",
       " ('fascism', 0.565),\n",
       " ('anarchist', 0.527),\n",
       " ('marxist', 0.518)]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.nearest_words(\"anarchism\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('ukraine', 0.671),\n",
       " ('russia', 0.629),\n",
       " ('poland', 0.55),\n",
       " ('belarus', 0.538),\n",
       " ('yugoslavia', 0.538),\n",
       " ('romania', 0.517),\n",
       " ('serbia', 0.507),\n",
       " ('austria', 0.5),\n",
       " ('hungary', 0.466),\n",
       " ('bulgaria', 0.43)]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.nearest_words(\"ussr\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('hop', 0.828),\n",
       " ('hip', 0.817),\n",
       " ('funk', 0.758),\n",
       " ('rock', 0.737),\n",
       " ('punk', 0.706),\n",
       " ('music', 0.676),\n",
       " ('pop', 0.666),\n",
       " ('scene', 0.659),\n",
       " ('band', 0.657),\n",
       " ('jazz', 0.625)]"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.nearest_words(\"rap\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "6. Calculate top 10 nearest neighbours with the corresponding cosine similarities for the words {matrix, multiplication, algorithm} and print them below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your solution is here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem 4 (eigenvalues)  45 pts\n",
    "\n",
    "### 1. Theoretical tasks (10 pts)\n",
    "\n",
    "* (5 pts) Prove that normal matrix is Hermitian iff its eigenvalues are real. Prove that normal matrix is unitary iff its eigenvalues satisfy $|\\lambda| = 1$. \n",
    "\n",
    "* (5 pts) The following problem illustrates instability of the Jordan form. Find theoretically the eigenvalues of the perturbed Jordan block:\n",
    "\n",
    "$$\n",
    "    J(\\varepsilon) = \n",
    "    \\begin{bmatrix} \n",
    "     \\lambda & 1 & & & 0 \\\\ \n",
    "     & \\lambda & 1 & & \\\\ \n",
    "     &  & \\ddots & \\ddots & \\\\ \n",
    "     & & & \\lambda & 1 \\\\ \n",
    "     \\varepsilon & & & & \\lambda  \\\\ \n",
    "    \\end{bmatrix}_{n\\times n}\n",
    "$$\n",
    "\n",
    "Comment how eigenvalues of $J(0)$ are perturbed for large $n$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your solution is here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. PageRank\n",
    "\n",
    "\n",
    "#### Damping factor importance\n",
    "\n",
    "* (5 pts) Write the function ```pagerank_matrix(G)``` that takes an adjacency matrix $G$ (in both sparse and dense formats) as an input and outputs the corresponding PageRank matrix $A$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# INPUT:  G - np.ndarray or sparse matrix\n",
    "# OUTPUT: A - np.ndarray (of size G.shape) or sparse matrix\n",
    "def pagerank_matrix(G):\n",
    "    # enter your code here\n",
    "    return A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* (3 pts) Find PageRank matrix $A$ that corresponds to the following graph: <img src=\"graph.png\" width='250'>\n",
    "What is its largest eigenvalue? What multiplicity does it have?\n",
    "\n",
    "\n",
    "* (5 pts) Implement the power method for a given matrix $A$, an initial guess $x_0$ and a number of iterations ```num_iter```. It should be organized as a function ```power_method(A, x0, num_iter)``` that outputs approximation to eigenvector $x$, eigenvalue $\\lambda$ and history of residuals $\\{\\|Ax_k - \\lambda_k x_k\\|_2\\}$. Make sure that the method converges to the correct solution on a matrix $\\begin{bmatrix} 2 & -1 \\\\ -1 & 2 \\end{bmatrix}$ which is known to have the largest eigenvalue equal to $3$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# INPUT:  A - np.ndarray (2D), x0 - np.ndarray (1D), num_iter - integer (positive)\n",
    "# OUTPUT: x - np.ndarray (of size x0), l - float, res - np.ndarray (of size num_iter + 1 [include initial guess])\n",
    "def power_method(A, x0, num_iter): # 5 pts\n",
    "    # enter your code here\n",
    "    return x, l, res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* (2 pts) Run the power method for the graph presented above and plot residuals $\\|Ax_k - \\lambda_k x_k\\|_2$ as a function of $k$ for ```num_iter=100``` and random initial guess ```x0```.  Explain the absence of convergence. \n",
    "\n",
    "\n",
    "* (2 pts) Consider the same graph, but with additional self loop at node 4 (self loop is an edge that connects a vertex with itself). Plot residuals as in the previous task and discuss the convergence. Now, run the power method with ```num_iter=100``` for 10 different initial guesses and print/plot the resulting approximated eigenvectors. Why do they depend on the initial guess?\n",
    "\n",
    "\n",
    "In order to avoid this problem Larry Page and Sergey Brin [proposed](http://ilpubs.stanford.edu:8090/422/1/1999-66.pdf) to use the following regularization technique:\n",
    "\n",
    "$$\n",
    "A_d = dA + \\frac{1-d}{N} \\begin{pmatrix} 1 & \\dots & 1 \\\\ \\vdots & & \\vdots \\\\ 1 & \\dots & 1 \\end{pmatrix},\n",
    "$$\n",
    "\n",
    "where $d$ is a small parameter in $[0,1]$ (typically $d=0.85$), which is called **damping factor**, $A$ is of size $N\\times N$. Now $A_d$ is the matrix with multiplicity of the largest eigenvalue equal to 1. \n",
    "Recall that computing the eigenvector of the PageRank matrix, which corresponds to the largest eigenvalue, has the following interpretation. Consider a person who stays in a random node of a graph (i.e. opens a random web page); at each step s/he follows one of the outcoming edges uniformly at random (i.e. opens one of the links). So the person randomly walks through the graph and the eigenvector we are looking for is exactly his/her stationary distribution â€” for each node it tells you the probability of visiting this particular node. Therefore, if the person has started from a part of the graph which is not connected with the other part, he will never get there.  In the regularized model, the person at each step follows one of the outcoming links with probability $d$ OR teleports to a random node from the whole graph with probability $(1-d)$.\n",
    "\n",
    "* (2 pts) Now, run the power method with $A_d$ and plot residuals $\\|A_d x_k - \\lambda_k x_k\\|_2$ as a function of $k$ for $d=0.97$, ```num_iter=100``` and a random initial guess ```x0```.\n",
    "\n",
    "* (5 pts) Find the second largest in the absolute value eigenvalue of the obtained matrix $A_d$. How and why is it connected to the damping factor $d$? What is the convergence rate of the PageRank algorithm when using damping factor?\n",
    "\n",
    "Usually, graphs that arise in various areas are sparse (social, web, road networks, etc.) and, thus, computation of a matrix-vector product for corresponding PageRank matrix $A$ is much cheaper than $\\mathcal{O}(N^2)$. However, if $A_d$ is calculated directly, it becomes dense and, therefore, $\\mathcal{O}(N^2)$ cost grows prohibitively large for  big $N$.\n",
    "\n",
    "\n",
    "* (2 pts) Implement fast matrix-vector product for $A_d$ as a function ```pagerank_matvec(A, d, x)```, which takes a PageRank matrix $A$ (in sparse format, e.g., ```csr_matrix```), damping factor $d$ and a vector $x$ as an input and returns $A_dx$ as an output. \n",
    "\n",
    "* (1 pts) Generate a random adjacency matrix of size $10000 \\times 10000$ with only 100 non-zero elements and compare ```pagerank_matvec``` performance with direct evaluation of $A_dx$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# INPUT:  A - np.ndarray (2D), d - float (from 0.0 to 1.0), x - np.ndarray (1D, size of A.shape[0/1])\n",
    "# OUTPUT: y - np.ndarray (1D, size of x)\n",
    "def pagerank_matvec(A, d, x): # 2 pts\n",
    "    # enter your code here\n",
    "    return y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### DBLP: computer science bibliography\n",
    "\n",
    "Download the dataset from [here](https://goo.gl/oZVxEa), unzip it and put `dblp_authors.npz`  and `dblp_graph.npz` in the same folder with this notebook. Each value (author name) from `dblp_authors.npz` corresponds to the row/column of the matrix from `dblp_graph.npz`. Value at row `i` and column `j` of the matrix from `dblp_graph.npz` corresponds to the number of times author `i` cited papers of the author `j`. Let us now find the most significant scientists according to PageRank model over DBLP data.\n",
    "\n",
    "* (4 pts) Load the weighted adjacency matrix and the authors list into Python using ```load_dblp(...)``` function. Print its density (fraction of nonzero elements). Find top-10 most cited authors from the weighted adjacency matrix. Now, make all the weights of the adjacency matrix equal to 1 for simplicity (consider only existence of connection between authors, not its weight). Obtain the PageRank matrix $A$ from the adjacency matrix and verify that it is stochastic.\n",
    " \n",
    " \n",
    "* (1 pts) In order to provide ```pagerank_matvec``` to your ```power_method``` (without rewriting it) for fast calculation of $A_dx$, you can create a ```LinearOperator```: \n",
    "```python\n",
    "L = scipy.sparse.linalg.LinearOperator(A.shape, matvec=lambda x, A=A, d=d: pagerank_matvec(A, d, x))\n",
    "```\n",
    "Calling ```L@x``` or ```L.dot(x)``` will result in calculation of ```pagerank_matvec(A, d, x)``` and, thus, you can plug $L$ instead of the matrix $A$ in the ```power_method``` directly. **Note:** though in the previous subtask graph was very small (so you could disparage fast matvec implementation), here it is very large (but sparse), so that direct evaluation of $A_dx$ will require $\\sim 10^{12}$ matrix elements to store - good luck with that (^_<).\n",
    "\n",
    "\n",
    "* (2 pts) Run the power method starting from the vector of all ones and plot residuals $\\|A_dx_k - \\lambda_k x_k\\|_2$  as a function of $k$ for $d=0.85$.\n",
    "\n",
    "\n",
    "* (1 pts) Print names of the top-10 authors according to PageRank over DBLP when $d=0.85$. Comment on your findings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy.sparse import load_npz\n",
    "import numpy as np\n",
    "def load_dblp(path_auth, path_graph):\n",
    "    G = load_npz(path_graph).astype(float)\n",
    "    with np.load(path_auth) as data: authors = data['authors']\n",
    "    return G, authors\n",
    "G, authors = load_dblp('dblp_authors.npz', 'dblp_graph.npz')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code is here"
   ]
  }
 ],
 "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.9"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}