diff --git a/Exercises/Session 2 - Tensor Decomposition (with solutions).ipynb b/Exercises/Session 2 - Tensor Decomposition (with solutions).ipynb deleted file mode 100644 index 0157287..0000000 --- a/Exercises/Session 2 - Tensor Decomposition (with solutions).ipynb +++ /dev/null @@ -1,379 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "id": "e4f44eb7", - "metadata": {}, - "source": [ - "$\\def\\tcoreleft{\\underset{\\tiny\\mid}{\\textcolor{MidnightBlue}{⦸}}}$\n", - "$\\def\\tcorecenter{\\underset{\\tiny\\mid}{\\textcolor{RedOrange}{⦿}}}$\n", - "$\\def\\tcoreright{\\underset{\\tiny\\mid}{\\textcolor{MidnightBlue}{\\oslash}}}$\n", - "

TMQS Workshop 2024 @ Zuse Institute Berlin

\n", - "

Summer School on Tensor Methods for Quantum Simulation

\n", - "

June 3 - 5, 2024

\n", - "

$\\tcoreleft - \\tcoreleft - \\tcoreleft - \\cdots - \\tcorecenter - \\cdots - \\tcoreright - \\tcoreright$

\n", - "
" - ] - }, - { - "cell_type": "markdown", - "id": "58a30939-f570-45cd-a736-d6f21aeb2a0c", - "metadata": {}, - "source": [ - "***" - ] - }, - { - "cell_type": "markdown", - "id": "6c3441f3-b5a9-42c4-ba21-2bc682b0d8ac", - "metadata": {}, - "source": [ - "## **Session 2 - Tensor Decomposition**" - ] - }, - { - "cell_type": "markdown", - "id": "6f6cc702", - "metadata": {}, - "source": [ - "***" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "id": "858a33fe-e16e-4dbb-b369-511d7c71fe71", - "metadata": {}, - "source": [ - "## Exercise 2.1\n", - "\n", - "In general, the question of whether a given tensor can be represented in a low-rank form is very challenging. The following (small!) example shall illustrate this.\n", - "\n", - "**a)**$\\quad$Consider the order-3 tensor $T \\in \\mathbb{R}^{2 \\times 2 \\times 2}$ given by\n", - "\n", - "$\\hspace{1.25cm}$$T_{:, :, 1} = \\begin{pmatrix} 1 & 0 \\\\ 0 & 1\\end{pmatrix}$ and $T_{:, :, 2} = \\begin{pmatrix} 0 & 1 \\\\ 1 & 0\\end{pmatrix}$.\n", - "\n", - "$\\hspace{0.75cm}$Find a canonical decomposition with minimum rank!\n", - "\n", - "$\\hspace{0.75cm}$*Hint:* Both matrices are circulant. That is, they share the same eigenvectors (Fourier modes!).\n", - "\n", - "**b)**$\\quad$*extra task:* Let $T \\in \\mathbb{R}^{n \\times n \\times n}$ with $T_{:,:,k} = S^{k-1}$, where $S = \\begin{pmatrix} 0 & & & 1 \\\\ 1 & 0 & & \\\\ & \\ddots & \\ddots & \\\\ 0 & & 1 & 0 \\end{pmatrix}$. What is the canonical rank of $T$?" - ] - }, - { - "cell_type": "markdown", - "id": "bc6abad4-3b47-447c-9e87-86b44381164b", - "metadata": {}, - "source": [ - "***" - ] - }, - { - "cell_type": "markdown", - "id": "11976a26-d5e7-4ca5-b02b-ce4709527c8c", - "metadata": {}, - "source": [ - "$\\textcolor{red}{\\textbf{SOLUTION:}}$" - ] - }, - { - "cell_type": "markdown", - "id": "742e700e-8485-4923-870b-d34994d2978b", - "metadata": {}, - "source": [ - "**a)**$\\quad$The (normalized) eigenvectors of $T_{:, :, 1}$ and $T_{:, :, 2}$ are given by $v_0 = \\frac{1}{\\sqrt{2}}\\begin{pmatrix}1\\\\1\\end{pmatrix}$ and $v_1 = \\frac{1}{\\sqrt{2}}\\begin{pmatrix}1\\\\-1\\end{pmatrix}$. \n", - "\n", - "$\\hspace{0.75cm}$The eigenvalues of $T_{:, :, 1}$ are $\\lambda_0 = \\lambda_1 =1$ and the eigenvalues of $T_{:, :, 2}$ are $\\mu_0 = 1$ and $\\mu_1 = -1$.\n", - "\n", - "$\\hspace{0.75cm}$That is, $T_{:, :, 1} = \\sum_{i=0}^1 \\lambda_i v_i v_i^\\dagger$ and $T_{:, :, 2} = \\sum_{i=0}^1 \\mu_i v_i v_i^\\dagger$.\n", - "\n", - "$\\hspace{0.75cm}$Since $v_i v_i^\\dagger = v_i \\otimes \\overline{v_i} = v_i \\otimes v_i$, we get\n", - "\n", - "$\\hspace{1.25cm}$$\\displaystyle T = \\begin{pmatrix} 1 & 0 \\\\ 0 & 1\\end{pmatrix} \\otimes \\begin{pmatrix} 1 \\\\ 0 \\end{pmatrix} + \\begin{pmatrix} 0 & 1 \\\\ 1 & 0\\end{pmatrix} \\otimes \\begin{pmatrix}0 \\\\ 1\\end{pmatrix}$\n", - "\n", - "$\\hspace{1.25cm}$$\\displaystyle\\phantom{T} = \\sum_{i=0}^1 \\lambda_i v_i \\otimes v_i \\otimes \\begin{pmatrix} 1 \\\\ 0 \\end{pmatrix} + \\sum_{i=0}^1 \\mu_i v_i \\otimes v_i \\otimes \\begin{pmatrix}0 \\\\ 1\\end{pmatrix}$\n", - "\n", - "$\\hspace{1.25cm}$$\\displaystyle\\phantom{T} = \\sum_{i=0}^1 v_i \\otimes v_i \\otimes \\begin{pmatrix} \\lambda_i \\\\ 0 \\end{pmatrix} + \\sum_{i=0}^1 v_i \\otimes v_i \\otimes \\begin{pmatrix}0 \\\\ \\mu_i\\end{pmatrix}$\n", - "\n", - "$\\hspace{1.25cm}$$\\displaystyle\\phantom{T} = \\sum_{i=0}^1 v_i \\otimes v_i \\otimes \\begin{pmatrix} \\lambda_i \\\\ \\mu_i \\end{pmatrix} $\n", - "\n", - "$\\hspace{1.25cm}$$\\displaystyle\\phantom{T} = \\frac{1}{2} \\begin{pmatrix} 1\\\\1 \\end{pmatrix} \\otimes \\begin{pmatrix} 1\\\\1 \\end{pmatrix} \\otimes \\begin{pmatrix} 1\\\\1 \\end{pmatrix} + \\frac{1}{2} \\begin{pmatrix} 1\\\\-1 \\end{pmatrix} \\otimes \\begin{pmatrix} 1\\\\-1 \\end{pmatrix} \\otimes \\begin{pmatrix} 1\\\\-1 \\end{pmatrix}$\n", - "\n", - "Note that there exists no rank-$1$ decomposition of $T$. This means the canonical rank of $T$ is equal to $2$.\n", - "\n", - "**b)**$\\quad$Again, all \"layers\" $T_{:,:,k}$ have the same set of (complex) eigenvectors $v_0, \\dots, v_{n-1}$ and can be written as \n", - "\n", - "$\\hspace{1.25cm}$$\\displaystyle T_{:,:,k} = \\sum_{i=0}^{n-1} \\lambda_i^{(k)} v_i \\otimes \\overline{v_i}$, \n", - "\n", - "$\\hspace{0.75cm}$where $\\lambda_i^{(k)}$ is the $i$th eigenvalue of $T_{:,:,k}$.\n", - "\n", - "$\\hspace{0.75cm}$Thus, we get \n", - "\n", - "$\\hspace{1.25cm}$$\\displaystyle T = \\sum_{i=0}^{n-1} v_i \\otimes \\overline{v_i} \\otimes \\begin{pmatrix} \\lambda_i^{(1)} \\\\ \\vdots \\\\ \\lambda_i^{(n-1)}\\end{pmatrix}$.\n" - ] - }, - { - "cell_type": "markdown", - "id": "57638421-d883-4f7f-ab0f-aae0505cffca", - "metadata": {}, - "source": [ - "***" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "id": "c7140595-9875-4201-ad30-e186a6722ff8", - "metadata": {}, - "source": [ - "## Exercise 2.2 \n", - "\n", - "We now want to consider a very well-known model in statistical physics: the Ising model. The Ising model is a simple mathematical model used to study phase transitions in magnetic systems. It consists of a lattice arrangement of spins, each of which can point \"up\" or \"down\", and interactions between neighboring spins.\n", - "\n", - "The Hamiltonian function (i.e., the energy of a configuration) for the Ising model can be expressed as follows:\n", - "\n", - "$\\hspace{1cm}$$ \\displaystyle H(\\sigma )=- J \\sum _{\\langle ij\\rangle }\\sigma _{i}\\sigma _{j}-h \\sum _{j} \\sigma _{j}$.\n", - "\n", - "Here, $\\sigma_i$ and $\\sigma_j$ represent the spins at lattice points $i$ and $j$, respectively, where $1 \\leq i,j \\leq d$. $J$ is the interaction constant between spins, $h$ represents the external magnetic field, and the summation is over neighboring lattice points $i$ and $j$.\n", - "\n", - "Since this is a system with nearest-neighbor interactions, we can represent it using a (non-cyclic) SLIM decomposition. \n", - "\n", - "**a)**$\\quad$Find the vectors $S$, $L$, $I$, and $M$ such that the tensor train defined by\n", - "\n", - "$\\hspace{1.25cm}$$\\displaystyle \\mathbf{H} = \\begin{bmatrix} S & L & I\\end{bmatrix} \\otimes \\begin{bmatrix} I & 0& 0\\\\ M & 0& 0\\\\ S & L & I \\end{bmatrix} \\otimes \\dots \\otimes \\begin{bmatrix} I & 0& 0\\\\ M & 0& 0\\\\ S & L & I \\end{bmatrix} \\otimes \\begin{bmatrix} I\\\\ M\\\\ S\\end{bmatrix}$\n", - "\n", - "$\\hspace{0.75cm}$satisfies $\\mathbf{H}_{\\nu_1, \\dots, \\nu_d} = H(\\sigma_1, \\dots, \\sigma_d)$, where\n", - "\n", - "$\\hspace{1.25cm}$$\\nu_i = \\begin{cases} 1 & \\text{if}~\\sigma_i = +1, \\\\ 2 & \\text{if}~\\sigma_i=-1.\\end{cases}$\n", - "\n", - "$\\hspace{0.75cm}$Compare your result with Ising model representation provided by Scikit-TT (use $J=h=1$ and $d=5$):\n", - "\n", - "$\\hspace{1.25cm}$1) Import ```scikit_tt.tensor_train``` as ```tt```.\n", - "\n", - "$\\hspace{1.25cm}$2) Use tt.build_core for constructing the cores of $\\mathbf{H}$, e.g., ```cores[i] = tt.build_core([[I, 0, 0], [M, 0, 0], [S, L, I]])``` for $i=1,2,3$.\n", - "\n", - "$\\hspace{1.25cm}$3) Create the class instance by ```H = TT(cores)```.\n", - "\n", - "$\\hspace{1.25cm}$4) Import ```scikit_tt.models``` as ```mdl``` and load the model from Scikit-TT by ```H_SKTT = mdl.ising(d,J,h)```.\n", - "\n", - "$\\hspace{1.25cm}$5) Compute the error in the Frobenius norm: ```(H-H_SKTT).norm()```.\n", - "\n", - "**b)**$\\quad$*extra task:* The configuration probability, i.e., the probability that the system is in a state $\\sigma$ in equilibrium, is given by\n", - "\n", - "$\\hspace{1.25cm}$$\\pi(\\sigma) = \\frac{e^{- \\beta H(\\sigma)}}{Z}$,\n", - "\n", - "$\\hspace{0.75cm}$where $Z = \\sum_{\\sigma} e^{-\\beta H(\\sigma)}$ is a normalization constant and $\\beta = 1/(k_B \\cdot T)$ is the inverse temperature with Boltzmann constant $k_B$. \n", - "\n", - "$\\hspace{0.75cm}$Suppose $\\beta=1$, can you derive a TT decomposition for $\\pi$?" - ] - }, - { - "cell_type": "markdown", - "id": "83e132fb-95b5-43a8-98ca-0ed7c397767d", - "metadata": {}, - "source": [ - "***" - ] - }, - { - "cell_type": "markdown", - "id": "1cb59312-1033-4dd0-886c-7a034eae661b", - "metadata": {}, - "source": [ - "$\\textcolor{red}{\\textbf{SOLUTION:}}$" - ] - }, - { - "cell_type": "markdown", - "id": "81208cd9-c190-46c9-95f8-d0d23778ebce", - "metadata": {}, - "source": [ - "**a)**$\\quad$The component $S$ represents the single-cell interactions, i.e., $S = -h\\begin{pmatrix}+1 \\\\-1\\end{pmatrix}$.\n", - "\n", - "$\\hspace{0.75cm}$$L \\otimes M$ represents the two-cell (nearest-neighbor) interactions, i.e., $L = -J\\begin{pmatrix}+1 \\\\-1\\end{pmatrix}$ and $M = \\begin{pmatrix}+1 \\\\-1\\end{pmatrix}$.\n", - "\n", - "$\\hspace{0.75cm}$$I$ is simply the vector $\\begin{pmatrix}1 \\\\1\\end{pmatrix}$.\n" - ] - }, - { - "cell_type": "code", - "execution_count": 20, - "id": "21650b1d-183d-403f-80d5-9f27c21e3c69", - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "\n", - "Tensor train with order = 10, \n", - " row_dims = [2, 2, 2, 2, 2, 2, 2, 2, 2, 2], \n", - " col_dims = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1], \n", - " ranks = [1, 3, 3, 3, 3, 3, 3, 3, 3, 3, 1] \n", - "\n", - "error: 1.9946147148160272e-13\n" - ] - } - ], - "source": [ - "import numpy as np\n", - "import scikit_tt.tensor_train as tt\n", - "from scikit_tt.tensor_train import TT\n", - "import scikit_tt.models as mdl\n", - "\n", - "J = 1\n", - "h = 1\n", - "d = 10\n", - "\n", - "S = -h *np.array([1,-1])\n", - "L = -J *np.array([1,-1])\n", - "I = np.array([1,1])\n", - "M = np.array([1,-1])\n", - "\n", - "cores = [None]*d\n", - "cores[0] = tt.build_core([[S, L, I]])\n", - "for i in range(1,d-1):\n", - " cores[i] = tt.build_core([[I, 0, 0], [M, 0, 0], [S, L, I]])\n", - "cores[-1] = tt.build_core([I, M, S])\n", - "\n", - "H = TT(cores)\n", - "print(H, '\\n')\n", - "\n", - "H_SKTT = mdl.ising(d,J,h)\n", - "\n", - "print('error:', (H-H_SKTT).norm())" - ] - }, - { - "cell_type": "markdown", - "id": "4b3280ac-c09b-4d73-a4f7-0e46ab49cbce", - "metadata": {}, - "source": [ - "**b)**$\\quad$The (unnormalized) stationary distribution can be expressed as \n", - "\n", - "$\\hspace{1.25cm}$$\\mathbf{\\Pi} = \\begin{bmatrix} \\begin{pmatrix} e^\\beta \\\\ 0 \\end{pmatrix} & \\begin{pmatrix} 0 \\\\ e^{-\\beta} \\end{pmatrix} \\end{bmatrix} \\begin{bmatrix} \\begin{pmatrix} e^{2\\beta} \\\\ 0 \\end{pmatrix} & \\begin{pmatrix} 0 \\\\ e^{-2\\beta} \\end{pmatrix} \\\\[0.4cm] \\begin{pmatrix} 1 \\\\ 0 \\end{pmatrix} & \\begin{pmatrix} 0 \\\\ 1 \\end{pmatrix} \\end{bmatrix} \\otimes \\ldots \\otimes \\begin{bmatrix} \\begin{pmatrix} e^{2\\beta} \\\\ 0 \\end{pmatrix} & \\begin{pmatrix} 0 \\\\ e^{-2\\beta} \\end{pmatrix} \\\\[0.4cm] \\begin{pmatrix} 1 \\\\ 0 \\end{pmatrix} & \\begin{pmatrix} 0 \\\\ 1 \\end{pmatrix} \\end{bmatrix} \\otimes \\begin{bmatrix} \\begin{pmatrix} e^{2\\beta} \\\\ e^{-2\\beta} \\end{pmatrix} \\\\[0.4cm] \\begin{pmatrix} 1 \\\\ 1 \\end{pmatrix}\\end{bmatrix}$." - ] - }, - { - "cell_type": "markdown", - "id": "179c2c04-c06d-4478-9eca-b662e7fa21ac", - "metadata": {}, - "source": [ - "***" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "id": "2e657e1a-f617-4540-89c8-3a51c6cc7ad6", - "metadata": {}, - "source": [ - "## Exercise 2.3\n", - "\n", - "You've already seen that any arbitrary tensor can be converted into the TT format. Let's try this out practically now. Create a random tensor with order $d$ and mode dimensions $n$:\n", - "\n", - "> T = np.random.random_sample(d*[n])\n", - "\n", - "First, we reshape $T$ into a matrix, thereby defining the \"residual tensor\" $T_{\\text{res}}$:\n", - "\n", - "> T_res = T.reshape([n, n**(d-1)])\n", - "\n", - "Now we need to write a for loop in which, at each step, one mode is separated, and the remaining matrix $T_\\text{res}$ is reshaped.\n", - "\n", - "Assume $T_\\text{res}$ is a $r*n \\times n^e$ matrix with $r \\in \\mathbb{N}$ and $e \\leq d-1$.\n", - "\n", - "To isolate the mode, we apply a (reduced) QR decomposition. The matrix $Q$ then forms the TT core for the separated mode, and the matrix $R$ represents the new residual tensor:\n", - "\n", - "> Q, R = np.linalg.qr(T_res)\n", - "> \n", - "> s = Q.shape[1]\n", - "> \n", - "> core = Q.reshape([r, n, 1, s])\n", - "> \n", - "> T_res = R.reshape([s * n, n**(e-1)])\n", - "\n", - "Try it for various values ​​of d and n, but initially don't choose them too large. What do you notice when you look at the TT ranks?" - ] - }, - { - "cell_type": "markdown", - "id": "a80b8228-1c42-4c43-8c32-6f5a0d6173c2", - "metadata": {}, - "source": [ - "***" - ] - }, - { - "cell_type": "markdown", - "id": "fb078531-7c5e-45f6-b9f2-631d27ad8b7d", - "metadata": {}, - "source": [ - "$\\textcolor{red}{\\textbf{SOLUTION:}}$" - ] - }, - { - "cell_type": "code", - "execution_count": 17, - "id": "b5325236-1967-4cae-8d61-4219d477fd2e", - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "\n", - "Tensor train with order = 10, \n", - " row_dims = [3, 3, 3, 3, 3, 3, 3, 3, 3, 3], \n", - " col_dims = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1], \n", - " ranks = [1, 3, 9, 27, 81, 243, 81, 27, 9, 3, 1]\n" - ] - } - ], - "source": [ - "import numpy as np\n", - "\n", - "d=10\n", - "n=3\n", - "cores = [None]*d\n", - "\n", - "T = np.random.random_sample(d*[n])\n", - "T_res = T.copy()\n", - "r = 1\n", - "for i in range(d-1):\n", - " e = d-i-1\n", - " T_res = T_res.reshape([r*n, n**e])\n", - " Q, R = np.linalg.qr(T_res)\n", - " s = Q.shape[1]\n", - " cores[i] = Q.reshape([r, n, 1, s])\n", - " T_res = R.reshape([s * n, n**(e-1)])\n", - " r = s\n", - "cores[-1] = R.reshape([r, n, 1, 1])\n", - "\n", - "T = TT(cores)\n", - "\n", - "print(T)" - ] - } - ], - "metadata": { - "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.8.10" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -}