lectures/lecture-20/lecture-20.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Lecture 20:  Tensors and tensor decompositions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Recap of the previous lecture\n",
    "- Matrix functions\n",
    "- Matrix exponential\n",
    "- Pade approximation and rational Krylov subspaces\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Today lecture \n",
    "\n",
    "- Tensor decompositions and applications\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Tensors\n",
    "\n",
    "By **tensor** we imply a **multidimensional array**:\n",
    "\n",
    "$$\n",
    "A(i_1, \\dots, i_d), \\quad 1\\leq i_k\\leq n_k,\n",
    "$$\n",
    "\n",
    "where $d$ is called dimension, $n_k$ - mode size.\n",
    "This is standard definition in applied mathematics community. For details see [[1]](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.153.2059&rep=rep1&type=pdf), [[2]](http://arxiv.org/pdf/1302.7121.pdf), [[3]](http://epubs.siam.org/doi/abs/10.1137/090752286).\n",
    "\n",
    "* $d=2$ (matrices) $\\Rightarrow$ classic theory (SVD, LU, QR, $\\dots$)\n",
    "\n",
    "* $d\\geq 3$ (tensors) $\\Rightarrow$ under development. Generalization of standard matrix results is not  **straightforward**.\n",
    "<img src=\"./tensor_dogs.png\">\n",
    "\n",
    "Picture is taken from [this presentation](http://users.cecs.anu.edu.au/~koniusz/tensors-cvpr17/present/anandkumar_anima_tmcv2017.pdf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## More formal definition\n",
    "\n",
    "Tensor is a multilinear form. When you fix the basis, you get a $d$-dimensional table."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Curse of dimensionality\n",
    "\n",
    "The problem with multidimensional data is that number of parameters grows <font color='red'> exponentially </font> with $d$:\n",
    "\n",
    "\n",
    "$$\n",
    "    \\text{storage} = n^d.\n",
    "$$\n",
    "For instance, for $n=2$ and $d=500$\n",
    "$$\n",
    "    n^d = 2^{500} \\gg 10^{83}  -  \\text{ number of atoms in the Universe}\n",
    "$$\n",
    "\n",
    "Why do we care? It seems that we are living in the 3D World :)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Applications\n",
    "\n",
    "#### Quantum chemistry\n",
    "\n",
    "Stationary Schroedinger equation for system with $N_{el}$ electrons\n",
    "\n",
    "$$\n",
    "    \\hat H \\Psi = E \\Psi,\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "\\Psi = \\Psi(\\{{\\bf r_1},\\sigma_1\\},\\dots, \\{{\\bf r_{N_{el}}},\\sigma_{N_{el}}\\})\n",
    "$$\n",
    "\n",
    "3$N_{el}$ spatial variables and $N_{el}$ spin variables. \n",
    "<img src=\"./large_mol.jpg\" width=600>\n",
    "\n",
    "* Drug and material design\n",
    "* Predicting physical experiments"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Uncertainty quantification\n",
    "\n",
    "Example: oil reservoir modeling. Model may depend on parameters $p_1,\\dots,p_d$ (like measured experimentally porocity or temperature) with uncertainty\n",
    "\n",
    "$$\n",
    "u = u(t,{\\bf r},\\,{p_1,\\dots,p_d})\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### And many more\n",
    "\n",
    "* Signal processing\n",
    "* Recommender systems\n",
    "* Neural networks\n",
    "* Language models\n",
    "* Financial mathematics\n",
    "* ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Working with many dimensions\n",
    "\n",
    "How to work with high-dimensional functions?\n",
    "\n",
    "- **Monte-Carlo**: class of methods based on random sampling. Convergence issues\n",
    "- **Sparse grids**: special types of grids with small number of grid points. Strong regularity conditions\n",
    "- **Best N-term approximation** : sparse expansions in certain basis.\n",
    "- Promising approach based on <font color='red'>tensor decompositions </font>\n",
    "- Neural networks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Tensor decompositions\n",
    "\n",
    "## 2D\n",
    "\n",
    "Skeleton decomposition:\n",
    "$$\n",
    "A = UV^T\n",
    "$$\n",
    "or elementwise:\n",
    "$$\n",
    "a_{ij} = \\sum_{\\alpha=1}^r u_{i\\alpha} v_{j\\alpha}\n",
    "$$\n",
    "leads us to the idea of **separation of variables.**\n",
    "\n",
    "**Properties:**\n",
    "* Not unique: $A = U V^T = UBB^{-1}V^T = \\tilde U \\tilde V^T$\n",
    "* Can be calculated in a stable way by **SVD**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Canonical decomposition\n",
    "\n",
    "The most straightforward way to generize separation of variables to many dimensions is the **canonical decomposition**: (CP/CANDECOMP/PARAFAC)\n",
    "\n",
    "$$\n",
    "a_{ijk} = \\sum_{\\alpha=1}^r u_{i\\alpha} v_{j\\alpha} w_{k\\alpha}\n",
    "$$\n",
    "\n",
    "minimal possible $r$ is called the **canonical rank**. Matrices $U$, $V$ and $W$ are called **canonical factors**. This decomposition was proposed in 1927 by Hitchcock, [link](https://onlinelibrary.wiley.com/doi/abs/10.1002/sapm192761164)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**Properties**:\n",
    "\n",
    "* For a $d$ dimensional tensor memory is $nrd$\n",
    "* Unique under mild conditions\n",
    "* Set of tensors with rank$\\leq r$ is not closed (by contrast to matrices): <br>\n",
    "  $a_{ijk} = i+j+k$, $\\text{rank}(A) = 3$, but\n",
    "  \n",
    "  $$a^\\epsilon_{ijk} = \\frac{(1+\\epsilon i)(1+\\epsilon j)(1+\\epsilon k) - 1}{\\epsilon}\\to i+j+k=a_{ijk},\\quad \\epsilon\\to 0 $$\n",
    "  \n",
    "  and $\\text{rank}(A^{\\epsilon}) = 2$\n",
    "* No stable algorithms to compute best rank-$r$ approximation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Alternating Least Squares algorithm\n",
    "\n",
    "0. Intialize random $U,V,W$\n",
    "1. fix $V,W$, solve least squares for $U$\n",
    "2. fix $U,W$, solve least squares for $V$\n",
    "3. fix $U,V$, solve least squares for $W$\n",
    "4. go to 2."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### DNN compression ([Lebedev, et. al 2015](https://arxiv.org/pdf/1412.6553.pdf))\n",
    "\n",
    "- Convolution is represented as 4D tensor\n",
    "- Apply CP decomposition to this tensor\n",
    "- Apply one-by-one convolutions with smaller kernels given by CP decomposition\n",
    "- Fine tune such model\n",
    "\n",
    "<img src=\"./cp_compress_dnn.png\" width=600>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Tucker decomposition\n",
    "\n",
    "Next attempt is the decomposition proposed by ([Tucker, 1963](https://link.springer.com/content/pdf/10.1007/BF02289464.pdf)) in [Psychometrika](https://link.springer.com/journal/11336):\n",
    "\n",
    "$$\n",
    "a_{ijk} = \\sum_{\\alpha_1,\\alpha_2,\\alpha_3=1}^{r_1,r_2,r_3}g_{\\alpha_1\\alpha_2\\alpha_3} u_{i\\alpha_1} v_{j\\alpha_2} w_{k\\alpha_3}.\n",
    "$$\n",
    "\n",
    "Here we have several different ranks. Minimal possible $r_1,r_2,r_3$ are called **Tucker ranks**.\n",
    "\n",
    "**Properties**:\n",
    "\n",
    "* For a $d$ dimensional tensor memory is <font color='red'>$r^d$</font> $+ nrd$. Still **curse of dimensionality**\n",
    "* Stable SVD-based algorithm:\n",
    "    1. $U =$ principal components of the unfolding `A.reshape(n1, n2*n3)`\n",
    "    2. $V =$ principal components of the unfolding `A.transpose([1,0,2]).reshape(n2, n1*n3)`\n",
    "    3. $W =$ principal components of the unfolding `A.transpose([2,0,1]).reshape(n3, n1*n2)`\n",
    "    4. $g_{\\alpha_1\\alpha_2\\alpha_3} = \\sum_{i,j,k=1}^{n_1,n_2,n_3} a_{ijk} u_{i\\alpha_1} v_{j\\alpha_2} w_{k\\alpha_3}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shape of core = (3, 4, 2)\n",
      "Shape of factors = (3, 3)\n",
      "Shape of factors = (4, 4)\n",
      "Shape of factors = (2, 2)\n",
      "Approximation error = 3.181718859384126e-14\n"
     ]
    }
   ],
   "source": [
    "tensor = tl.tensor(np.arange(24).reshape((3, 4, 2)))\n",
    "core, factors = tldec.tucker(tensor, ranks=[3, 4, 2])\n",
    "print(\"Shape of core = {}\".format(core.shape))\n",
    "for f in factors:\n",
    "    print(\"Shape of factors = {}\".format(f.shape))\n",
    "print(\"Approximation error = {}\".format(np.linalg.norm(tensor - tl.tucker_to_tensor(core, factors))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Application in recommender systems [(Frolov, Oseledets 2016)](https://dl.acm.org/citation.cfm?id=2959170)\n",
    "\n",
    "- Represent User-Item **matrix** as binary User-Item-Rating **tensor**\n",
    "- After Tucker decomposition of this tensor we have three factors\n",
    "- Two of them represent projections on the space of ratings and space of items\n",
    "- This approach takes into account negative feedback to items from some users"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## CP and Tucker decompositions implementations\n",
    "\n",
    "- Matlab: [Tensorlab](https://www.tensorlab.net/) and [Tensor Toolbox](https://www.sandia.gov/~tgkolda/TensorToolbox/index-2.5.html)\n",
    "- Python: [TensorLy](http://tensorly.org/stable/home.html) (I don't like it) and [Scikit-tensor](https://github.com/mnick/scikit-tensor)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Tensor Train decomposition\n",
    "\n",
    "* Calculation of the canonical decomposition is unstable\n",
    "* Tucker decomposition suffers from the curse of dimensionality\n",
    "\n",
    "Tensor Train (**TT**) decomposition ([Oseledets, Tyrtyshnikov 2009](http://www.mat.uniroma2.it/~tvmsscho/papers/Tyrtyshnikov4.pdf) and [Oseledets, 2011](https://epubs.siam.org/doi/abs/10.1137/090752286)) is both stable and contains linear in $d$ number of parameters:\n",
    "\n",
    "$$\n",
    "a_{i_1 i_2 \\dots i_d} = \\sum_{\\alpha_1,\\dots,\\alpha_{d-1}} \n",
    "g_{i_1\\alpha_1} g_{\\alpha_1 i_2\\alpha_2}\\dots g_{\\alpha_{d-2} i_{d-1}\\alpha_{d-1}} g_{\\alpha_{d-1} i_{d}}\n",
    "$$\n",
    "\n",
    "or in the matrix form\n",
    "\n",
    "$$\n",
    "    a_{i_1 i_2 \\dots i_d} = G_1 (i_1)G_2 (i_2)\\dots G_d(i_d)\n",
    "$$\n",
    "\n",
    "* The storage is $\\mathcal{O}(dnr^2)$ \n",
    "* Stable TT-SVD algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**Example**\n",
    "$$a_{i_1\\dots i_d} = i_1 + \\dots + i_d$$\n",
    "Canonical rank is $d$. At the same time TT-ranks are $2$:\n",
    "$$\n",
    "i_1 + \\dots + i_d = \\begin{pmatrix} i_1 & 1 \\end{pmatrix} \n",
    "\\begin{pmatrix} 1 & 0 \\\\ i_2 & 1 \\end{pmatrix}\n",
    "\\dots\n",
    "\\begin{pmatrix} 1 & 0 \\\\ i_{d-1} & 1 \\end{pmatrix}\n",
    "\\begin{pmatrix} 1  \\\\  i_d \\end{pmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Implementations\n",
    "\n",
    "- Matlab: [TT Toolbox](https://github.com/oseledets/TT-Toolbox)\n",
    "- Python: [ttpy](https://github.com/oseledets/ttpy)\n",
    "- TensorFlow: [t3f](https://github.com/Bihaqo/t3f)\n",
    "- PyTorch: [tt-pytorch](https://github.com/KhrulkovV/tt-pytorch) and [tntorch](https://github.com/rballester/tntorch)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Using TT in Riemannien optimization (example is from t3f examples)\n",
    "\n",
    "- Given some TT tensor $A$ with large tt-ranks\n",
    "- We would like to find a tensor $X$ (with small prescribed tt-ranks $r$) which is closest to $A$ (in the sense of Frobenius norm):\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{aligned}\n",
    "& \\underset{X}{\\text{minimize}} \n",
    "& & \\frac{1}{2}\\|X - A\\|_F^2 \\\\\n",
    "& \\text{subject to} \n",
    "& & \\text{tt_rank}(X) = r\n",
    "\\end{aligned}\n",
    "\\end{equation*}\n",
    "\n",
    "- It is known that the set of TT tensors with elementwise fixed TT ranks forms a manifold. \n",
    "- Thus we can solve this problem using the so called Riemannian gradient descent. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Riemannian gradient descent\n",
    "\n",
    "- Given some functional $F$ on a manifold $\\mathcal{M}$ Riemannian gradient descent is defined as\n",
    "\n",
    "$$\\hat{x}_{k+1} = x_{k} - \\alpha P_{T_{x_k}\\mathcal{M}} \\nabla F(x_k),$$\n",
    "\n",
    "$$x_{k+1} = \\mathcal{R}(\\hat{x}_{k+1})$$\n",
    "\n",
    "with $P_{T_{x_k}\\mathcal{M}}$ being the projection onto the **tangent space** of $\\mathcal{M}$ at the point $x_k$ and $\\mathcal{R}$ being a **retraction** - an operation which projects points to the manifold, and $\\alpha$ is the learning rate.\n",
    "\n",
    "- We can implement this in `t3f` using the `t3f.riemannian` module. As a retraction it is convenient to use the rounding method (`t3f.round`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import t3f\n",
    "import tensorflow as tf\n",
    "tf.set_random_seed(0)\n",
    "np.random.seed(0)\n",
    "sess = tf.InteractiveSession()\n",
    "\n",
    "# Initialize A randomly, with large tt-ranks\n",
    "shape = 10 * [2]\n",
    "init_A = t3f.random_tensor(shape, tt_rank=16)\n",
    "A = t3f.get_variable('A', initializer=init_A, trainable=False)\n",
    "\n",
    "# Create an X variable and compute the gradient of the functional. Note that it is simply X - A.\n",
    "\n",
    "init_X = t3f.random_tensor(shape, tt_rank=2)\n",
    "X = t3f.get_variable('X', initializer=init_X)\n",
    "\n",
    "gradF = X - A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Let us compute the projection of the gradient onto the tangent space at X\n",
    "riemannian_grad = t3f.riemannian.project(gradF, X)\n",
    "\n",
    "# Compute the update by subtracting the Riemannian gradient\n",
    "# and retracting back to the manifold\n",
    "alpha = 1.0\n",
    "\n",
    "train_step = t3f.assign(X, t3f.round(X - alpha * riemannian_grad, max_tt_rank=2))\n",
    "\n",
    "# let us also compute the value of the functional\n",
    "# to see if it is decreasing\n",
    "F = 0.5 * t3f.frobenius_norm_squared(X - A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "81.62205\n",
      "58.53462\n",
      "56.27\n",
      "56.08319\n",
      "51.73249\n",
      "50.77673\n",
      "50.77672\n",
      "50.77671\n",
      "50.77672\n",
      "50.77672\n"
     ]
    }
   ],
   "source": [
    "sess.run(tf.global_variables_initializer())\n",
    "\n",
    "log = []\n",
    "for i in range(100):\n",
    "    F_v, _ = sess.run([F, train_step.op])\n",
    "    if i % 10 == 0:\n",
    "        print (F_v)\n",
    "    log.append(F_v)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "It is intructive to compare the obtained result with the quasioptimum delivered by the TT-round procedure. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "52.40742\n"
     ]
    }
   ],
   "source": [
    "quasi_sol = t3f.round(A, max_tt_rank=2)\n",
    "\n",
    "val = sess.run(0.5 * t3f.frobenius_norm_squared(quasi_sol - A))\n",
    "print (val)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We see that the value is slightly bigger than the exact minimum, but TT-round is faster and cheaper to compute, so it is often used in practice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0xd261fa4a8>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZkAAAEKCAYAAADAVygjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8VOXZ//HPlT2EEEjYA0pQFlkTJKCCSute0Fp3qrVuULUu1bZP7VOtS+vz68+f2qq1ta7Ypxaq4galalXQSmtZFAREcQEhgIBBwhqyXb8/ZjINIcskzGTIzPf9eo3hnDnnPteZE+fKfd/n3Le5OyIiItGQFOsAREQkfinJiIhI1CjJiIhI1CjJiIhI1CjJiIhI1CjJiIhI1CjJiIhI1CjJiIhI1CjJiIhI1KTEOoBY69q1q/fr1y/WYYiItCuLFy/+0t27NbddwieZfv36sWjRoliHISLSrpjZ5+Fsp+YyERGJGiUZERGJGiUZERGJmoTvkxGJhMrKSkpKSigvL491KCIRlZGRQZ8+fUhNTW3V/koyIhFQUlJCdnY2/fr1w8xiHY5IRLg7paWllJSUUFBQ0Koy1FwmEgHl5eXk5eUpwUhcMTPy8vIOqIauJCMSIUowEo8O9PdaSaaV/vr+RqbNXx3rMEREDmpKMq302spN/OGtz2IdhkhIcnIyhYWFDBs2jNNPP51t27YBsGHDBs4555wYR7e/n//857z22muxDqNBHTt2BA78s/vNb37D7t27m91u3rx5TJo0qdXHOVDbtm3jd7/7XVTKVpJppUE9s9lYVk7Z7spYhyICQGZmJkuWLGH58uXk5uby4IMPAtC7d2+effbZGEe3vzvuuIMTTzyxzY5XXV3d4n0O9LMLN8nEmpLMQWhQj2wAPtq0I8aRiOzv6KOPZv369QCsWbOGYcOGAYEv2h//+McUFxczYsQI/vCHPwCBv6SPP/54zjvvPAYOHMhNN93EU089xZgxYxg+fDiffvopALNmzWLs2LEUFRVx4oknsmnTJgBuu+02LrvsMiZMmED//v25//77Q8c+4ogjmDJlCkOHDuXkk09mz549AFxyySWhL/A77riD4uJihg0bxtSpU3F3ACZMmMBPfvITxowZw8CBA/nHP/6x37nW1NRw9dVXM3ToUCZNmsQ3vvGNULn9+vXjjjvuYPz48TzzzDM88sgjFBcXM3LkSM4+++xQAli9ejVHH300xcXF3HLLLaGyw/3sJkyYwDnnnMPgwYO58MILcXfuv/9+NmzYwNe+9jW+9rWv7Rf3yy+/zODBgxk/fjzPPfdcaP2uXbu47LLLKC4upqioiBdffBGAFStWMGbMGAoLCxkxYgQff/wxAH/84x8ZMWIEI0eO5Dvf+Q4AW7Zs4eyzz6a4uJji4mLmz5/f5HW66aab+PTTTyksLOTHP/5xc79eLaJbmFtpUM//JJkxBbkxjkYOJrfPWsEHG7ZHtMwhvTtx6+lDw9q2urqa119/ncsvv3y/9x577DFycnJYuHAhe/fuZdy4cZx88skALF26lJUrV5Kbm0v//v254oorWLBgAffddx8PPPAAv/nNbxg/fjzvvPMOZsajjz7KXXfdxT333APAhx9+yNy5c9mxYweDBg3iqquuAuDjjz9m+vTpPPLII5x33nnMnDmTiy66aJ+4rrnmGn7+858D8J3vfIfZs2dz+umnA1BVVcWCBQuYM2cOt99++35NbM899xxr1qxh2bJlbN68mSOOOILLLrss9H5GRgZvv/02AKWlpUyZMgWAm2++mccee4xrr72W66+/nquuuoqLL744VANsyWf33nvvsWLFCnr37s24ceOYP38+1113Hffeey9z586la9eu+5RVXl7OlClTeOONNzj88MM5//zzQ+/deeedfP3rX+fxxx9n27ZtjBkzhhNPPJGHHnqI66+/ngsvvJCKigqqq6tZsWIFd955J/Pnz6dr165s3boVgOuvv54bbriB8ePHs3btWk455RRWrlzZ6HX61a9+xfLly1myZEmD534glGRaqVdOBtkZKXz0RWS/TERaa8+ePRQWFrJmzRqOPPJITjrppP22efXVV3n//fdDf+mXlZXx8ccfk5aWRnFxMb169QLgsMMOC32BDh8+nLlz5wKB54HOP/98Nm7cSEVFxT7PTkycOJH09HTS09Pp3r17qJZTUFBAYWEhAEceeSRr1qzZL665c+dy1113sXv3brZu3crQoUNDSeass85qct+3336bc889l6SkJHr27LlfraHuF/jy5cu5+eab2bZtGzt37uSUU04BYP78+cycORMIJLmf/OQnLfrsxowZQ58+fQBC12D8+PH7lVHrww8/pKCggAEDBgBw0UUX8fDDD4eO89JLL3H33XcDgYS0du1ajj76aO68805KSko466yzGDBgAG+88QbnnHNOKInl5gb+4H3ttdf44IMPQsfbvn07O3YEWl0au07REpdJxsz6Az8Dctw9Kj2eZsbgntl89IWay2Rf4dY4Iq22T6asrIxJkybx4IMPct111+2zjbvzwAMPhL5ca82bN4/09PTQclJSUmg5KSmJqqoqAK699lpuvPFGzjjjDObNm8dtt90W2qfu/snJyaF96q+vbS6rVV5eztVXX82iRYvo27cvt9122z7PZdTuX7fM+ufUlKysrNC/L7nkEl544QVGjhzJtGnTmDdvXui95m7VDfezayzO+ho7nrszc+ZMBg0atM/6I444grFjx/LXv/6VU045hUcffRR3b7Ccmpoa/vWvf5GZmbnfe62J9UBEtU/GzG4wsxVmttzMpptZRivLedzMNpvZ8gbeO9XMPjKzT8zsJgB3/8zd928riLCBPbL58Isdzf6Si7SlnJwc7r//fu6++24qK/e9MeWUU07h97//fWj9qlWr2LVrV9hll5WVkZ+fD8CTTz4ZkXhrE0rXrl3ZuXNnizvax48fz8yZM6mpqWHTpk37JI76duzYQa9evaisrOSpp54KrR83bhwzZswA2Gd9Xa357LKzs0M1iLoGDx7M6tWrQ31d06dP3+c4DzzwQOh75b333gPgs88+o3///lx33XWcccYZvP/++5xwwgk8/fTTlJaWAoSay04++WR++9vfhspsrhmssTgjIWpJxszygeuA0e4+DEgGLqi3TXczy6637vAGipsGnNrAMZKBB4HTgCHAZDMbEpETCMPgntnsKK/ii+0ar0oOLkVFRYwcOTL0xVnriiuuYMiQIYwaNYphw4bxve99r0V/yd52222ce+65HHvssfv1M7RW586dmTJlCsOHD+fMM8+kuLi4RfufffbZ9OnTJ3Q+Y8eOJScnp8Ftf/GLXzB27FhOOukkBg8eHFp/33338eCDD1JcXExZWVmD+7bms5s6dSqnnXbafk14GRkZPPzww0ycOJHx48dz6KGHht675ZZbqKysZMSIEQwbNix0I8Jf/vIXhg0bRmFhIR9++CEXX3wxQ4cO5Wc/+xnHH388I0eO5MYbbwTg/vvvZ9GiRYwYMYIhQ4bw0EMPNRlnXl4e48aNY9iwYRHv+Mfdo/IC8oF1QC6BZrnZwMn1tjkXeAPICC5PAeY0Ul4/YHm9dUcDr9RZ/inw0zrLzzYX55FHHumt9e/PSv3Qn8z2Nz7c1OoyJD588MEHsQ4hoe3YscPd3b/88kvv37+/b9y4McYRxZeGfr+BRR5GLohaTcbd1wN3A2uBjUCZu79ab5tngJeBGWZ2IXAZcF4LDlObyGqVAPlmlmdmDwFFZvbThnY0s9PN7OHG/moJR+g2ZvXLiMTUpEmTKCws5Nhjj+WWW26hZ8+esQ5JgqLW8W9mXYBvAgXANuAZM7vI3f9Udzt3v8vMZgC/Bw5z950tOUwD69zdS4Erm9rR3WcBs0aPHj2lBcfbR06HVHp2ylCSEYmxpvphJLai2fF/IrDa3be4eyXwHHBM/Y3M7FhgGPA8cGsLj1EC9K2z3AfY0LpwW2eQ7jATEWlUNJPMWuAoM+tggXvsTgBW1t3AzIqARwjUeC4Fcs3sly04xkJggJkVmFkagRsLXopI9GEa3DObT7bspKq6pi0PKyLSLkSzT+bfwLPAu8Cy4LEerrdZB+Bcd//U3WuA7wKf1y/LzKYD/wIGmVmJmV0ePEYVcA3wCoEE9rS7r4jSKTVoYI9sKqpqWFMa/m2gIiKJIqoPY7r7rTTRBObu8+stVxKo2dTfbnITZcwB5hxAmAekdniZD7/YweHds5vZWkQksWiAzAN0ePeOJCcZq9QvIzFUWlpKYWEhhYWF9OzZk/z8fIYPH05hYSFDhgwJravdpqKiImax1g6jD7Bx48b9hri//vrryc/Pp6bmP03Qs2fP5tZbW9plKwcDJZkDlJGaTL+8DnyoJCMxlJeXx5IlS1iyZAlXXnklN9xwA8uWLWPJkiV88MEHoXW126SlpTVaVrSHGanr3nvvDQ1YCYHhUJ5//nn69u3LW2+9FVo/ceJEXnrppXYxbL7sS0kmAgb1zNaQ/9Ku3XbbbUydOpWTTz6Ziy++mPLyci699FKGDx9OUVFRaIDMadOmcc0114T2mzRpUuj24Y4dO/Kzn/2MkSNHctRRR4UGXmxsGH2AmTNncuqp/xnMY+7cuQwbNoyrrrpqn6FWzIwJEyYwe/bsaH0EEiVKMhFweLeOrNu6m4oq3WEm7dfixYt58cUX+fOf/xwa7n7ZsmVMnz6d7373u/sMWtmQXbt2cdRRR7F06VKOO+44Hnkk0L1aO4z+woUL93lIcvXq1XTp0mWfARunT5/O5MmT+da3vsXs2bP3GXtt9OjRDc4nIwe3uByFua0VdMuixmHt1t0c3r1j8ztI3Js3bx5vvvlmaLm2Saj2ixfg+OOPZ8KECdxzzz3s3Bl4BrlXr15MnTqVWbNm8e6774a2vfHGG8nOju6NJWeccUZo1N63336ba6+9FggM5njooYeyatWqJvdPS0sL9a8ceeSR/P3vfwcaH0Z/48aNdOvWLbR/RUUFc+bM4de//jXZ2dmMHTuWV199lYkTJwLQvXt3Nmxo08fgJAKUZCKgX15gKPHVX+5SkhEgMKPjhAkT9lvfUOf1D3/4w/3WnX766aH5VKLhwQcfDCW8OXMCN2fWHRLfGxlZPCUlZZ8O+bq1m9TU1NCw8/WHkG9oOPrMzMx99n/55ZcpKytj+PDhAOzevZsOHTqEkkx5eXmDQ9fLwU3NZRFQ0DXwP+eaL/WsjLQP3//+90M3AfTu3Xu/94877rjQkPerVq1i7dq1DBo0iH79+rFkyRJqampYt24dCxYsaPZYjQ2jP3DgwH0mIZs+fTqPPvooa9asYc2aNaxevZpXX3011Nm/atWq0FTI0n4oyURA5w5pdOmQymdKMhInrr76aqqrqxk+fDjnn38+06ZNIz09nXHjxlFQUMDw4cP50Y9+xKhRo5otq7Fh9LOysjjssMP45JNP2L17N6+88kqo1lL7/vjx45k1axYQuCmg7vvSPlhj1eJEMXr0aF+0aNEBl3PW7+aTlpLEjKlHRyAqaW9WrlzJEUccEesw2p3nn3+exYsX88tfNj2a1KZNm/j2t7/N66+/3kaRSV0N/X6b2WJ3H93cvqrJREhB146s+VL38Iu0xLe+9S369evX7HZr167lnnvuiX5AEnHq+I+Qgq4dmPluObsrquiQpo9VJFxXXHFFs9u0dLZMOXioJhMhBV0Dd5WpNpO4Er3pWeLTgf5eK8lESL+uHYDAbcySeDIyMigtLVWikbji7pSWlpKRkdHqMtSuEyG1z8poyP/E1KdPH0pKStiyZUusQxGJqIyMDPr06dPq/ZVkIiQrPYUendL5bIuSTCJKTU2loKAg1mGIHHTUXBZBBV2zVJMREalDSSaCCrpmqU9GRKQOJZkIKuiaxdZdFZTtrmx+YxGRBKAkE0GhgTLVZCYiAijJRFT/brWjMe+McSQiIgcHJZkI6pvbgSSD1XogU0QEUJKJqPSUZPK7ZKrzX0QkSEkmwgIDZSrJiIiAkkzEFeR1UE1GRCRISSbCeuRksHNvFXsqqmMdiohIzCnJRFheVhoApbv2xjgSEZHYU5KJsC4dAklm666KGEciIhJ7SjIRltdRSUZEpJaSTITlZqUDSjIiItDEUP9mtgxoaAYmA9zdR0QtqnYsN0s1GRGRWk3NJzOpzaKII50yUkhJMkqVZEREGk8y7v55WwYSL8yMLllpfKUkIyLSfJ+MmR1lZgvNbKeZVZhZtZltb4vg2qu8rDTVZERECK/j/7fAZOBjIBO4AnggmkG1d7lZaeqTEREhzLvL3P0TINndq939CeBr0Q2rfeuiJCMiAjTd8V9rt5mlAUvM7C5gI5AV3bDatzwlGRERILyazHeC210D7AL6AmdHM6j2LjcrjbI9lVRW18Q6FBGRmGq2JlPnLrNy4PbohhMfascv+2p3Bd2zM2IcjYhI7IRzd9k4M/u7ma0ys89qX20RXHvVpTbJ7KqMcSQiIrEVTp/MY8ANwGJA49eHIXefkZizYxuMiEgMhZNkytz9b1GPJI7kafwyEREgvCQz18z+H/AcEJokxd3fjVpU7ZzGLxMRCQgnyYwN/hxdZ50DX498OPGhc4dUQElGRCScu8v04GULpSYnkZOZqiQjIgkvnLvLcszsXjNbFHzdY2Y5bRFce6bxy0REwnsY83FgB3Be8LUdeCKaQR0oM+tvZo+Z2bOxiiE3K42tO5VkRCSxhZNkDnP3W939s+DrdqB/czuZ2SAzW1Lntd3MftCaIM3scTPbbGbLG3jvVDP7yMw+MbObAIJxXt6aY0VKl6w0vtqtJCMiiS2cJLPHzMbXLpjZOGBPczu5+0fuXujuhcCRwG7g+brbmFl3M8uut+7wBoqbBpxaf6WZJQMPAqcBQ4DJZjak2TNqA2ouExEJL8lcBTxoZmvM7HMCQ/9f2cLjnAB82sBEaMcDL5pZBoCZTQHur7+zu78FbG2g3DHAJ8GaSwUwA/hmC2OLitzgxGXuDc1gLSKSGJpNMu6+xN1HAiOA4e5e5O5LW3icC4DpDZT9DPAyMMPMLgQuI9DvE658YF2d5RIg38zyzOwhoMjMftrQjmZ2upk9XFZW1oLDhS83K42qGmd7eVVUyhcRaQ8avYXZzC5y9z+Z2Y311gPg7veGc4DgNAFnAA1+2bv7XWY2A/g9gf6fnWHGDmANF+mlNFPbcvdZwKzRo0dPacHxwlb3gcyczNRoHEJE5KDXVE2mds6Y7AZeHVtwjNOAd919U0NvmtmxwDAC/TW3tqBcCNRc+tZZ7gNsaGEZUfGfJLO3mS1FROJXozUZd/9D8J+vufv8uu8FO//DNZkGmsqC5RQBjwATgdXAn8zsl+5+c5hlLwQGmFkBsJ5As9y3WxBb1NSOX1aq25hFJIGF0/H/QJjr9mNmHYCTCIx71pAOwLnu/qm71wDfBerfHICZTQf+BQwysxIzuxzA3asITKb2CrASeNrdV4QTW7R1yQo0kek2ZhFJZE31yRwNHAN0q9cv0wlIDqdwd98N5DXx/vx6y5UEajb1t5vcRBlzgDnhxNOWQjUZ3cYsIgmsqbHL0gj0vaSw76Qo24FzohlUPMhMSyYzNVlP/YtIQmuqT+ZN4E0zm9bA8y0ShtysNA2SKSIJLZw+mUfNrHPtgpl1MbNXohhT3MjNSmOr+mREJIGFk2S6uvu22gV3/wroHr2Q4odqMiKS6MJJMjVmdkjtgpkdSmDSMmlGXlaabmEWkYQWzsyYPwPeNrM3g8vHAVOjF1L8yNVIzCKS4MKZGfNlMxsFHEVgGJcb3P3LqEcWB3I7prG7opo9FdVkpoV117eISFwJp7kMIJ3AKMhlwBAzOy56IcWPnp0yAPhie3mMIxERiY1mazJm9n+B84EVQE1wtQNvRTGuuNC7cyYA67/aQ0HXrGa2FhGJP+H0yZwJDHJ3jfTYQvnBJLNhW7NzvImIxKVwmss+AzRWfSv06JSBGaxXkhGRBBVOTWY3sMTMXgdCtRl3vy5qUcWJtJQkumenqyYjIgkrnCTzUvAlrdC7cyYbypRkRCQxhXML85NtEUi8yu+cyYoN22MdhohITIRzd9lqGnjC3937RyWiOJPfOZNXP9iEu4emrhYRSRThNJeNrvPvDOBcIDc64cSf3p0zqaiqoXRXBV07psc6HBGRNtXs3WXuXlrntd7dfwN8vQ1iiwu9dRuziCSwcJrLRtVZTCJQs8luZHOpp3fnwFP/67/aw4g+nZvZWkQkvoTTXHZPnX9XAauB86ITTvypfSBTz8qISCJqNMmY2fXufh9wi7u/3YYxxZWczFQ6pCWzYZvGLxORxNNUn8ylwZ/3t0Ug8crMyO+cqT4ZEUlITTWXrTSzNUA3M3u/znoD3N1HRDWyOKIHMkUkUTWaZNx9spn1BF4Bzmi7kOJP786ZrNhQFuswRETaXJMd/+7+BTCyjWKJW/mdM/hyZwXlldVkpGryMhFJHOFOWiYHoPZZmY1l6vwXkcSiJNMG6k5eJiKSSMJOMmamqR1bSZOXiUiiajbJmNkxZvYBsDK4PNLMfhf1yOKIJi8TkUQVTk3m18ApQCmAuy8FjotmUPEmLSWJHtkZqsmISMIJq7nM3dfVW1UdhVjiWu/OGXpWRkQSTjhJZp2ZHQO4maWZ2Y8INp1J+Hp3ztTQMiKScMJJMlcC3wfygRKgMLgsLZDfOZP12/bgvt/8byIicSuc6Ze/BC5sg1jiWu3kZZt37KVHp4xYhyMi0ibCmU+mGzAF6Fd3e3e/LHphxZ/ifoHJROcs28il4wpiHI2ISNsIp7nsRSAHeA34a52XtMCQ3p0Y0SeHGQvWqclMRBJGOJOWdXD3n0Q9kgRwQfEh/Pfzy3hv3TZGHdIl1uGIiERdODWZ2Wb2jahHkgDOKOxNh7RkZixYG+tQRETaRKNJxsx2mNl24HoCiWaPmW2vs15aqGN6CqeP6M2spRvZUV4Z63BERKKu0STj7tnu3in4M8ndM+ssd2rLIOPJBWP6sqeympeWboh1KCIiURfO2GWvh7NOwlPYtzODe2YzY0H9QRREROJPU81lGWaWB3Q1sy5mlht89QN6t1WA8cbMuKC4L8vWlzFn2cZYhyMiElVN3V32PeAHBBLKYsCC67cDD0Y5rrh2wZhDeHHpBn4wYwmdO6RyzGFdYx2SiEhUNNUnc5+7FwA/cvf+7l4QfI1099+2YYxxJyM1mScuKebQvA5M/eNilq8vi3VIIiJR0WyfjLs/0BaBJJrOHdL44+Vj6JSRwiVPLGD1l7tiHZKISMRp+uUY6pWTyR8vH0uNw7cfeYd1W3fHOiQRkYhqquN/XPBnetuFk3gO796RP10+lt0V1Ux+5B3NnikicaWpmsz9wZ//aotAEtmQ3p340+VjKdtTybcfeYfN2zXvjIjEh6aSTKWZPQHkm9n99V9tFWCiGN4nhycvG8OWHXuZ8r+LKa/U5KMi0v41lWQmAa8A5QRuYa7/kggbdUgX7rugiPdLtvGjZ5ZqtGYRafcafU4mOFnZDDNb6e5L2zCmhHbSkB785NTB/OpvHzKgezbXnzgg1iGJiLRaOHeXlZrZ82a22cw2mdlMM+sT9cgS2PeO68/Zo/rw69dW8eaqLbEOR0Sk1cJJMk8ALxF48j8fmBVcJ1FiZvzPWcPo0Smd//3X57EOR0Sk1cJJMt3d/Ql3rwq+pgHdohxXwktPSeabhfnM+2gzX+2qiHU4IiKtEk6S2WJmF5lZcvB1EVAa7cAEvlnYm6oa568aSFNE2qlwksxlwHnAF8BG4JzgOomyIb06MbBHR154b32sQxERaZVwxi5b6+5nuHs3d+/u7me6+0HdUWBm/c3sMTN7NtaxHAgz48yifBZ9/hVrSzXkjIi0P1Edu8zMOpvZs2b2oZmtNLOjW1nO48G725Y38N6pZvaRmX1iZjcBuPtn7n75gcZ/MDhjZGDqnheXqDYjIu1PtAfIvA942d0HAyOBlXXfNLPuZpZdb93hDZQzDTi1/kozSyYwt81pwBBgspkNiUzoB4c+XTowpiCXF5as18OZItLuRC3JmFkn4DjgMQB3r3D3bfU2Ox540cwygvtM4T9jpoW4+1vA1gYOMwb4JFhzqQBmAN8MM77TzezhsrKDfy6XbxXl8+mWXSxfvz3WoYiItEizScbMegT7N/4WXB5iZuE0RfUHtgBPmNl7ZvaomWXV3cDdnwFeJjCywIX85yaDcOUD6+oslxAYay3PzB4Ciszspw3t6O6z3H1qTk5OCw4XG98Y1ovUZFOTmYi0O+HUZKYRGMOsd3B5FYFpmZuTAowCfu/uRcAu4Kb6G7n7XQTGR/s9cIa77wyj7FrWwDp391J3v9LdD3P3/9OC8g5KOR1SOW5AN+Ys26gmMxFpV8JJMl3d/WmgBsDdq4BwhgguAUrc/d/B5WcJJJ19mNmxwDDgeeDWcIKud4y+dZb7ABtaWEa78I3hvdhQVs576+q3OIqIHLzCSTK7zCwPcAAzOwpotiPD3b8A1pnZoOCqE4AP6m5jZkXAIwT6US4Fcs3sl+GHz0JggJkVmFkacAGBIXDizolDepCabMx5Xw9mikj7EU6SuZHAF/dhZjYf+CNwbZjlXws8ZWbvA4XA/9R7vwNwrrt/6u41wHeB/Z7BMbPpBCZPG2RmJbV9QsFa1TUEmvNWAk+7+4owY2tXcjJTOXZAN/62/As1mYlIu9HoUP+13P1dMzseGESgD+Qjd68Mp3B3XwKMbuL9+fWWKwnUbOpvN7mJMuYAc8KJp72bOLwXb3y4mSXrtlF0SJdYhyMi0qxmk4yZXVxv1Sgzw93/GKWYpBG1TWZ/fX+jkoyItAvhNJcV13kdC9wGnBHFmKQRajITkfYmnOayffpfzCwH+N+oRSRN+oaazESkHWnNE/+7Ac0JHCMnBZvMnvznmliHIiLSrHD6ZGYRvH2ZQFIaAjwdzaCkcTmZqVx5/GE88MYnjDu8K+eO7tv8TiIiMdJskgHurvPvKuBzdy+JUjwShh+cOJBFa77ilheXM7xPDoN7dop1SCIiDQpnPpk367zmK8HEXnKScd/kQrIzUrn6T++yc29VrEMSEWlQo0nGzHaY2fa7xOvWAAAQ2klEQVQGXjvMTMMBx1j37AwemFzEmtJdTHlyEVt3VcQ6JBGR/TSaZNw92907NfDKdne1zxwEjuqfxz3njWTx2q84/YG3WVZy8E9bICKJJey7y4ITjB1S+4pmUBK+bxX1YeaVxwBw9kP/ZNr81VRV18Q4KhGRgHDmkznDzD4GVgNvAmuAv0U5LmmB4X1ymHXteI7qn8dtsz5g4v1v84+Pt8Q6LBGRsGoyvwCOAla5ewGB0ZTnN72LtLXcrDSevLSYhy4axe7KKr7z2AIum7aQD79Q95mIxE44SabS3UuBJDNLcve5BEZUloOMmXHqsF68duPx/PS0wSxcs5XT7vsHP3x6KSVf7Y51eCKSgMJ5TmabmXUE3iIwbP9mAs/LyEEqPSWZ7x1/GOcX9+V38z5l2j/X8NLS9ZxZmM+VEw7jsG4dYx2iiCQIa26gRTPLIjA9sgEXAjnAU8HaTbs3evRoX7RoUazDiKoN2/bwhzc/ZcbCdVRU13DykB5cPr4/xf26YNbQDNYiIk0zs8Xu3uhULqHtGksyZvZb4M/u/s9IB3cwSYQkU+vLnXt5Yv5q/vTOWsr2VDK0dycuOaYfE0f0okNaOJVaEZGASCSZ6wlMZ9wL+AswPTgJWVxJpCRTa09FNc+/t54n5q/m48076ZiewsThvThndB+OPKQLSUmq3YhI0w44ydQp6FACyeYCIAOYDsxw91WRCDTWEjHJ1HJ3Fq75imcWreOvyzayu6Kanp0yOG14TyYO70XRIV1IVsIRkQZELMnUK7QIeBwY4e7JBxDfQSORk0xdu/ZW8fcPNvHXZRt5c9UWKqpqgpOkdeW4gd04un8efbpkqg9HRIDwk0w4Q/2nAqcSqMmcQOCBzNsPOEI5qGSlp3BmUT5nFuWzo7ySeR9t4a1VW3hz1RZmv78RgJ6dMhjdrwsj+3RmaH4nhvbOISczNcaRi8jBrKk+mZOAycBEYAEwA3jB3Xe1XXjRp5pM09ydVZt2smB1KQvWfMWiNVvZWFYeer9Hp3T6d+1I/25ZHJLbgd6dM+ndOYPu2RnkZqXRIS35gGo/NTXOtj2VlO7cy5ade9myI/jauZetOyvYuquCsj2VoQmPks3o3imdPl060Dc3k68P7k6vnMwD/BREpL5IdPzPBf4MzHT3rRGO76ChJNNypTv3smLDdlZs2M4nm3fy2Zc7+XTzTraX7//4VEZqEp0yUumYnkLHjBQyUpNJS04iNdn26e+pqnEqqmrYW1XDnopqduytZEd5FTvKq6iu2f93NC05idysNHKz0sjJTA2VVVldwxfby9mwbQ+V1U6SwXEDu3H+6L6cNKQHKcmtmQxWROqLSp9MPFKSiZzt5ZVs3FbOhrI9fLljL6W7AjWN7Xsq2bm3ip17qyivrKay2qmsrqGq2qmt5CQnGekpSaSlJJGRkkx2RgrZGankZKaS1zGNvI7pdM1Ko3undLp1zKBTZkqTNaSaGmd16S6ef3c9zy4u4Yvt5YwpyOW3k4vo3imjjT4RkfilJBMmJZn4V13jzHy3hFtfXEFWegoPTC7i6MPyYh2WSLsWbpJR24HEveQk47zRfXnxmnF0ykzhwkff4cUl62MdlkhCUJKRhDGwRzYvXTOeIw/tws3PL2dj2Z5YhyQS95RkJKF0TE/h7nNHUlXj/Pdzy0j05mKRaFOSkYRzaF4WPz5lEHM/2sJz76rZTCSalGQkIV1yTD9GH9qF22etYPP28uZ3EJFWUZKRhJSUZNx1zgj2VtVwz6txMQyfyEFJSUYSVv9uHTlrVB9eXLqest2VsQ5HJC4pyUhCu+ioQyivrGHmuyWxDkUkLinJSEIb2juHokM689S/P9edZiJRoCQjCe+isYfy6ZZd/OuzuJhRXOSgoiQjCW/iiF507pDKU++sjXUoInFHSUYSXkZqMuce2YdXVnyh25lFIkxJRgT49thDqapx/rJwXaxDEYkrSjIiQEHXLMYW5PLS0g2xDkUkrijJiASdNqwnH2/eyadbdsY6FJG4oSQjEnTy0J4AvLLiixhHIhI/lGREgnp3zmRknxxeWbEp1qGIxA0lGZE6ThnWk6XrtmmuGZEIUZIRqeOUYJPZq6rNiESEkoxIHYd168iA7h15ebn6ZUQiQUlGpJ5ThvZkwZqtbN1VEetQRNo9JRmRek4d1pPqGue1lWoyEzlQSjIi9Qzt3Yn8zplqMhOJACUZkXrMjEkje/HWqi18uXNvrMMRadeUZEQacFZRH6pqnFkaZkbkgCjJiDRgUM9shvbuxPPvrY91KCLtWkqsA4i1HTt2cPvtt4eWp0yZAsAjjzwSWnf88cczYcIE7rnnHnbuDIxr1atXL6ZOncqsWbN49913Q9veeOONbNiwgRkzZoTWTZo0iSOPPHKf4wwcOJDJkyczffp0Vq1aFVp/6623snjxYmbPnh1ad8EFF9C7d2/uvffe0LpRo0Zx+umn8/DDD7Nx40YAOnbsyA9/+EPmzZvHm2++qXM6wHMaA6ytzuGTzSNZ+PrsuDineLxOOqfYnFO4LNGnnB09erQvWrQo1mHIQWjzjnKO/j9v8L3j+vNfpw6OdTgiBxUzW+zuo5vbTs1lIo3onp3BsQO68sJ766mpSew/xkRaS0lGpAnfKspnQ1k576wujXUoIu2SkoxIE04e0pOO6SnMXKwbAERaQ0lGpAmZacmcNSqfF5asZ/n6sliHI9LuKMmINOPGkwbSpUMaP31uGVXVNbEOR6RdUZIRaUbnDmncfsZQlq0v44n5a2Idjki7oiQjEoZvDO/JiUd0556/f8Ta0t2xDkek3VCSEQmDmfGLM4eRkpTE9X95TzNnioRJSUYkTL1yMvnV2cNZuXE7J937Fn9653M9PyPSjIQfVkakJSaN6M2I/M7c9Nz73PzCcv7877WMH9CVwr6dGdKrE7kd08hOT8HMYh2qyEFBSUakhQ7J68BTV4zlmcUl/Pnfa5k2fw0Vde46S04yOqankJqcRGqykZJsJFngZQDB/GPQaDKKVIpqqp6lNCi/u3AUA3pkR/UYSjIirWBmnDe6L+eN7sveqmpWbtzBqk07KNtdSdmeSnaUV1JZ41RW1VBV47g7NQ41wbECPfSf/XmTqaEVsTaQTiJ9DGmf0lOSo34MJRmRA5Sekkxh384U9u0c61BEDjrq+BcRkaiJyyRjZv3N7DEzezbWsYiIJLKoJhkzW2Nmy8xsiZm1etIWM3vczDab2fIG3jvVzD4ys0/M7CYAd//M3S8/kNhFROTAtUVN5mvuXtjQ5DZm1t3MsuutO7yBMqYBpzawfzLwIHAaMASYbGZDIhK1iIgcsFg3lx0PvGhmGQBmNgW4v/5G7v4WsLWB/ccAnwRrLhXADOCb4RzYzE43s4fLyjSyrohItEQ7yTjwqpktNrOp+73p/gzwMjDDzC4ELgPOa0H5+cC6OsslQL6Z5ZnZQ0CRmf20wcDcZ7n71JycnBYcTkREWiLatzCPc/cNZtYd+LuZfRislYS4+11mNgP4PXCYu+9sQfkNPU/m7l4KXNn6sEVEJBKiWpNx9w3Bn5uB5wk0b+3DzI4FhgXfv7WFhygB+tZZ7gNsaFWwIiISceYenSd/zSwLSHL3HcF//x24w91frrNNETAdmAisBv4EfObuNzdQXj9gtrsPq7MuBVgFnACsBxYC33b3FS2IcwvweYtPMKAr8GUr923PEvG8E/GcITHPOxHPGVp+3oe6e7fmNopmc1kP4Png2EwpwJ/rJpigDsC57v4pgJl9F7ikfkFmNh2YAHQ1sxLgVnd/zN2rzOwa4BUgGXi8JQkGIJwPqTFmtqihu+biXSKedyKeMyTmeSfiOUP0zjtqScbdPwNGNrPN/HrLlcAjDWw3uYky5gBzWhmmiIhEUaxvYRYRkTimJHNgHo51ADGSiOediOcMiXneiXjOEKXzjlrHv4iIiGoyIiISNUoyrdTQwJzxxsz6mtlcM1tpZivM7Prg+lwz+7uZfRz82SXWsUaamSWb2XtmNju4XGBm/w6e81/MLC3WMUaamXU2s2fN7MPgNT863q+1md0Q/N1ebmbTzSwjHq91Q4MMN3ZtLeD+4Hfb+2Y26kCOrSTTCgk0MGcV8EN3PwI4Cvh+8DxvAl539wHA68HleHM9sLLO8v8Ffh0856+AeBzl+z7gZXcfTODO0JXE8bU2s3zgOmB08Pm7ZOAC4vNaT2P/QYYbu7anAQOCr6kERmNpNSWZ1mn1wJztibtvdPd3g//eQeBLJ5/AuT4Z3OxJ4MzYRBgdZtaHwAPCjwaXDfg6UDs/UTyecyfgOOAxAHevcPdtxPm1JvAYR2bwwe4OwEbi8Fo3MshwY9f2m8AfPeAdoLOZ9WrtsZVkWqfBgTljFEubCI64UAT8G+jh7hshkIiA7rGLLCp+A/wXUBNczgO2uXtVcDker3d/YAvwRLCZ8NHgSB1xe63dfT1wN7CWQHIpAxYT/9e6VmPXNqLfb0oyrdPgwJxtHkUbMbOOwEzgB+6+PdbxRJOZTQI2u/viuqsb2DTerncKMAr4vbsXAbuIo6axhgT7IL4JFAC9gSwCTUX1xdu1bk5Ef9+VZFonYQbmNLNUAgnmKXd/Lrh6U231Ofhzc6zii4JxwBlmtoZAM+jXCdRsOgebVCA+r3cJUOLu/w4uP0sg6cTztT4RWO3uW4KjjTwHHEP8X+tajV3biH6/Kcm0zkJgQPAulDQCnYUvxTimiAv2RTwGrHT3e+u89RLw3eC/vwu82NaxRYu7/9Td+7h7PwLX9Q13vxCYC5wT3CyuzhnA3b8A1pnZoOCqE4APiONrTaCZ7Cgz6xD8Xa8957i+1nU0dm1fAi4O3mV2FFBW26zWGnoYs5XM7BsE/sKtHZjzzhiHFHFmNh74B7CM//RP/DeBfpmngUMI/I96rrs3NHNpu2ZmE4AfufskM+tPoGaTC7wHXOTue2MZX6SZWSGBmx3SgM+ASwn8IRq319rMbgfOJ3An5XvAFQT6H+LqWtcdZBjYRGBalRdo4NoGE+5vCdyNthu41N0XtfrYSjIiIhItai4TEZGoUZIREZGoUZIREZGoUZIREZGoUZIREZGoUZIRiRAz2xn82c/Mvh3hsv+73vI/I1m+SLQoyYhEXj+gRUkmOLJ3U/ZJMu5+TAtjEokJJRmRyPsVcKyZLQnOV5JsZv/PzBYG5+f4HgQe9gzO1/NnAg+8YmYvmNni4BwnU4PrfkVgpOAlZvZUcF1trcmCZS83s2Vmdn6dsufVmR/mqeBDdiJtKqX5TUSkhW4iOFIAQDBZlLl7sZmlA/PN7NXgtmOAYe6+Orh8WfCp60xgoZnNdPebzOwady9s4FhnAYUE5n/pGtznreB7RcBQAuNOzScwLtvbkT9dkcapJiMSfScTGAtqCYEhefIITAgFsKBOggG4zsyWAu8QGKRwAE0bD0x392p33wS8CRTXKbvE3WuAJQSa8UTalGoyItFnwLXu/so+KwNjo+2qt3wicLS77zazeUBGGGU3pu54W9Xo/3eJAdVkRCJvB5BdZ/kV4KrgtAmY2cDghGD15QBfBRPMYAJTXteqrN2/nreA84P9Pt0IzG65ICJnIRIB+stGJPLeB6qCzV7TgPsINFW9G+x830LDU/q+DFxpZu8DHxFoMqv1MPC+mb0bnHqg1vPA0cBSAhNL/Ze7fxFMUiIxp1GYRUQkatRcJiIiUaMkIyIiUaMkIyIiUaMkIyIiUaMkIyIiUaMkIyIiUaMkIyIiUaMkIyIiUfP/AcVa3qdMXJ8cAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.semilogy(log, label='Riemannian gradient descent')\n",
    "plt.axhline(y=val, lw=1, ls='--', color='gray', label='TT-round(A)')\n",
    "plt.xlabel('Iteration')\n",
    "plt.ylabel('Value of the functional')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Exponential machines ([Novikov et al 2017](https://arxiv.org/pdf/1605.03795.pdf))\n",
    "\n",
    "- Linear model depends on vector \n",
    "\n",
    "$$ \\hat{y}_l(x) = \\langle w, x \\rangle + b$$ \n",
    "\n",
    "- Exponential machine depends on tensor $W$\n",
    "\n",
    "$$\\hat{y}_{expm}(x) = \\sum_{i_1 = 0}^1 \\ldots \\sum_{i_d = 0}^1 W_{i_1, \\ldots, i_d} \\prod_{k=1}^d x_k^{i_k}$$\n",
    "\n",
    "- Represent tensor $W$ in TT-format and control TT-rank \n",
    "- Learning of this model is loss minimization with constraint on the tensor rank - riemannien optimization again"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Quantized Tensor Train\n",
    "\n",
    "Consider a 1D array $a_k = f(x_k)$, $k=1,\\dots,2^d$ where $f$ is some 1D function calculated on grid points $x_k$.\n",
    "\n",
    "Let $$k = {2^{d-1} i_1 + 2^{d-2} i_2 + \\dots + 2^0 i_{d}}\\quad i_1,\\dots,i_d = 0,1 $$ \n",
    "be binary representation of $k$, then\n",
    "\n",
    "$$\n",
    "    a_k = a_{2^{d-1} i_1 + 2^{d-2} i_2 + \\dots + 2^0 i_{d}} \\equiv \\tilde a_{i_1,\\dots,i_d},\n",
    "$$\n",
    "\n",
    "where $\\tilde a$ is nothing, but a reshaped tensor $a$. TT decomposition of $\\tilde a$ is called **Quantized Tensor Train (QTT)** decomposition. \n",
    "\n",
    "Interesting fact is that the QTT decomposition has relation to wavelets, more details see [here](https://epubs.siam.org/doi/abs/10.1137/100811647).\n",
    "\n",
    "Contains <font color='red'>$\\mathcal{O}(\\log n r^2)$</font> elements!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Cross approximation method\n",
    "\n",
    "- If decomposition of a tensor is given, then there is no problem to do basic operations fast. \n",
    "\n",
    "- However, the question is if it is possible to find decomposition taking into account that typically tensors even can not be stored. \n",
    "\n",
    "- **[Cross approximation](https://ac.els-cdn.com/S0024379509003747/1-s2.0-S0024379509003747-main.pdf?_tid=0739e049-650f-4fed-8087-61563927f6ac&acdnat=1544040611_b7cd24f6dd7a48d85665d6b562129ab4)** method allows to find the decomposition using only few of its elements (more details in the blackboard)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Tensor networks\n",
    "\n",
    "- Universal way for writing tensor decompositions\n",
    "- Can be interpreted as special neural network architecture\n",
    "- Some theoretical results on DNN expressive power can be derived with tensor network framework, see an example [here](https://arxiv.org/pdf/1711.00811.pdf)\n",
    "- [tntorch](https://github.com/rballester/tntorch) combines tensor networks and PyTorch\n",
    "- More details see [here](https://arxiv.org/pdf/1609.00893.pdf) and [here](https://arxiv.org/abs/1708.09165)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Summary \n",
    "\n",
    "* Tensor decompositions - useful tool to work with multidimensional data\n",
    "* Canonical, Tucker, TT, QTT decompositions\n",
    "* Cross approximation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Next week\n",
    "\n",
    "- Exam: Monday and Tuesday"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>\n",
       "<style>\n",
       "    @font-face {\n",
       "        font-family: \"Computer Modern\";\n",
       "        src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "    }\n",
       "    div.cell{\n",
       "        /*width:80%;*/\n",
       "        /*margin-left:auto !important;\n",
       "        margin-right:auto;*/\n",
       "    }\n",
       "    h1 {\n",
       "        font-family: 'Alegreya Sans', sans-serif;\n",
       "    }\n",
       "    h2 {\n",
       "        font-family: 'Fenix', serif;\n",
       "    }\n",
       "    h3{\n",
       "\t\tfont-family: 'Fenix', serif;\n",
       "        margin-top:12px;\n",
       "        margin-bottom: 3px;\n",
       "       }\n",
       "\th4{\n",
       "\t\tfont-family: 'Fenix', serif;\n",
       "       }\n",
       "    h5 {\n",
       "        font-family: 'Alegreya Sans', sans-serif;\n",
       "    }\t   \n",
       "    div.text_cell_render{\n",
       "        font-family: 'Alegreya Sans',Computer Modern, \"Helvetica Neue\", Arial, Helvetica, Geneva, sans-serif;\n",
       "        line-height: 1.2;\n",
       "        font-size: 120%;\n",
       "        /*width:70%;*/\n",
       "        /*margin-left:auto;*/\n",
       "        margin-right:auto;\n",
       "    }\n",
       "    .CodeMirror{\n",
       "            font-family: \"Source Code Pro\";\n",
       "\t\t\tfont-size: 90%;\n",
       "    }\n",
       "/*    .prompt{\n",
       "        display: None;\n",
       "    }*/\n",
       "    .text_cell_render h1 {\n",
       "        font-weight: 200;\n",
       "        font-size: 50pt;\n",
       "\t\tline-height: 110%;\n",
       "        color:#CD2305;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "    }\t\n",
       "    .text_cell_render h5 {\n",
       "        font-weight: 300;\n",
       "        font-size: 16pt;\n",
       "        color: #CD2305;\n",
       "        font-style: italic;\n",
       "        margin-bottom: .5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "    }\n",
       "    \n",
       "    li {\n",
       "        line-height: 110%;\n",
       "    }\n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "        }  \n",
       "\n",
       "</style>\n",
       "\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"]\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.core.display import HTML\n",
    "def css_styling():\n",
    "    styles = open(\"./styles/custom.css\", \"r\").read()\n",
    "    return HTML(styles)\n",
    "css_styling()"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.13"
  },
  "nav_menu": {},
  "rise": {
   "scroll": true,
   "theme": "sky",
   "transition": "zoom"
  },
  "toc": {
   "navigate_menu": true,
   "number_sections": false,
   "sideBar": true,
   "threshold": 6,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}