{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n# Non-negative CP decomposition in Tensorly >=0.6\nExample and comparison of Non-negative Parafac decompositions.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Introduction\nSince version 0.6 in Tensorly, several options are available to compute\nnon-negative CP (NCP), in particular several\nalgorithms:\n\n1. Multiplicative updates (MU) (already in Tensorly < 0.6)\n2. Non-negative Alternating Least Squares (ALS) using Hierarchical ALS (HALS)\n\nNon-negativity is an important constraint to handle for tensor decompositions.\nOne could expect that factors must have only non-negative values after it is\nobtained from a non-negative tensor.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import numpy as np\nimport tensorly as tl\nfrom tensorly.decomposition import non_negative_parafac, non_negative_parafac_hals\nfrom tensorly.decomposition._nn_cp import initialize_nn_cp\nfrom tensorly.cp_tensor import CPTensor\nimport time\nfrom copy import deepcopy"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Create synthetic tensor\nThere are several ways to create a tensor with non-negative entries in Tensorly.\nHere we chose to generate a random from the sequence of integers from 1 to 24000.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Tensor generation\ntensor = tl.tensor(np.arange(24000).reshape((30, 40, 20)), dtype=tl.float32)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Our goal here is to produce an approximation of the tensor generated above\nwhich follows a low-rank CP model, with non-negative coefficients. Before\nusing these algorithms, we can use Tensorly to produce a good initial guess\nfor our NCP. In fact, in order to compare both algorithmic options in a\nfair way, it is a good idea to use same initialized factors in decomposition\nalgorithms. We make use of the ``initialize_cp`` function to initialize the\nfactors of the NCP, and transform these factors (and factors weights) into\nan instance of the CPTensor class:\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"weights_init, factors_init = initialize_nn_cp(tensor, init='random', rank=10)\n\ncp_init = CPTensor((weights_init, factors_init))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Non-negative Parafac\nFrom now on, we can use the same ``cp_init`` tensor as the initial tensor when\nwe use decomposition functions. Now let us first use the algorithm based on\nMultiplicative Update, which can be called as follows:\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"tic = time.time()\ntensor_mu, errors_mu = non_negative_parafac(tensor, rank=10, init=deepcopy(cp_init), return_errors=True)\ncp_reconstruction_mu = tl.cp_to_tensor(tensor_mu)\ntime_mu = time.time()-tic"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here, we also compute the output tensor from the decomposed factors by using\nthe cp_to_tensor function. The tensor cp_reconstruction_mu is therefore a\nlow-rank non-negative approximation of the input tensor; looking at the\nfirst few values of both tensors shows that this is indeed\nthe case but the approximation is quite coarse.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"print('reconstructed tensor\\n', cp_reconstruction_mu[10:12, 10:12, 10:12], '\\n')\nprint('input data tensor\\n', tensor[10:12, 10:12, 10:12])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Non-negative Parafac with HALS\nOur second (new) option to compute NCP is the HALS algorithm, which can be\nused as follows:\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"tic = time.time()\ntensor_hals, errors_hals = non_negative_parafac_hals(tensor, rank=10, init=deepcopy(cp_init), return_errors=True)\ncp_reconstruction_hals = tl.cp_to_tensor(tensor_hals)\ntime_hals = time.time()-tic"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Again, we can look at the reconstructed tensor entries.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"print('reconstructed tensor\\n',cp_reconstruction_hals[10:12, 10:12, 10:12], '\\n')\nprint('input data tensor\\n', tensor[10:12, 10:12, 10:12])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Non-negative Parafac with Exact HALS\nFrom only looking at a few entries of the reconstructed tensors, we can\nalready see a huge gap between HALS and MU outputs.\nAdditionally, HALS algorithm has an option for exact solution to the non-negative\nleast squares subproblem rather than the faster, approximate solution.\nNote that the overall HALS algorithm will still provide an approximation of\nthe input data, but will need longer to reach convergence.\nExact subroutine solution option can be used simply choosing exact as True\nin the function:\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"tic = time.time()\ntensorhals_exact, errors_exact = non_negative_parafac_hals(tensor, rank=10, init=deepcopy(cp_init), return_errors=True, exact=True)\ncp_reconstruction_exact_hals = tl.cp_to_tensor(tensorhals_exact)\ntime_exact_hals = time.time()-tic"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Comparison\nFirst comparison option is processing time for each algorithm:\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"print(str(\"{:.2f}\".format(time_mu)) + ' ' + 'seconds')\nprint(str(\"{:.2f}\".format(time_hals)) + ' ' + 'seconds')\nprint(str(\"{:.2f}\".format(time_exact_hals)) + ' ' + 'seconds')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As it is expected, the exact solution takes much longer than the approximate\nsolution, while the gain in performance is often void. Therefore we recommend\nto avoid this option unless it is specifically required by the application.\nAlso note that on appearance, both MU and HALS have similar runtimes.\nHowever, a closer look suggest they are indeed behaving quite differently.\nComputing the error between the output and the input tensor tells that story better.\nIn Tensorly, we provide a function to calculate Root Mean Square Error (RMSE):\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from tensorly.metrics.regression import RMSE\nprint(RMSE(tensor, cp_reconstruction_mu))\nprint(RMSE(tensor, cp_reconstruction_hals))\nprint(RMSE(tensor, cp_reconstruction_exact_hals))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"According to the RMSE results, HALS is better than the multiplicative update\nwith both exact and approximate solution. In particular, HALS converged to a\nmuch lower reconstruction error than MU. We can better appreciate the difference\nin convergence speed on the following error per iteration plot:\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\ndef each_iteration(a,b,c,title):\n fig=plt.figure()\n fig.set_size_inches(10, fig.get_figheight(), forward=True)\n plt.plot(a)\n plt.plot(b)\n plt.plot(c)\n plt.title(str(title))\n plt.legend(['MU', 'HALS', 'Exact HALS'], loc='upper left')\n\n\neach_iteration(errors_mu, errors_hals, errors_exact, 'Error for each iteration')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In conclusion, on this quick test, it appears that the HALS algorithm gives\nmuch better results than the MU original Tensorly methods. Our recommendation\nis to use HALS as a default, and only resort to MU in specific cases (only\nencountered by expert users most likely).\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## References\n\nGillis, N., & Glineur, F. (2012). Accelerated multiplicative updates and\nhierarchical ALS algorithms for nonnegative matrix factorization.\nNeural computation, 24(4), 1085-1105. (Link)\n\n\n"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"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.12"
}
},
"nbformat": 4,
"nbformat_minor": 0
}