{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Computing Tucker via the HOSVD\n",
    "\n",
    "```\n",
    "Copyright 2022 National Technology & Engineering Solutions of Sandia,\n",
    "LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the\n",
    "U.S. Government retains certain rights in this software.\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Higher-order Singular Value Decomposition (HOSVD) and Sequentially-truncased HOSVD (ST-HOSVD)\n",
    "The HOSVD computes a Tucker decomposition of a tensor via a simple process. For each mode $k$, it computes the $R_k$ leading left singular values of the matrix unfolding and stores those as factor matrix $U_k$. Then it computes a `ttm` of the original tensor and all the factor matrices to yield the core with shape $(R_1, R_2, \\ldots, R_d)$. The core and factor matrices are used to form the `ttensor`. The values of $R_k$ that lead to a good approximation can be computed automatically to yield a specified error tolerance; this is recommended and the default in our code. The ST-HOSVD is an improvement on the HOSVD that does a TTM in *each* mode before moving on to the next mode. This has the advantage of shrinking the tensor at each step and reducing subsequent computations. ST-HOSVD is the default in the `hosvd` code.\n",
    "\n",
    "* L. R. Tucker, Some mathematical notes on three-mode factor analysis, Psychometrika, 31:279-311, 1966, http://dx.doi.org/10.1007/BF02289464\n",
    "* L. D. Lathauwer, B. D. Moor and J. Vandewalle, A multilinear singular value decomposition, SIAM J. Matrix Analysis and Applications, 21(4):1253-1278, 2000, http://dx.doi.org/10.1137/S0895479896305696\n",
    "* N. Vannieuwenhoven, R. Vandebril and K. Meerbergen, A New Truncation Strategy for the Higher-Order Singular Value Decomposition, SIAM J. Scientific Computing, 34(2):A1027-A1052, 2012, http://dx.doi.org/10.1137/110836067"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyttb as ttb\n",
    "import numpy as np\n",
    "\n",
    "eps_machine = np.finfo(float).eps  # gets machine epsilon for floats, will be used later"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simple example of usage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create random tensor with shape (10, 20, 30) with core with shape (2, 3, 4)\n",
    "def generate_sample_tensor():\n",
    "    np.random.seed(0)\n",
    "    noise = 0.01\n",
    "    core = ttb.tensor(np.random.rand(2, 3, 4), shape=(2, 3, 4))  # The core tensor.\n",
    "    U = [\n",
    "        np.random.rand(10, 2),\n",
    "        np.random.rand(20, 3),\n",
    "        np.random.rand(30, 4),\n",
    "    ]  # The factor matrices.\n",
    "    Soln = ttb.ttensor(core, U)  # Create the solution ttensor.\n",
    "    Z = Soln.full()\n",
    "    Rdm = ttb.tenrand((10, 20, 30))\n",
    "\n",
    "    Data = Z + noise * Z.norm() * Rdm / Rdm.norm()\n",
    "\n",
    "    return Soln, Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "S, X = generate_sample_tensor()\n",
    "\n",
    "# Compute HOSVD with desired relative error = 0.01\n",
    "T = ttb.hosvd(input_tensor=X, tol=0.01)\n",
    "\n",
    "# Check shape of core\n",
    "coreshape = T.core.shape\n",
    "print(f\"Shape of core: {coreshape}\")\n",
    "\n",
    "# Check relative error\n",
    "relerr = (X - T.double()).norm() / X.norm()\n",
    "print(f\"Relative error: {relerr}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate a core with different accuracies for different shapes\n",
    "We will create a core `tensor` that has is nearly block diagonal. The blocks are expontentially decreasing in norm, with the idea that we can pick off one block at a time as we increase the prescribed accuracy of the HOSVD. To do this, we define and use a function `tenrandblk()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def tenrandblk(verbose=True):\n",
    "    np.random.seed(0)\n",
    "    # Block shapes (need not be cubic). Number of rows is the number\n",
    "    # of levels and number of columns is the order of the tensor.\n",
    "    bsz = np.array([[3, 2, 1], [2, 2, 2], [2, 3, 4]])\n",
    "\n",
    "    # Squared norm of each block. Must be length L and sum to <= 1\n",
    "    bns = np.array([0.9, 0.09, 0.009])\n",
    "\n",
    "    # Find the end of each block\n",
    "    bend = np.cumsum(bsz, 0)\n",
    "    # Extract shape: D = # dimensions, L = # levels\n",
    "    D, L = np.shape(bsz)\n",
    "    # Final shape\n",
    "    gsz = bend[-1, :]\n",
    "\n",
    "    ## Create tensor\n",
    "    # Figure out norm of off-block-diagonal\n",
    "    dltnrmsqr = 1 - np.sum(bns)\n",
    "    # Create pattern for off-block-diagonal to be modified as we go\n",
    "    dltpattern = np.ones(tuple(gsz))\n",
    "    # Create tensor to fill in\n",
    "    G = ttb.tenzeros(tuple(gsz))\n",
    "\n",
    "    # Create random entries to use\n",
    "    A = np.sign(np.random.randn(*gsz))\n",
    "    B = 0.1 * np.random.rand(*gsz) + 0.9\n",
    "    Grnd = ttb.tensor(\n",
    "        np.sign(np.random.randn(*gsz)) * (0.1 * np.random.rand(*gsz) + 0.9)\n",
    "    )\n",
    "\n",
    "    # Loop through and create blocks\n",
    "    for i in range(L):\n",
    "        # Figure out ith block pattern\n",
    "        blkrange = []\n",
    "        for k in range(D):\n",
    "            if i == 0:\n",
    "                blkrange.append(np.arange(bend[i, k]))\n",
    "            else:\n",
    "                blkrange.append(np.arange(bend[i - 1, k], bend[i, k]))\n",
    "\n",
    "        # Create pattern that has ones for the block\n",
    "        pattern = np.zeros(tuple(gsz))\n",
    "        ix = np.ix_(blkrange[2], blkrange[0], blkrange[1])\n",
    "        pattern[ix] = 1\n",
    "\n",
    "        # Zero out block in the off-diagonal pattern\n",
    "        dltpattern[ix] = 0\n",
    "\n",
    "        # Randomly fill delta-pattern and rescale\n",
    "        block = Grnd * pattern\n",
    "        sse = (block**2).collapse()\n",
    "        block *= np.sqrt(bns[i] / sse)\n",
    "\n",
    "        # Add to main tensor\n",
    "        G += block\n",
    "\n",
    "        # Verbose output\n",
    "        if verbose:\n",
    "            print(\n",
    "                f\"Created block with shape {tuple(bsz[i,:])} with norm ({block.norm()})^2 = {block.norm()**2}\"\n",
    "            )\n",
    "\n",
    "    if dltnrmsqr > 0:\n",
    "        # final pattern\n",
    "        block = Grnd * dltpattern\n",
    "        sse = (block**2).collapse()\n",
    "        block *= np.sqrt(dltnrmsqr / sse)\n",
    "        G += block\n",
    "        if verbose:\n",
    "            print(\n",
    "                f\"Created tensor with shape {tuple(gsz)} with off-block-diagonal norm ({block.norm()})^2 = {block.norm()**2}\"\n",
    "            )\n",
    "\n",
    "    return G, bsz, bns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create core tensor with given block structure and norm 1\n",
    "G, bsz, bns = tenrandblk()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Shape of G: {G.shape}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate data `tensor` with core described above\n",
    "We take the core `G` and embed into into a larger tensor `X` by using orthogonal transformations. The true rank of this tensor is equal to the shape of `G`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import ortho_group\n",
    "\n",
    "\n",
    "def mat_rand_orth(N, seed=0):\n",
    "    # Generates random n x n orthogonal real matrix.\n",
    "    return ortho_group.rvs(N, seed)\n",
    "\n",
    "\n",
    "# shape of X\n",
    "xsz = np.array([20, 20, 20])\n",
    "\n",
    "# Create orthogonal matrices\n",
    "U = []\n",
    "for k in np.arange(3):\n",
    "    V = mat_rand_orth(xsz[k])[:, : G.shape[k]]\n",
    "    U.append(V)\n",
    "\n",
    "# Create X\n",
    "X = ttb.ttensor(G, U).full()\n",
    "\n",
    "# The norm should be unchanged\n",
    "print(f\"||X|| = {X.norm()}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute (full) HOSVD\n",
    "We compute the ST-HOSVD using the `hosvd` method. We specify the tolerance close to machine precision. Ideally, it finds a core that is the same shape as G."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"ST-HOSVD...\")\n",
    "T = ttb.hosvd(X, 2 * np.sqrt(eps_machine))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute low-rank HOSVD approximation\n",
    "The norm squared of the first two blocks of `G` is 0.99, so specifying an error of `1e-2` *should* yield a core with shape $(4, 4, 3)$(mileage may vary). However, the conservative nature of the algorithm means that it may pick something larger. We can compensate by specifying a larger tolerance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using 1e-2 exactly is potentially too conservative...\n",
    "print(\"Result with tol = sqrt(1e-2):\")\n",
    "T = ttb.hosvd(X, np.sqrt(1e-2), verbosity=1)\n",
    "\n",
    "# But a small multiple (i.e., |ndims(X)|) usually works...\n",
    "print(\"\\nResult with tol = sqrt(3e-2):\")\n",
    "T = ttb.hosvd(X, np.sqrt(3e-2), verbosity=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly, the norm squared of the first block of `G` is 0.9, so specifying an error of `1e-1` *should* result in a core with shape $(3, 2, 1)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using 1e-1 exactly is potentially too conservative...\n",
    "print(\"Result with tol = sqrt(1e-1):\")\n",
    "T = ttb.hosvd(X, np.sqrt(1e-1), verbosity=1)\n",
    "\n",
    "# But a small multiple (i.e., |ndims(X)|) usually works...\n",
    "print(\"\\nResult with tol = sqrt(3e-1):\")\n",
    "T = ttb.hosvd(X, np.sqrt(3e-1), verbosity=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Verbosity - Getting more or less information.\n",
    "Setting the verbosity to zero suppresses all output. Cranking up the verbosity gives some insight into the decision-making process..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = ttb.hosvd(X, tol=np.sqrt(3e-1), verbosity=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = ttb.hosvd(X, tol=3 * np.sqrt(eps_machine), verbosity=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Specify the ranks\n",
    "If you know the rank you want, you can specify it. But there's no guarantee that it will satisfy the specified tolerance. In such cases, the method will throw a warning."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Rank is okay\n",
    "T = ttb.hosvd(X, tol=np.sqrt(3e-1), ranks=bsz[0, :].tolist())\n",
    "\n",
    "# Rank is too small for the specified error\n",
    "T = ttb.hosvd(X, tol=np.sqrt(3e-1), ranks=[1, 1, 1])\n",
    "\n",
    "# But you can set the error to the tensor norm to make the warning go away\n",
    "T = ttb.hosvd(X, tol=X.norm(), ranks=[1, 1, 1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Specify the mode order\n",
    "It's also possible to specify the order of the modes. The default is `np.arange(X.ndims)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reverse_dimorder = np.arange(X.ndims).tolist()[::-1]\n",
    "T = ttb.hosvd(X, tol=np.sqrt(3e-1), dimorder=reverse_dimorder)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate bigger data tensor with core described above\n",
    "Uses the same procedure as before, but now the shape is bigger."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# shape of X\n",
    "xsz = np.array([100, 100, 100])\n",
    "\n",
    "# Create orthogonal matrices\n",
    "U = []\n",
    "for k in np.arange(3):\n",
    "    V = mat_rand_orth(xsz[k])[:, : G.shape[k]]\n",
    "    U.append(V)\n",
    "\n",
    "# Create X\n",
    "Y = ttb.ttensor(G, U).full()\n",
    "\n",
    "# The norm should be unchanged\n",
    "print(f\"||Y|| = {Y.norm()}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ST-HOSVD compared to HOSVD\n",
    "The answers are essentially the same for the sequentially-truncated HOSVD and the HOSVD..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"ST-HOSVD...\")\n",
    "T = ttb.hosvd(Y, tol=2 * np.sqrt(eps_machine))\n",
    "print(\"\\nHOSVD...\")\n",
    "T = ttb.hosvd(Y, tol=2 * np.sqrt(eps_machine), sequential=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But ST-HOSVD may be slightly faster than HOSVD for larger tensors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Time for 5 runs of ST-HOSVD:\")\n",
    "%timeit -n5 T = ttb.hosvd(Y, tol=2 * np.sqrt(eps_machine), verbosity=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"\\nTime for 5 runs of HOSVD:\")\n",
    "%timeit -n5 T = ttb.hosvd(Y, tol=2 * np.sqrt(eps_machine), sequential=False, verbosity=0)"
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 2
}