{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Alternating least squares for Canonical Polyadic (CP) Decomposition\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": [
    "The function `cp_als` computes an estimate of the best rank-$R$ CP model of a tensor $\\mathcal{X}$ using the well-known alternating least-squares algorithm (see, e.g., Kolda and Bader, SIAM Review, 2009, for more information). The input $\\mathcal{X}$ can be almost any type of tensor including a `tensor`, `sptensor`, `ktensor`, or `ttensor`. The output CP model is a `ktensor`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "import pyttb as ttb\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pick the rank and shape\n",
    "R = 3\n",
    "tensor_shape = (6, 8, 10)\n",
    "\n",
    "# Set seed for reproducibility\n",
    "np.random.seed(0)\n",
    "\n",
    "# Create a low-rank dense tensor from a Kruskal tensor\n",
    "factor_matrices = [np.random.rand(s, R) for s in tensor_shape]\n",
    "M_true = ttb.ktensor(factor_matrices)\n",
    "X = M_true.to_tensor()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic call to the method, specifying the data tensor and its rank\n",
    "This uses a *random* initial guess. At each iteration, it reports the *fit* `f` which is defined as \n",
    "```\n",
    "f = 1 - ( X.norm()**2 + M.norm()**2 - 2*<X,M> ) / X.norm()\n",
    "``` \n",
    "and is loosely the proportion of the data described by the CP model, i.e., a fit of 1 is perfect."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute a solution with final ktensor stored in M1\n",
    "np.random.seed(1)  # Set seed for reproducibility\n",
    "few_iters = 10  # Cut off solve early for demo\n",
    "M1, M1_init, M1_info = ttb.cp_als(X, R, maxiters=few_iters)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `cp_als` method returns a the following objects:\n",
    "1. `M1`: the solution as a `ktensor`. \n",
    "2. `M1_init`: the initial guess as a `ktensor` that was generated at runtime since no initial guess was provided. \n",
    "3. `M1_info`: a dictionary containing runtime information with keys:\n",
    "    * `params`: parameters used by `cp_als`\n",
    "    * `iters`: number of iterations performed\n",
    "    * `normresidual`: the norm of the residual `X.norm()**2 + M.norm()**2 - 2*<X,M>`\n",
    "    * `fit`: the fit `f` described above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"M1_info:\")\n",
    "for k, v in M1_info.items():\n",
    "    print(f\"\\t{k}: {v}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run again with a different initial guess, output the initial guess"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute a solution with final ktensor stored in M2\n",
    "np.random.seed(2)  # Set seed for reproducibility\n",
    "M2, M2_init, _ = ttb.cp_als(\n",
    "    X, R, maxiters=few_iters\n",
    ")  # Will not use third output, so leaving unassigned"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Increase the maximum number of iterations\n",
    "Note that the previous run kicked out at only 10 iterations, before reaching the specified convegence tolerance. Let's increase the maximum number of iterations and try again, using the same initial guess."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "more_iters = 10 * few_iters\n",
    "M2_better, _, _ = ttb.cp_als(X, R, maxiters=more_iters, init=M2_init)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compare the two solutions\n",
    "Use the `ktensor` `score()` member function to compare the two solutions. A score of 1 indicates a perfect match."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "score_M2 = M2.score(M_true)\n",
    "score_M2_better = M2_better.score(M_true)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, `score()` returns a tuple with the score as the first element:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Score of M2 and M_true: {score_M2[0]}\")\n",
    "print(f\"Score of M2_better and M_true: {score_M2_better[0]}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See the `ktensor` documentation for more information about the return values of `score()`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rerun with same initial guess\n",
    "Using the same initial guess (and all other parameters) gives the exact same solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M2_rerun, _, _ = ttb.cp_als(X, R, maxiters=few_iters, init=M2_init)\n",
    "score_M2_rerun = M2.score(M2_rerun)  # Score of 1 indicates the same solution\n",
    "print(f\"Score: {score_M2_rerun[0]}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Changing the output frequency\n",
    "Using the `printitn` option to change the output frequency."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = ttb.cp_als(X, R, maxiters=few_iters, printitn=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Suppress all output\n",
    "Set `printitn` to zero to suppress all output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = ttb.cp_als(X, R, maxiters=few_iters, printitn=0)  # No output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Use HOSVD initial guess\n",
    "Use the `\"nvecs\"` option to use the leading mode-$n$ singular vectors as the initial guess. The initialization process will require more computation in general, but it may lead to better solutions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M3, _, _ = ttb.cp_als(X, R, maxiters=few_iters, init=\"nvecs\", printitn=0)\n",
    "score_M3 = M3.score(M_true)\n",
    "print(f\"Score of M2 and M_true: {score_M2[0]}\")\n",
    "print(f\"Score of M3 and M_true: {score_M3[0]}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Change the order of the dimensions in CP\n",
    "Changing the order of the dimensions in which `cp_als` iterates over the input tensor can lead to a different solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M4, _, M4_info = ttb.cp_als(\n",
    "    X, R, maxiters=few_iters, init=\"nvecs\", printitn=0, dimorder=[2, 1, 0]\n",
    ")\n",
    "score_M4 = M4.score(M_true)\n",
    "print(f\"Score of M3 and M_true: {score_M3[0]}\")\n",
    "print(f\"Score of M4 and M_true: {score_M4[0]}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the last example, we also collected the third output argument `M4_info` which has runtime information in it. The field `M4_info[\"iters\"]` has the total number of iterations. The field `M4_info[\"params\"]` has the information used to run the method. Unless the initialization method is `\"random\"`, passing the parameters back to the method will yield the exact same results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M4_rerun, _, M4_rerun_info = ttb.cp_als(X, R, init=\"nvecs\", **M4_info[\"params\"])\n",
    "score_M4_rerun = M4.score(M4_rerun)\n",
    "print(f\"Score of M4 and M4_rerun: {score_M4_rerun[0]}\")\n",
    "\n",
    "print(\"M4_info:\")\n",
    "for k, v in M4_info.items():\n",
    "    print(f\"\\t{k}: {v}\")\n",
    "\n",
    "print(\"M4_rerun_info:\")\n",
    "for k, v in M4_rerun_info.items():\n",
    "    print(f\"\\t{k}: {v}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Change the stopping tolerance\n",
    "It's also possible to loosen or tighten the stopping tolerance on the change in the fit. Note that you may need to increase the number of iterations for it to converge."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M5 = ttb.cp_als(X, R, init=\"nvecs\", maxiters=1000, stoptol=1e-12, printitn=100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Control sign ambiguity of factor matrices\n",
    "The default behavior of `cp_als` is to make a call to `fixsigns()` to fix the sign ambiguity of the factor matrices. You can turn off this behavior by passing the `fixsigns` parameter value of `False` when calling `cp_als`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create rank-2 tensor\n",
    "X2 = ttb.ktensor(\n",
    "    factor_matrices=[\n",
    "        np.array([[1.0, 1.0], [-1.0, -10.0]]),\n",
    "        np.array([[1.0, 1.0], [-2.0, -10.0]]),\n",
    "    ],\n",
    "    weights=np.array([1.0, 1.0]),\n",
    ")\n",
    "print(f\"X2=\\n{X2}\\n\")\n",
    "\n",
    "M_fixsigns, _, _ = ttb.cp_als(X2, 2, printitn=0, init=ttb.ktensor(X2.factor_matrices))\n",
    "print(f\"M_fixsigns=\\n{M_fixsigns}\\n\")  # default behavior, fixsigns called\n",
    "\n",
    "M_no_fixsigns, _, _ = ttb.cp_als(\n",
    "    X2, 2, printitn=0, init=ttb.ktensor(X2.factor_matrices), fixsigns=False\n",
    ")\n",
    "print(f\"M_no_fixsigns=\\n{M_no_fixsigns}\\n\")  # fixsigns not called"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Recommendations\n",
    "* Run multiple times with different guesses and select the solution with the best fit.\n",
    "* Try different ranks and choose the solution that is the best descriptor for your data based on the combination of the fit and the interpretation of the factors, e.g., by visualizing the results."
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 2
}