{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Alternating least squares for Canonical Polyadic (CP) Decomposition\n",
    "\n",
    "```\n",
    "Copyright 2025 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": [
    "from __future__ import annotations\n",
    "\n",
    "import os\n",
    "import sys\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "import pyttb as ttb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate low-rank data tensor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Choose the rank and shape\n",
    "R = 2\n",
    "tensor_shape = (3, 4, 5)\n",
    "\n",
    "# Set the random seed for reproducibility\n",
    "np.random.seed(0)\n",
    "\n",
    "# Create a low-rank dense tensor from a Kruskal tensor (i.e., ktensor)\n",
    "factor_matrices = [np.random.randn(s, R) for s in tensor_shape]\n",
    "M_true = ttb.ktensor(factor_matrices)  # true solution\n",
    "X = M_true.to_tensor()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run `cp_als` using default parameters\n",
    "The `cp_als` method has two required inputs: a data tensor (`X`) and the rank of the CP model (`R`) to compute.\n",
    "\n",
    "By default, `cp_als` uses a *random* initial guess. At each iteration, it reports the *fit* `f` which is defined for a data tensor `X` and CP model `M` 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 the random seed for reproducibility\n",
    "M1, M1_init, M1_info = ttb.cp_als(X, R)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `cp_als` method returns the following:\n",
    "1. `M1`: the solution as a `ktensor`. \n",
    "2. `M1_init`: the initial guess as a `ktensor` that was used in computing the solution. \n",
    "3. `M1_info`: a `dict` 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 `cp_als` using different initial guesses"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Different random initial guesses can be generated and used by setting different random seeds (via `numpy`). You can also explicitly set the `init` parameter to `\"random\"` to make it clear that the default random initialization is being used, although this is not necessary as illustrated above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(2)  # Set seed for reproducibility\n",
    "M2, M2_init, _ = ttb.cp_als(X, R, init=\"random\")  # leaving third output unassigned"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A specific `ktensor` can also be used as an initial guess. As an example, using the same initial guess (and all other parameters) as the previous run of `cp_als` gives the exact same solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M2_rerun, _, _ = ttb.cp_als(X, R, init=M2_init)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `\"nvecs\"` initialization option computes the initial guess using the eigenvectors of the outer product of the matricized data tensor. This initialization process will require more computation than the default random initialization, but it can often lead to better solutions in fewer iterations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M2_nvecs, _, _ = ttb.cp_als(X, R, init=\"nvecs\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Evaluate and compare the outputs of `cp_als`\n",
    "Use the `ktensor.score()` method to compare outputs of `cp_als` to the true solution if known. In the examples above, the data tensor was generated from a `ktensor`, `M_true`, which can be used to evaluate solutions computed using `cp_als`. A score of 1 indicates a perfect match.\n",
    "\n",
    "\n",
    "Note that the `ktensor.score()` method returns a tuple with the score as the first element and other information related to the score computation as the remaining elements. See the `ktensor` documentation for more information about the return values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "score_M1 = M1.score(M_true)  # not a good solution\n",
    "score_M2 = M2.score(M_true)  # a better solution\n",
    "score_M2_nvecs = M2_nvecs.score(M_true)  # an even better solution\n",
    "\n",
    "print(f\"Score of M1 and M_true: {score_M1[0]}\")\n",
    "print(f\"Score of M2 and M_true: {score_M2[0]}\")\n",
    "print(f\"Score of M2_nvecs and M_true: {score_M2_nvecs[0]}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When two solutions match, as is the case with `M2` and `M2_rerun`, the score will be 1 (up to floating point error)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "score_M2_rerun = M2.score(M2_rerun)\n",
    "\n",
    "print(f\"Score of M2 and M2_rerun: {score_M2_rerun[0]}\")"
   ]
  },
  {
   "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 = 100\n",
    "M2_better, _, _ = ttb.cp_als(X, R, maxiters=more_iters, init=M_true)"
   ]
  },
  {
   "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, printitn=5)"
   ]
  },
  {
   "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, printitn=0)  # No output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Use initial guess based \n",
    "Use the `\"nvecs\"` option to initialize `cp_als` with the leading mode-$n$ singular vectors of the input tensor. This initialization process will require more computation than the default `\"random` initialization in general, but it can lead to better solutions."
   ]
  },
  {
   "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": [
    "M3, _, M3_info = ttb.cp_als(X, R, init=\"nvecs\", printitn=0, dimorder=[2, 1, 0])\n",
    "score_M3 = M3.score(M_true)\n",
    "print(f\"Score of M3 and M_true: {score_M2[0]}\")\n",
    "print(f\"Score of M4 and M_true: {score_M3[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": [
    "M3_rerun, _, M3_rerun_info = ttb.cp_als(X, R, init=\"nvecs\", **M3_info[\"params\"])\n",
    "score_M3_rerun = M3.score(M3_rerun)\n",
    "print(f\"Score of M3 and M3_rerun: {score_M3_rerun[0]}\")\n",
    "\n",
    "print(\"M3_info:\")\n",
    "for k, v in M3_info.items():\n",
    "    print(f\"\\t{k}: {v}\")\n",
    "\n",
    "print(\"M3_rerun_info:\")\n",
    "for k, v in M3_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, using the `stoptol` parameter. Note that you may need to increase the number of iterations for it to converge. Note that you may need to increase the maimum number of iterations allowed, `maxiters`, in order to compute a decomposition whose fit is below the requested stopping tolerance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M4 = 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",
    "        np.array([[1.0, 1.0], [-3.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
}