\section{INTRODUCTION}
Traditional agriculture across the world is challenged by population pressure, land scarcity, as well as abnormal weather conditions \cite{gebbers2010precision}. One of the more important agricultural advancements in the last 50 years is the concept of the plant factory. Plant factories can be built anywhere, and by stacking upwards, increases the agricultural output of a plot of land many times over traditional farms. Carefully controlled environments through the use of artificial lighting and air-conditioning allow crop production to thrive all year round, with reduced reliance on the outdoors climate when compared to traditional farming methods \cite{brandonnext}.

Smart sensors and machines are increasingly being integrated into plant factories with the development of Internet-of-Things technology. Through analysing the data collected from sensors, control algorithms can be implemented to automate processes and ultimately reduce the cost of production and to improve yield. Control algorithms for plant factories depend on the type of crops being grown and can range from simple rule-based control to the more complex optimization algorithms for decision making. \cite{chieochan2017iot} utilized humidity sensors and the NodeMCU microcontroller to track and control the humidity in a Lingzhi farm according to a rule-based flow chart. \cite{cruz2017design} developed a system to optimize the water and electrical resources of a aquaculture setup based on fuzzy logic control over a wireless sensor network with water level, pH and EC sensors.

A major challenge of growing leafy vegetables in plant factories will be to manage the quantity and quality of artificial lighting to the crops \cite{kim2004green}, while balancing electricity cost. To optimize both plant growth and electricity cost, knowledge of how a plant develop with respect to different lighting or environmental conditions, as well as the power requirements of the artificial lighting, are required. Experiments will have to be conducted to collect such information if they are unavailable beforehand. An example is \cite{morimoto1996intelligent}, where the authors collected data on different tomato plants grown under different conditions, before compressing the information in a feed-forward neural network and optimizing for nutrient delivery schedule using a genetic algorithm.

Also, in recent years, countries like the US \cite{salari2016residential} and Japan \cite{murakami2018energy} are introducing dynamic electricity pricing to the electricity grid in an attempt to smooth out electricity demand, adding complexity to the electrical cost optimization problem for smart plant factories.

As seen in the literature reviewed, deployment of sensing technology and control algorithms are conditional on the optimization objectives, crop grown, and farming setup. In this paper, we would like to explore the effects of red and blue LED lighting on lettuce (Latuca Sativa) grown in a Nutrient Film Technique (NFT) \cite{graves1983nutrient} hydroponics setup. Ultimately, we would also like to optimize the electrical cost of operating a lettuce farm subjected to the plant model and dynamic electrical pricing. Our methodology is summarized in three steps: raw plant growth data collection, plant growth modelling, and finally lighting cost optimization.

The paper made the following contributions:
\begin{itemize}
    \item For plant growth data collection, we proposed an experimental IOT setup and image processing framework to monitor the leaf area growth of lettuce
    \item For plant growth modelling, we modelled the effects of LED lightning on lettuce in a neural network
    \item For lighting cost optimization, we optimized the lighting schedule for operation based on plant model and dynamic electricity pricing
\end{itemize}

While the this paper focuses on optimizing the lightning control of indoor plant factories, the same methodology can be applied to other operational factors like temperature, humidity, and nutrients control. The rest of the paper is organized as such: Section II discusses the experimental setup and data processing methodology, Section III introduces the plant growth model, Section IV describes the lighting cost optimization algorithm and presents the outcome of the optimization. Finally Section V is the conclusion and discusses future work.

\section{SETUP AND IMAGE PROCESSING FRAMEWORK}
\subsection{Experimental Setup}
The experimental setup shown in Fig. \ref{setup} is housed in an air conditioned room where the temperature is kept at 25 $^{\circ}$deg Celsius. It consists of three levels of NFT hydroponics grow-bed, with the lighting conditions of each level individually controlled with a Raspberry Pi 3 micro-controller. The light source for each level is a dimmable panel with integrated red (650-670nm) and blue (450-465nm) LEDs mounted at 40cm above the grow-bed. Each grow-bed can support up to 20 lettuce plants, sharing the same nutrient solution with the other levels. The nutrient solution is circulated through the grow-beds with a pump from a 120 litre tank. The aim of the experiment is to collect information on how the lettuce plants will react to changes in lighting conditions. Hence, the values for nutrient concentration as well as the pH, as measured by the electro-conductivity (EC) sensor and pH sensor sourced from Atlas Scientific, are kept at consistent between 1600-2000 $\mu S/cm$ and $6.4-6.7$. Once the EC, pH, or water-level values dip below the desired threshold, a notification to top up the nutrient solution will be sent out by the micro-controller.
Lighting, temperature and humidity conditions on each test bed are monitored by an array of sensors placed on each test bed. The raw values for the ISL29125 lux sensors used to monitor lighting were calibrated to the standard $\mu mol/m^2s$ units used to measure photosynthetic flux density (PPFD). Temperature and humidity under each were measured by DHT22 sensors. Data from the experimental setup are uploaded to a web server hourly. Prior to seeding of the lettuce plants, the power requirement of the LED panel at each red and blue PPFD settings were quantified in Watts using a power meter. The power requirements are shown in Fig. \ref{power_req}.

\begin{figure}[htb!]\centering \footnotesize
	\includegraphics[width=80mm]{setup.png}
	\caption{Experimental setup of NFT hydroponics system}
	\label{setup}
\end{figure}

\begin{figure}[htb!]\centering \footnotesize
	\includegraphics[width=80mm]{power_req.png}
	\caption{Power requirements of LED lighting panel}
	\label{power_req}
\end{figure}

\subsection{Image Processing}
A top mounted camera is utilized to monitor the plant area throughout the course of each experiment in an attempt to quantify plant growth, as there is a strong correlation between plant leaf area and plant biomass  \cite{weraduwage2015relationship}. Pictures on the growth of the plant were also sent to a web server hourly and the images batch-processed at the end of the experiment. However, it is a challenge for computers to recognize the leaf area for plants in an image. Leaf area segmentation techniques include simple thresholding like Otsu thresholding \cite{kurita1992maximum}, to state of the art convolutional neural network (CNN) based approaches \cite{atanbori2018towards}. In our use case, simple thresholding does not work well with our images as they have too much variations in color, while CNN based techniques require too much effort to retrain the network for our setup. To avoid the time consuming task of manually labelling each image, we modified the K-means algorithm in \cite{dhanachandra2015image} to come up with a helper algorithm in Algorithm \ref{segmentation} for plant segmentation.

\begin{algorithm}
Undistort raw images via homography estimation \cite{dubrofsky2009homography}\\
Cluster images using K-means algorithm via their pixel brightness and distance from center of each pot\\
Estimate the threshold to determine which clusters are plant clusters and which are non-plant clusters \linebreak
(I) Sample a few images in chronological order from a experimental sequence \linebreak
(II) Manually determine which clusters are part of the plant and find the centroid of the cluster which can still be classified as a plant cluster \linebreak
(III) (Assuming plant grow exponentially) Fit an exponential rule to predict the growth of the maximum centroid from data gather in (I) and (II) \\
Apply the exponential rule to all images in the growth stage to obtain segmented images. Final results are shown in Fig. \ref{seg_result}
\caption{Batch Image Segmentation of Leaf Area}
\label{segmentation}
\end{algorithm}

In Fig. \ref{seg_result}, under step 2 to step 3, we show the results of applying Algorithm \ref{segmentation} to images taken from day 3 to day 13 after the lettuce plants are transplants into the setup.

\begin{figure}[htb!]\centering \footnotesize

    \includegraphics[width=85mm]{segmentation.png}
	\caption{Results of Algorithm \ref{segmentation}}
	\label{seg_result}
\end{figure}

\section{PLANT GROWTH MODELLING}%

\begin{table}[]
\caption{Settings for different experiments}
\label{exp_set}
\resizebox{\columnwidth}{!}{
\begin{tabular}{@{}llll@{}}
\toprule
Experiment &
  \begin{tabular}[c]{@{}l@{}}Avg. Red PPFD On\\ $(\mu mol/m^2s)$ \\\end{tabular} &
  \begin{tabular}[c]{@{}l@{}}Avg. Blue PPFD On\\ $(\mu mol/m^2s)$\end{tabular} &
  \begin{tabular}[c]{@{}l@{}}Duty Cycle \\ from 0000hrs\end{tabular} \\ \midrule
1 & 166 & 18 & 18hrs On 6hrs Off  \\
2 & 144 & 54 & 18hrs On 6hrs Off  \\
3 & 160 & 54 & 9hrs On 3hrs Off   \\
4 & 87 & 37 & 13hrs On 11hrs Off \\ \bottomrule
\end{tabular}}
\end{table}

\subsection{Plant Growth Equation}
The relative growth rate (RGR) of a plant's biomass typically follows an exponential relationship \cite{klassen2004real}, as in Eq. \ref{eq1}.

\begin{equation} \label{eq1}
RGR  = ln (M_2 / M_1)/\triangle t
\end{equation}
where $M_2$ is the final dry mass, $M_1$ is the initial dry mass and $\triangle t$ is the change in time in days.
In our experiment, we assumed that the \textbf{increase in leaf area} of a plant is related exponentially to a non-linear function that is dependent on the averaged light energy given to the plant over $\triangle t$, and external variables like EC and pH. The initial leaf area $L_1$ and final leaf area $L_2$ takes the place of $M_1$ and $M_2$, resulting in Eq. \ref{eq2}. $t$ is the time in days since transplanting.

\begin{equation} \label{eq2}
\begin{split}
& f(PPFD_{red_{avg}},PPFD_{blue_{avg}},EC,pH,t) \\
& = ln (L_2 - L_1)/\triangle t
\end{split}
\end{equation}

\subsection{Plant Model}
A neural network is a universal function approximator \cite{sifaoui2008use}, and will allow us to approximate the relationship between variables in Eq. \ref{eq2}. Data from the 4 experiments are split into testing and training data. Data from experiments 1,2,4 goes into the training set, while data from experiment 3 are used for testing. The design of the neural network, as well as the variables used are shown in Fig. \ref{nn}. Variables are min-maxed to 0-1 according to the range displayed, again in Fig. \ref{nn}. The neural network is a two layer feed forward network with 'relu' neurons. Dropout of probability 0.5 is used for regularization to improve the accuracy of the network. The network was trained with back-propagation using the Adam optimizer with Mean Squared Error as the loss function for 2000 epochs.

\begin{figure}[htb!]\centering \footnotesize
    \includegraphics[width=80mm]{neuralnetwork.png}
	\caption{Structure of Neural Network}
	\label{nn}
\end{figure}

A simple linear regression is also used to approximate $f(PPFD_{red_{avg}},PPFD_{blue_{avg}},EC,pH,t)$ for comparison with the neural network. The testing and training errors for both methods are summarized in Table \ref{error_tab} in terms of Mean Squared Error (MSE) and $R^2$ values to determine the goodness of fit of the model. The model is trained on the training data, and tested with either the training data or the testing data.

\begin{table}[htb!]
\caption{Training and testing error for plant model}
\label{error_tab}
\begin{tabular}{@{}lll@{}}
\toprule
\textbf{\begin{tabular}[c]{@{}l@{}}Predicted Leaf Area vs.\\ Actual Leaf Area\end{tabular}} &
  \textbf{\begin{tabular}[c]{@{}l@{}}Linear Regression\\ (MSE)\end{tabular}} &
  \textbf{\begin{tabular}[c]{@{}l@{}}Neural Network\\ (MSE)\end{tabular}} \\ \midrule
\textbf{Training} & 1513.919 & \textbf{126.3719} \\
\textbf{Testing}  & 2017.402 & \textbf{1807.548} \\ \midrule
\textbf{\begin{tabular}[c]{@{}l@{}}Predicted Leaf Area vs.\\ Actual Leaf Area\end{tabular}} &
  \textbf{\begin{tabular}[c]{@{}l@{}}Linear Regression\\ ($R^2$)\end{tabular}} &
  \textbf{\begin{tabular}[c]{@{}l@{}}Neural Network\\ ($R^2$)\end{tabular}} \\ \midrule
\textbf{Training} & 0.8601   & \textbf{0.9902}   \\
\textbf{Testing}  & 0.8783   & \textbf{0.8864}   \\ \bottomrule
\end{tabular}
\end{table}

The neural network model performed much better than that of the linear regression by having a lower MSE and high $R^2$ value on both the training and testing datasets. However the difference in accuracy is lesser in the testing dataset as compared to the training dataset due to a slight distribution shift in the testing dataset.

Fig. \ref{pred_result} shows the shift in distribution and results of the eventual leaf area prediction. With more data, the neural network is expected to perform better than that of the linear regression, with its ability to model the non linear behaviour that may arise with lighting and plant growth interaction .

\begin{figure}[htb!]\centering \footnotesize
    \includegraphics[width=72mm]{pred_result.png}
	\caption{Leaf area prediction}
	\label{pred_result}
\end{figure}

\subsection{Sensitivity Analysis}
Once we have the plant model, we can track the sensitivity of the increase in leaf area to changes in lighting conditions. Eq. \ref{eq2} can be manipulated as such to Eq. \ref{eq3}:

\begin{equation} \label{eq3}
\begin{split}
L_2-L_1 = e^{f(PPFD_{red_{avg}},PPFD_{blue_{avg}},EC,pH,t) \cdot \triangle t}
\end{split}
\end{equation}
Assuming $EC=1800$, $pH=6.5$, and $\triangle t=1 hr$, we investigated the increase in leaf area ($L_2-L_1$) when $t$= $1^{st} day$, $5^{th} day$, and $10^{th} day$ across the range of the light settings. The results are shown in Fig. \ref{growth}.

\begin{figure}[htb!]\centering \footnotesize
    \includegraphics[width=75mm]{growth.png}
	\caption{Growth rate of lettuce with different lighting settings at different stages}
	\label{growth}
\end{figure}

The brighter portions of the heatmap in Fig. \ref{growth} denotes the lighting settings with higher growth rates. As we see, at the early stages of lettuce growth, the plant prefers a mixture of red and blue lighting for best growth. After some time, the preference of the plant switches over to red light. This is consistent in the findings in \cite{zhang2018effects}, who recommends mixing red and blue light in the ratio of 2.2:1 for lettuce growth, in the case of non adjustable LEDs. This ratio might not be the most cost effective option as blue leds will require more energy to provide the same irradiance level as red leds (Fig. \ref{power_req}). For LEDs with adjustable levels, we could leverage on the natural preference of the plant and adjust the lighting conditions for maximum growth and savings at each stage. For best growth on the 5th and 10th day, the ideal red to blue ratio is 14:1 and 20:1 respectively, as suggested by the plant model.

\section{LIGHTING COST OPTIMIZATION}
 One of the concerns for indoor plant factories is electricity pricing, due to the reliance on artificial lighting to drive plant growth. In countries where there are dynamic electricity prices, this can be a bane or a boon. For the growth of vegetative plants which usually do well under a 24 hr lighting schedule, the lighting cycle will have to be planned away from periods when electricity cost are high to prevent huge running cost. At the same time if plant factories can take advantage of the periods of lower-ed electricity costs, they could experience savings in electricity bills. The problem will be even more interesting if solar photovoltics are involved, although this aspect is not investigated in this paper. Prerequisites to optimizing the cost of electricity in a plant factories while maximizing the biomass of the plant depends on three factors: knowledge of the electricity prices, knowledge of how plants react to lighting, and knowledge of the energy consumption of the lighting system.
 In our case study, we attempt to optimize the energy consumption and maximize the biomass of lettuces plants in our setup for a period of 15 days.
 The electricity price plan for reference is taken from the dynamic pricing scheme of Tokyo Electric Power Company \cite{ida2015electricity} and listed in Table \ref{pricing}. The plant model is described in Eq. \ref{eq2} in Section. III and the operating characteristic of the LED lighting panel is illustrated in Fig. \ref{power_req}. This results in the following optimization problem to maximise the profit margin of operation for a single lighting panel and grow-bed:
\begin{equation}
\label{eq4}
\resizebox{80mm}{!}{$
\begin{aligned}
&\underset{s}{max}& P.L_{360}.n - \sum_{\delta=0}^{\delta=360} E(\delta,P(s_\delta)) \\
&\text{s.t.} & L_{\delta+1} - L_{\delta}& = \text{Plant Model in Eq. \ref{eq3}}  \\
& & E(\delta,P(s_\delta)) & = \text{Elec. cost in Table \ref{pricing}} \\
& & P(s_\delta) & =  \text{LED power in Fig. \ref{power_req}} \\
& & L_{360} & > 400
\end{aligned}$}
\end{equation}

In this problem, the time step $\delta$ represents 1 hr in a day, resulting in 360 total timesteps over a period of 15 days.The control sequence $s$ is a 2 by 360 vector, representing the red and blue lighting level at each hour for a period of 15 days.  $L_t$ is the average leaf area of a lettuce at time $\delta$ and $n$ is the number of lettuce cultivated per panel, which is 20. $P$ is the cost of lettuce per unit area per lettuce and will determine the output of the optimization. If P is very low, the optimal behavior of the system will be to not turn on any lights, which is detrimental to the growth of the plants. To prevent this, we also set a constraint on the size of the lettuce in the optimization problem to ensure it is at least 400 $cm^2$ at the end of the 15 days.

\begin{table}[h]
\caption{Dynamic Electricity Pricing in Japan}
\label{pricing}
\begin{tabular}{@{}l|lllll@{}}
\toprule
Time of Day (hrs) &
  \begin{tabular}[c]{@{}l@{}}0000-\\ 0700\end{tabular} &
  \begin{tabular}[c]{@{}l@{}}0700-\\ 1000\end{tabular} &
  \begin{tabular}[c]{@{}l@{}}1000-\\ 1700\end{tabular} &
  \begin{tabular}[c]{@{}l@{}}1700-\\ 2300\end{tabular} &
  \begin{tabular}[c]{@{}l@{}}2300-\\ 0000\end{tabular} \\ \midrule
Cents (Yen) / kWh &
  12 &
  25 &
  38 &
  25 &
  12 \\ \bottomrule
\end{tabular}
\end{table}

\subsection{Genetic Algorithm}
The optimization problem in Eq. \ref{eq4} is a non-linear one due to the presence of non-linear constraints. Among the many other ways to solve non-linear optimization problems, Genetic Algorithm (GA) is reputed as a fast and efficient solver if the fitness function can be quickly evaluated, making it a great algorithm to solve Eq. \ref{eq4}. The first step in a genetic algorithm is to represent the solution of the problem in a chromosome. For our optimization problem, the chromosome is a 360 long vector with each bin taken up by a tuple, $(r_{\delta},b_{\delta})$. Both $(r_{\delta},b_{\delta}) \in Z:\{1 ... 10\}$, with 1 representing 10\% of maximum PPFD and 10 representing 100\% of max PPFD. The maximum PPFD for red LEDs in our system is 200 $\mu mol/m^2s$, while for the blue LEDs it is at 100 $\mu mol/m^2s$. A $(r_{\delta}=3,b_{\delta}=3)$ tuple will mean that the red LEDs are outputting PPFD at 60 $\mu mol/m^2s$ while the blue LEDs are outputting at 30 $\mu mol/m^2s$. The fitness function is the objective function in Eq. \ref{eq4}. Beyond that, the steps in GA like population generation, crossover, and mutation follows that of the canonical version stated in \cite{whitley1994genetic}. Parameters relevant to GA are listed in Table \ref{ga}.

\begin{table}[h]
\caption{GA parameters}
\label{ga}
\begin{tabular}{@{}lllll@{}}
\toprule
\begin{tabular}[c]{@{}l@{}}Population\\ Size\end{tabular} &
  Parents &
  \begin{tabular}[c]{@{}l@{}}Crossover\\ Points\end{tabular} &
  \begin{tabular}[c]{@{}l@{}}Mutations \\ per child\end{tabular} &
  Generations \\ \midrule
100 &
  50 &
  6 &
  20 &
  300 \\ \bottomrule
\end{tabular}
\end{table}

\subsection{Results and Discussions}
Different values of $P$ were tested, $P$ = 0, $P$ = 0.001, and $P$ = 0.01. The generated control sequences from GA were also compared to a $(r_{\delta} = 7,b_{\delta} =7)$ 24hr always-on baseline sequence, as recommended by \cite{zhang2018effects}.

Fig. \ref{convergence} shows the convergence of the fitness function in various cases of $P$.

\begin{figure}[htb!]\centering \footnotesize
    \includegraphics[width=75mm]{convergence.png}
	\caption{Convergence characteristics of GA for different $P$ values}
	\label{convergence}
\end{figure}

When $P$ = 0, the fitness function is dominated by the cost of electricity over the period of 15 days. Hence, this will lead to the control sequences being optimized for a low electricity cost, subjected to constraint that the final leaf area of the plant must be more than 400 $cm^2$. When $P$ is low at 0.001, the revenue from selling the plant will be approximately equal to the cost of electricity, this will induce the GA to search for control sequences that balances the cost of electricity with the final leaf area. Finally when $P$ is high at 0.01, the GA will tend to search for control sequence that favor a larger final leaf area over lowering the cost of electricity.

\begin{figure}[htb!]\centering \footnotesize
    \includegraphics[width=75mm]{control_seq.png}
	\caption{Control sequences for different $P$ values (first 24 hours)}
	\label{control_seq}
\end{figure}

The control sequence for the $1^{st}$ 24 hour sequence in $s$ for various $P$ values are shown in Fig. \ref{control_seq}. We see that when the $P$ is highest, the GA optimizes the growth of the plant by offering more light even when the cost of electricity is high. When $P$ is the lowest, the GA choose to only operate the lights at the portion when electricity prices are the lowest. For $P$ = 0.001, the GA balances the growth of the plant with the cost of electricity.

The final leaf area and electricity cost are summarized in the Table. \ref{comparision}. Profit is the objective function of the problem defined in Eq. \ref{eq4}. When compared to the baseline, the GA allows the system to achieve between 40-52\% savings in electricity cost. Meanwhile, the leaf area output can be improve by up to 6\%.
\begin{table*}[]
\caption{Leaf area and electricity cost comparison for GA with different $P$ values}
\label{comparision}
\resizebox{160mm}{!}{
\begin{tabular}{@{}lllllllllll@{}}
\toprule
\begin{tabular}[c]{@{}l@{}}Plant Price\\ (cents/$cm^2$\\ /plant)\end{tabular} &
  \begin{tabular}[c]{@{}l@{}}Baseline\\ Leaf Area\\ ($cm^2$)\end{tabular} &
  \begin{tabular}[c]{@{}l@{}}Baseline\\ Elec. Cost\\ (cents)\end{tabular} &
  \begin{tabular}[c]{@{}l@{}}Proposed\\ Leaf Area\\ ($cm^2$)\end{tabular} &
  \begin{tabular}[c]{@{}l@{}}Proposed \\ Elec. Cost\\ (cents)\end{tabular} &
  \begin{tabular}[c]{@{}l@{}}Proposed \\ revenue\\ from plants\\ (cents)\end{tabular} &
  \begin{tabular}[c]{@{}l@{}}Profit-Proposed\\ (cents)\end{tabular} &
  \begin{tabular}[c]{@{}l@{}}Profit-Baseline\\ (cents)\end{tabular} &
  \textbf{\begin{tabular}[c]{@{}l@{}}\%\\ Improvement\\ in Leaf Area\end{tabular}} &
  \textbf{\begin{tabular}[c]{@{}l@{}}\%\\ Improvement\\ in Elec. Cost\end{tabular}} &
  \textbf{\begin{tabular}[c]{@{}l@{}}\%\\ Improvement\\ in Profits\end{tabular}} \\ \midrule
\textit{P=0} & 458.0541 & 27.9206 & 442.3375 & 13.5139 & 0        & -13.5139 & -27.9206 & \textbf{-3.4311} & \textbf{51.5985} & -                \\
P=0.001      & 458.0541 & 27.9206 & 452.5732 & 13.5576 & 10.8617  & -2.6958  & -16.9273 & \textbf{1.1965}  & \textbf{51.4422} & -                \\
P=0.01       & 458.0541 & 27.9206 & 486.6086 & 16.5461 & 116.7860 & 100.2399  & 82.0128  & \textbf{6.2338}  & \textbf{40.7387} & \textbf{22.2247} \\ \bottomrule
\end{tabular}}
\end{table*}

\section{CONCLUSION AND FUTURE WORK}
In all, we have observed the growth of lettuce (Latuca Sativa) in our IOT-based experimental setup, and collected the data for analysis. A image processing framework based on the K-means method is used to extract the leaf area information from images captured during the experimental phase. Using the leaf area information and sensor data from the experiments, we modelled the growth of the plant with respect to changes in lighting conditions, pH and EC factors. Finally, we proposed an optimization problem to find the best lighting schedule to reduce cost of production based on our plant growth model. Genetic Algorithm is used to find solutions to the optimization problem. The optimized lighting system is able to achieve between 40-52\% savings in electricity cost while improving the leaf-area of the plants grown by 6\%, in the simulation setting with the plant growth model.

For future work, the optimization problem can be extended into a model-predictive-control problem in a real world experiment. Our framework of collecting data, modelling the phenomena, and optimizing the various factors based on the developed model is generally flexible. We hope to extend this framework into optimizing the other aspects of the plant factory, which includes air-conditioning and nutrient delivery. It will also be interesting to consider other ancillary systems like solar PV or energy storage within the plant factory in the optimization problem.


\section{Introduction}
Extra-large volumes of data are constantly being generated nowadays in areas such as personalized medicine, biology, space, nuclear science, climate and in many other fields. A common way to store and use such data is to first reduce its size without losing important information. Unfortunately, the classical compression techniques that are searching for repeated patterns often cannot be used when the data comes from high-performance computing (HPC) simulations that require a relatively high-precision and need a high accuracy of prediction. Such simulation data is usually high-dimensional and  is naturally represented by tensors, i.e., as multi-dimensional arrays. Tensors in such applications typically represent multiple concurrent \textit{latent} (not-directly observable) processes (e.g., explosions, earth-quakes, etc.) of the simulated phenomenon, imprinting their signatures in various simulated variables (such as temperature, pressure, density, etc.) in different dimensions (e.g., space/time).
Despite the differences, most of the simulated phenomena share the property that the information content of the generated data is quite low, i.e., it can be represented by a low number of parameters, which are however not known a priori and are implicit in the data.
Classical tensor decomposition techniques such as Tucker Decomposition \cite{tucker1966some} and Canonical Polyadic Decompositions (CPD) \cite{hitchcock1927expression}, can extract the latent structures and parameters describing the data, which allows a new type of compression.
For instance, if we consider a $d$-dimensional tensor with $n$ elements in each dimension and a \textit{tensor rank} $r$, CPD allows the representation with the smallest number of parameters, $\mathcal{O}(dnr)$, among all tensor decomposition methods. Unfortunately, finding the canonical rank is an NP-hard problem \cite{haastad1990tensor} and, moreover, the approximation with a fixed rank can be ill-posed \cite{de2008tensor}. For the same tensor, the Tucker decomposition is stable and will have $\mathcal{O}(dnr+r^d)$ number of parameters. Tucker decomposition was previously proposed for compression of large scientific data and it demonstrated an excellent compressibility \cite{austin2016parallel}.  However, for  tensors  of  a  higher  dimension, $d$ (say 10 or more), Tucker decomposition is not feasible since the memory and amount of operations grow exponentially with the tensor dimension $d$.
\begin{figure}[htp]
    \centering
    \raisebox{0.4cm}{\includegraphics[width=0.22\linewidth]{images/tensors.pdf}}\hfill
    \includegraphics[width=0.7\linewidth]{images/tensor2.pdf}
    \caption{(left) A graphical representation of a TT for a 4-dimensional tensor. A core is represented by a circle, while arms indicate the modes of the tensor and rank indices. The first and last TT cores are matrices, whereas the rest are 3-dimensional tensors. (right) Example of a TT for a $5\times 4\times 5\times 6$ tensor with \textit{ tensor train rank} ($r_1=4, r_2=3, r_3=2$).}
    \label{fig:tensors}
\end{figure}

\textit{Tensor train} (TT) has been introduced in \cite{oseledets1} as a method for decomposition of extra large and high-dimensional tensors targeting application in the solution of partial differential equations. TT is a state-of-the-art tensor network based on Singular Value Decomposition (SVD) process that parameterizes the initial high-dimensional tensor by a network of three-dimensional tensors \cite{oseledets2011tensor}. Figure~\ref{fig:tensors} represents a TT decomposition for a four-dimensional tensor. TT (also known in other areas as a matrix product state \cite{li2018shortcut}) is stable and requires only a linear storage in $d$ with $\mathcal{O}(dnr^2)$ parameters, making it suitable for compression of high-dimensional tensors \cite{wang2019principal}.

Many real-world types of data, such as density, temperature, population, probability, etc., are non-negative and hence algorithms that preserve the non-negativity are preferred in order to retain the interpretability and meaning of the compressed data. Nonnegative factorization is used as a model for recovering latent structures in such data and it has some additional useful features.
For instance, with interpretation of the elements of the data as conditional probabilities, one can make use of the duality between latent graphical models and tensor networks \cite{ishteva2015tensors,robeva2019duality} to obtain a graphical model representation of the data, thereby allowing suitable algorithms for graphical models to be applied for tensor analysis. Graphical models have many applications in finance, machine learning, computer vision, speech recognition, and bioinformatics. Moreover, in contrast with the general case, the nonnegative best rank approximation of a tensor \textit{always} exists \cite{lim2009nonnegative} and it is almost always \textit{unique} \cite{qi2016semialgebraic}, which are useful properties for tensor data analysis.

Due to the large input size and the additional computational efforts needed to ensure the nonnegativity of the resulted low-dimensional factors, nonnegative tensor train (nTT) cannot be applied to data coming from many real-life applications, unless massive parallelism is used. Unfortunately, no distributed algorithm for computing an nTT decomposition has been designed yet.
In this paper, we introduce such distributed  nTT algorithm
and analyze its performance.
Our algorithm is based on distributed implementation of nonnegative matrix factorization (NMF) and tensor unfolding and reshaping operations. In the next sections, we describe details of the algorithm, its implementation, and analyze its compression ability and scalability on synthetic and real-world datasets.

\section{Background}
Tensor train (TT) networks were first introduced by Oseledets in \cite{oseledets1}.
Unlike other tensor decompositions, TT, being based on a sequence of SVD operations, is easy to compute and does not suffer from the curse of dimensionality. Based on Eckart–Young theorem \cite{eckart1936approximation}, SVD can be used to provide the best low-rank approximation; hence, we expect that TT (which ignores the nonnegativity) will lead to excellent compression results. From the number of parameters that Tucker decomposition and TT require, we can also conclude that the compresibility of TT will be better that the compresibility of Tucker decomposition. In fact, Figure~\ref{fig:Compr-Tuck-TT} shows this phenomenon, where we present the compression versus relative error results for TT, nTT, Tucker and nonnegative Tucker decompositions, performed for the same synthetically generated tensor calculated with~\cite{novikov2020tensor, kossaifi2019tensorly}. Importantly, although with the smallest number of parameters, CP decomposition uses a single rank for all modes, which makes it less robust compared to TT and Tucker decompositions.

\begin{figure}[h]
    \centering
    \includegraphics[width=\linewidth]{images/compression_vs_error_ttdgen_tt.png}
    \caption{Compression versus relative error for various algorithms on a synthetic data of dimensions 32 x 32 x 32 x 32.}
    \label{fig:Compr-Tuck-TT
\end{figure}

The main factor to be taken into consideration when computing tensor train decompositions is the tradeoff between relative error and the tensor train rank \cite{ding2020tensor}.
Decompositions of specified ranks could have either a large computational complexity, or may result in a large error. A number of algorithms to compute tensor trains exists  (each addressing this issue in a different way), such as, the sequential approach known as TT-SVD \cite{oseledets2011tensor}, Alternating Least Squares (ALS) algorithms \cite{holtz2012alternating}, and Density Matrix Renormalization Group (DMRG) algorithms from quantum physics \cite{huckle2013computations,marti2010complete}.

 A number of distributed TT algorithms \cite{fonal2019distributed}, \cite{wang2019distributed}, \cite{wang2020adtt} exist in the literature. An application of a new out-of-memory algorithm for truncated SVD computation for use in the TT algorithms, targeting tensors so large that they cannot fit in the memory, was also proposed in \cite{9110583}. Recently, a MapReduce-based distributed TT algorithm for large-scale dimension reduction and classification  was introduced in \cite{fonal2019distributed}. However, the proposed framework is unable to achieve significant scalability compared to non-distributed implementation. Next,~\cite{wang2019distributed} used a blocking strategy for distributed TT, where the reshape operation of the tensor block is done locally before decomposition.  They compute the relative error between the TT cores calculated in a distributed manner and the non-distributed approach. In such cases, the paper fails to evaluate the correctness of the framework for large-sized tensors where the non-distributed approach cannot be applied. Recently, \cite{wang2020adtt} proposed a distributed TT for processing Internet of things (IoT) data. Similar to \cite{wang2019distributed}, this framework has the limitations on the ability to compute the final reconstruction error from the factors of a large scale tensor.

Nonnegative tensor train is a much less explored topic,
with only a handful of published works \cite{lee2016nonnegative, shcherbakova2019nonnegative, shcherbakova2019nonnegative2}.

In \cite{lee2016nonnegative}, Lee \textit{et. al.} proposed a nonnegative tensor train and provide an NTT-HALS algorithm to compute it. NTT-HALS aims to minimize the Frobenius norm reconstruction error for a given tensor train rank by using a Hierarchical Alternating Least Squares (HALS) procedure. The authors demonstrated that the storage cost of their  NTT-HALS decomposition is significantly lower than a corresponding nonnegative Tucker decomposition. Due to its computational complexity, the NTT-HALS algorithm is ill suited to large high-dimensional tensors. In \cite{shcherbakova2019nonnegative2} the authors propose an algorithm called NTTF, which relies on successive unfoldings and NMF approximations. That algorithm can accommodate a desired relative error by ensuring that each NMF approximation achieves a certain accuracy. This is computationally efficient approach, but it can result in unbalanced and large tensor train ranks. Finally, in \cite{shcherbakova2019nonnegative}, Shcherbakova proposes a NTT-MU algorithm similar in goal to \cite{lee2016nonnegative}, but based on a DMRG algorithm to additionally minimize the tensor train ranks. DMRG-based algorithms simultaneously optimize over consecutive pairs of tensors in the train, this allows the tensor train rank between the pair to be adjusted. However, the proposed algorithm utilizes multiplicative update to ensure nonnegativity, which suffers from the inability to converge to a high accuracy.

\begin{table}[htp]\footnotesize
	\centering
	\caption{Notation} \label{tab:notations}
	\begin{tabular}{|c|c|l|}
		\toprule
		\textbf{Notation} & \textbf{Dimensions} & \textbf{Description} \\
		\midrule
		$\ten{A}$ & $n_1 \times \hdots \times n_{d}$ & Tensor of interest\\
		$\ten[i][]{G}$ & $r_{i-1} \times n_i \times r_i$ & Tensor train core\\
		$p_i$ & scalar & Processor count along each dimension \\
		$\ten[i_1,\dots,i_d]{A}$ & $\frac{n_1}{p_1} \times \hdots \times \frac{n_{d}}{p_{d}}$ & Distributed block of the tensor  \\
		$m$   & scalar & Number of rows of a matrix \\
		$n$   & scalar & Number of columns of a matrix \\
		$\mat{X}$ & $m \times n$ & Input matrix to NMF \\
		$p_r$ & scalar & Proc. count along rows of matrix\\
		$p_c$ & scalar & Proc. count along columns of matrix\\
		$\mat[i,j]{X}$ & $\frac{m}{p_r} \times \frac{n}{p_c}$ & Distributed block of matrix \\

		$r$   & scalar & Low rank \\
		$\mat{W}$ & $m \times r$ & Left low rank factor \\
		$\mat{H}$ & $r \times n$ & Right low rank factor \\
		$p$ & scalar & Total processor count, $p = \prod_{i=1}^{d} p_i $ \\
		$\mat{(W^i)^{j}}$ & $\frac{m}{p} \times r$ & Left low rank factor on $(i,j)^{th}$ \\
		&&processor\\
		$\mat{(H^j)^{i}}$ & $r \times \frac{n}{p}$ & Right low rank factor on $(i,j)^{th}$  \\
		&&processor\\

		$\mat[i]{W}$ & $\frac{m}{p_r} \times r$ & Block of factor W  \\
		&& corresponding  to block $\mat[i,j][]{X}$ \\
		$\mat[j]{H}$ & $r \times \frac{n}{p_c}$ & Block of factor H  \\
		&& corresponding to block $\mat[i,j][]{X}$ \\
		\bottomrule
	\end{tabular}
\end{table}

\section{Algorithmic details}
The tensor-train format is an efficient representation of a higher dimensional tensor in  terms of storage requirements and computational robustness \cite{oseledets1},\cite{oseledets2}. Tensor train representation achieves such performance metrics by combining the major advantages of the  Canonical format  and the Tucker format \cite{kolda2009tensor}. Tensor train decomposes a $d$-dimensional tensor $\ten{A}\in \mathbb{R}^{n_1 \times \hdots \times n_d}$ into \textit{d} 3-dimensional tensors $\ten[i]{G} \in \mathbb{R}^{r_{i-1}\times n_{i}\times r_{i}}$, where $r_{0} = r_{d} = 1$ (so $\ten[1]{G}$ and $\ten[d]{G}$ are actually matrices) and $r_{k} \geq 1$ for $k = 1,...,d-1$,  such that

\begin{equation}
\begin{aligned}
    \ten{A} = \ten[1]{G} \circ \ten[2]{G} \circ \hdots \circ \ten[d]{G},
    \end{aligned}
\end{equation}
where $\ten{A} \circ \ten{B}$ is used to define tensor-tensor multiplication with contraction along the last axis of $\ten{A}$ and the first axis of $\ten{B}$,
\linebreak
$(\ten{A} \circ \ten{B})_{i_1,i_2,\hdots,i_{d-1},j_2,j_3,\hdots,j_d} = \sum_{k} \ten[][i_1,\hdots,i_{d-1},k]{A} \ten[][k,j_1,j_2,\hdots,j_d]{B}$. Here, the 3-dimensional tensors $\ten[i]{G}$ are called  \textit{TT cores} and the numbers $r_{1},\dots,r_{d-1}$ are called  \textit{TT ranks}.
Based on the above, any element of the tensor $\ten{A}$ can be further be represented as

\begin{equation}
\begin{aligned}
\ten[][i_{1},\dots,i_{d}]{A} = \sum_{k_1, \hdots, k_{d-1}}^{r_1, \hdots, r_{d-1}} \ten[1][i_1,k_1]{G} \ten[2][k_1,i_2,k_2]{G} \hdots \ten[d][k_{d-1},i_d]{G}.
\end{aligned}
\end{equation}

Table~\ref{tab:notations} shows the notations used in the paper.

\begin{figure}[htp]
    \centering
    \includegraphics[scale=0.6]{images/diagram3.pdf}
    \caption{The decomposition procedure of a tensor into a TT format iteratively isolates the next mode and applies NMF to get the next factor. The TT cores are represented by circles, while the remaining unfactored matrix is represented by a bar.}
    \label{fig:illustration}
\end{figure}

\subsection{Non-Negative Tensor Train}

\begin{algorithm}[htp]
    \caption{$\mat[i,j]{X}$ = $\operatorname{distReshape} (\ten[i_i,\hdots,i_d]{A},[m,n], [p_{i},p_{j}])$ -- Distributed reshaping of tensor $\ten{A}$ into $\mat{X}$}. \label{alg:distreshape}
    \begin{algorithmic}[1]
    \Require: Tensor $\ten[i_1,i_2,\hdots,i_d]{A}\in \mathbb{R}^{\frac{n_{1}}{p_{1}} \times \hdots \times \frac{n_{d}}{p_d}}$,target matrix shape $m \times n$ and target processor grid size $p_{r} \times p_{c}$.
    \State  Each MPI rank writes a block of $\ten[i_1,i_2,\hdots,i_d]{A}$.
    \State Perform global reshaping of the tensor from tensor dimensions $n_{1} \times n_{2} \times \hdots \times n_{d}$ into $m \times n$  i.e $\ten{X}=reshape(\ten{A}, [m_{1},m_{2},...m_{k}]$).
    \State Compute in-memory data $\mat[i,j]{X}$ for each MPI rank from the reshaped tensor $\mat{X}$ where size of $\mat[i,j]{X} = \frac{m}{p_r} \times \frac{n}{p_c}$
    \Ensure  $\mat{X} \in \mathbb{R}^{ m \times n }$ and $\mat[i,j]{X} \in \mathbb{R}^{\frac{m}{p_r} \times \frac{n}{p_c}}$
     \end{algorithmic}
\end{algorithm}

\begin{algorithm}[htp]
    \caption{$\ten[][1]{G}, \ten[][2]{G}, \hdots, \ten[][d]{G}$ = $\operatorname{distnTT}(\ten[i_1,i_2,\hdots,i_d][]{A},\epsilon)$ -- Distributed non negative tensor train algorithm}\label{alg:distTT}
	\begin{algorithmic}[1]
    \Require  Tensor $\ten[i_1,i_2,\hdots,i_d]{A}\in \mathbb{R}^{\frac{n_{1}}{p_{1}} \times \hdots \times \frac{n_{d}}{p_d}}$ and threshold $\epsilon$.
    \State Set $r_{0}$ = $r_{d}$ = 1
    \State Set $\mat[i,j]{A} = \operatorname{distReshape}$ $(\ten[i_1,i_2,\hdots,i_d]{A},[n_1,\frac{S }{n_1}],[p_1,p/p_1])$  \Comment{$S=n_1 \times n_2 \hdots \times n_d$ and $p=p_1 \times p_2 \hdots \times p_d$}
    \For{$l$ in 1 to $d-1$}

    \State $\mat[i,j]{X} = \operatorname{distReshape}$ $(\mat[i,j]{A},[r_{l-1} n_{l},\frac{S}{r_{l-1} n_{l}}],[p_1,p/p_1])$ \Comment{$S=n_l \times n_{l+1} \hdots \times n_d$}
    \State $\mat[i,j]{U} \mat{\Sigma} \mat{V^{(i,j)T}}$ =  $\operatorname{distSVD}( \mat[i,j]{X})$  \Comment{$\Sigma \in \mathbb{R}^{N \times N}$}
    \State Choose $r_l$ smallest $k$ such that $\frac{\sqrt{\sigma^{2}_{k+1}+ ... + \sigma^{2}_{N}}}{\sqrt{\sigma^{2}_{1}+ ... + \sigma^{2}_{N}}} \leq \epsilon$
    \State $\mat{(W^i)^{j}}$, $\mat{(H^j)^{i}} =\operatorname{distBCDnmf}(\mat[i,j]{X}, r_{l})
    \State $\ten{G}^{(l)}=\operatorname{reshape}(\operatorname{all\_gather}(\mat{(W^i)^{j}}), [r_{l-1},n_{l},r_{l}])$ \Comment{Perform global all\_gather on $\mat{(W^i)^{j}}$ along row-wise distribution. $\mathbb{R}^{r_{l-1}\times n_{l}\times r_{l}}$}
    \State $\mat[i,j]{A}=\mat{(H^j)^{i}}$ \Comment{$\in \mathbb{R}^{r_{l} \times n_{l+1} ... n_{d}}$  is 1D-distributed}
    \EndFor
    \State Set d-th core to $\ten[d]{G}_{:,:,1} = \operatorname{all\_gather}(\mat[i,j]{A})$ \Comment{Perform global all\_gather on $\mat[i,j]{A}$ along column-wise distribution}
    \Ensure Approximation $\ten[1]{G} \circ \ten[2]{G} \circ \hdots \circ \ten[d]{G}$ of $\ten{A}$ in TT format with cores $\ten[1]{G}, \ten[2]{G}, \hdots, \ten[d]{G}$ and TT ranks $r_{0},r_{1},....,r_{d}$
    \end{algorithmic}

\end{algorithm}
Figure~\ref{fig:illustration} illustrates the construction of a $d$-dimensional TT. Given a $d$-dimensional tensor $\ten{A} \in \mathbb{R}^{n_1 \times n_2 \times ... \times n_d}$, the first step of TT is to isolate the first mode of $\ten{A}$ with left unfolding to produce a matrix $\mat{X} \in \mathbb{R}^{n_1 \times  n_2 n_3  \dots n_d}$. NMF is used to reduce $\mat[][]{X}$ to its corresponding factors $\mat{W} \in \mathbb{R}^{n_1 \times r_1}$ and $\mat{H} \in \mathbb{R}^{r_1 \times n_2 n_3  ... n_d}$. The rank of this NMF decomposition, $r_1$, is selected by an SVD-based heuristic, namely, $r_1$ is selected by choosing it to be equal to the smallest $k$ such that $\sqrt{\sigma^{2}_{k+1}+ ... + \sigma^{2}_{N}}/\sqrt{\sigma^{2}_{1}+ ... + \sigma^{2}_{N}} \leq \epsilon$, where $\sigma_{i}$ is the $i^{th}$ singular value of $\mat{X}$ and $N = \min(n_1,n_2 n_3  \dots n_d )$. The singular values decrease and reach a plateau, separating a cluster of non-latent dimensions ($r \geq r_{s+1}$) to that of latent ones ($r \leq r_{s}$). The number of these larger singular values is used to calculate the effective rank of a  matrix \cite{roy2007effective} as well in Bayesian PCA \cite{bishop1999bayesian}. There are other sophisticated approaches to identify latent dimensionality of NMF \cite{distnmfk}, which could be included in our method. The left NMF factor ${W}$ is the first core of the non-negative tensor train decomposition of \ten[][]{X}. The right NMF factor $\mat{H}$ is further reshaped to $\mat{X} \in \mathbb{R}^{r_1n_2 \times n_3 \times n_4 \times... \times n_d}$. The procedure is repeated until the final NMF factor $\mat{H} \in \mathbb{R}^{r_{d-1} \times n_d}$, i.e., core $\ten[d][]{G}$ of nTT is obtained.

\begin{algorithm}[htp]\footnotesize
    \caption{$\mat{(W^i)^{j}}$, $\mat{(H^j)^{i}}$ = $\operatorname{distBCDnmf}(\mat[i,j]{X}, r)$ -- Distributed BCD NMF algorithm}\label{alg:distbcd}
	\begin{algorithmic}[1]
    \Require  $\mat[i,j]{X} \in \mathbb{R}_{+}^{\frac{m}{p_r} \times \frac{n}{p_c}}$ and desired rank $r$.

    \State Initialize $\mat{(W^i)^{j}}$, $\mat{(H^j)^{i}}$ = $\operatorname{rand}(\frac{m}{p}, r)$, $\operatorname{rand}(r, \frac{n}{p})$

    \State $\mat{(W^i)^{j}}_m$, $\mat{(H^j)^{i}}_m$ = $\frac{\mat{(W^i)^{j}}}{||\mat[]{W}||}\sqrt{||\mat[]{X}||}$, $\frac{\mat{(H^j)^{i}}}{||\mat[]{H}||}\sqrt{||\mat[]{X}||}$ \Comment{Normalize}
    \State $\mat{{HH^T}}$, $((\mat{{XH^T}})^i)^j$ =$\operatorname{distMM^T}(\mat{(H^j)^{i}_m)}, \operatorname{distXH^T}(\mat[i,j]{X},\mat{(H^j)^{i}_m)}$

    \State $t$, $obj$ = $1.0$, $\frac{1}{2}||X||$ \Comment{Correction of $\mat{W}$ and $\mat{H}$}
    \For{$l$ in $max\_iters$}
        \Statex \textbf{/* Update W given H */}

        \State $\mat{{(W_{m}{HH^T}^i)^j}}$ = $\mat{(W^i)^{j}}_m \mat[]{{HH^T}}$ \Comment{$i^{th}$ proc $\mathbb{R}_{+}^{\frac{m}{p} \times r}$}
        \State $\mat{(G_{W}^i)^j}$  = $\mat{{(W_{m}{HH^T}^i)^j}}$ -- $((\mat{{XH^T}})^i)^j$ \Comment{$\mat{(G_{W}^i)^j} \in \mathbb{R}_{+}^{\frac{m}{p} \times r}$}
        \State $\mat{(W^i)^{j}}$ = $\operatorname{max}(0, \frac{\mat{(W^i)^{j}_m}}{||\mat{{HH^T}}||})$
        \State $\mat{(W^i)^{j}} /= ||\mat{W}||_1$ \Comment{Normalize with $L_1$ norm }
        \State $\mat{W^TW}$ = $\operatorname{distMM^T}(\mat{((W^i)^{j}})^T)$
        \Statex \textbf{/* Update H given W*/}

        \State $\mat{(W^TWH_m^j)^i}$ = $\mat{W^TW} \mat{(H^j)^{i}}_m$ \Comment{$\mat{(W^TWH_m^j)^i} \in \mathbb{R}_{+}^{r \times \frac{n}{p} }$}
        \State $((\mat{{W^TX}})^j)^i$ = $\operatorname{distWTX}(\mat[i,j]{X}, \mat{(W^i)^{j}})$
         \State $\mat{(G_H^j)^i}$  = $\mat{(W^TWH_m^j)^i}$ -- $((\mat{{W^TX}})^j)^i$ \Comment{$\mat{(G_H^j)^i} \in \mathbb{R}_{+}^{r \times \frac{n}{p}}$}
        \State $\mat{(H^j)^{i}}$ = $\operatorname{max}(0,\mat{(H^j)^{i}_m} - \frac{\mat{(G_H^j)^i}}{ \sqrt{WTW}})$
        \State $\mat{{HH^T}}$ = $\operatorname{distMM^T}(\mat{(H^j)^{i}})$
        \State $((\mat{{XH^T}})^i)^j$ = $\operatorname{distXH^T}(\mat[i,j]{X},\mat{(H^j)^{i}})$
        \If{$\frac{1}{2}||\mat{X}-\mat{WH}||^2 >= obj$} \Comment{\textbf{/* Correction */}}\label{bcd:start_correct}
            \State $Initialize \; \mat{(W^i)^{j}}$ and $\mat{(H^j)^{i}}$
            \State $\mat{{HH^T}}$ = $\operatorname{distMM^T}(\mat{(H^j)^{i}})$
            \State $((\mat{{XH^T}})^i)^j$ = $\operatorname{distXH^T}(\mat[i,j]{X},\mat{(H^j)^{i}})$ \label{bcd:end_correct}
        \Else{} \Comment{\textbf{/* Extrapolation */}}\label{bcd:start_extrapolate}
            \State $w$ = $(t-1)/\frac{1+\sqrt{1+4t^2}}{2}$
            \State $w_W$ = $\operatorname{min}(w, \delta \sqrt{||\mat{HH^T_{(l-1)})}||/||\mat{{HH^T}}||})$
            \State $w_H$ = $\operatorname{min}{w, \delta \sqrt{||\mat{W^TW_{(l-1)}}||/||\mat{{W^TW}}||}}$
            \State $\mat{(W^i)^{j}}_m$ = $\mat{(W^i)^{j}}$ + $w_W$*($\mat{(W^i)^{j}}$--$\mat{(W^i)^{j}_{l-1}}$)
            \State $\mat{(H^j)^{i}}_m$ = $\mat{(H^j)^{i}}$ + $w_H$*($\mat{(H^j)^{i}}$--$\mat{(H^j)^{i}_{l-1}}$)
            \State  $t,\; obj$ = $\frac{1+\sqrt{1+4t^2}}{2},\; \frac{1}{2}||\mat{X}-\mat{WH}||^2$ \label{bcd:end_extrapolate}
        \EndIf
    \EndFor
    \State \Return $\mat{(W^i)^{j}}$, $\mat{(H^j)^{i}}$
    \Ensure $\mat{(W^i)^{j}} \in \mathbb{R}_{+}^{\frac{m}{p} \times r}$ and $\mat{(H^j)^{i}} \in \mathbb{R}_{+}^{r \times \frac{n}{p}}$ and $\mat{W},\mat{H} \approx \operatorname*{argmin}_{\widetilde{\mat{W}} \geq 0,\widetilde{\mat{H}} \geq 0} ||\mat{X}-\widetilde{\mat{W}} \widetilde{\mat{H}}||_{F}^{2}$
    \end{algorithmic}
\end{algorithm}

\begin{figure*}
    \centering
    \includegraphics[width=\textwidth]{images/final_overview3.pdf}
    \caption{Overview of distributed tensor-train decomposition of a 4D tensor with a processor grid of size 2x2x2x2.}
    \label{fig:overviewl}
\end{figure*}

\subsection{Distribution Strategy}
Figure~\ref{fig:overviewl} shows the distributed TT for a 4-dimensional tensor $\ten{A}$ into four cores $\ten[1][]{G}$, $\ten[2][]{G}$ , $\ten[3][]{G}$ and  $\ten[4][]{G}$. We choose a 4-dimensional processor grid of size $2\times 2\times 2\times 2$ that divides each mode of the tensor. Considering the size of tensor $\ten{A}$ is $n_{1}\times n_{2}\times n_{3}\times n_{4}$, each distributed block of tensor $\ten[i_1,i_2,i_3,i_4][]{A}$ will have a size of  $\frac{n_{1}}{2} \times \frac{n_{2}}{2} \times \frac{n_{3}}{2} \times \frac{n_{4}}{2}$. We first  perform a distributed reshaping/unfolding of the tensor $\ten{A}$ into a matrix $\mat{X}$ of size $n_{1} \times n_{2}n_{3}n_{4}$ using Zarr and Dask packages. We use a Zarr shared file system to store the tensor and the intermediate factors. Dask operates on the Zarr file object for reshaping. Dask first performs global tensor reshape operation via a lazy evaluation/call-by-need approach and then each MPI rank computes the in-memory chunk of data afterward.  The algorithm for distributed reshape is presented in Algorithm~\ref{alg:distreshape}. Each block of reshaped  matrix $\mat{X}$ is of size $\frac{n_{1}}{2} \times \frac{n_{2}}{2} \frac{n_{3}}{2} \frac{n_{4}}{2}$. Next, we carry out the distributed SVD to find the TT rank. Then, we perform the distributed NMF algorithm with block co-ordinate descent (BCD) optimization \cite{xu2013block} presented in Algorithm \ref{alg:distbcd} on the distributed matrix $\mat{X}$. The distribution strategy of the matrix and factors for the distributed NMF is following  \cite{distnmfk}. The dist-NMF produces distributed factors matrices $\mat{W}$ and $\mat{H}$. The matrix $\mat{W}$ is gathered across the processor cores to obtain the first core tensor $\ten[1][]{G}$. Next, the \mbox{1-d} distributed factor matrix $\mat{H}$ undergoes a distributed reshape as per Algorithm~\ref{alg:distreshape} to obtain a 2D-distributed matrix $\mat{X}$. The decomposition is carried out as described above, followed by obtaining factors $\mat{W}$ and $\mat{H}$, which are reshaped accordingly, to obtain cores $\ten[i][]{G}$ and matrices $\mat{X}$, respectively, until the final decomposition is constructed.

Algorithm~\ref{alg:distTT} describes the distributed nonnegative tensor train. Along with the distributed reshape and the SVD, the next key component is the distributed BCD algorithm (see Algorithm~\ref{alg:distbcd}). The inputs to $\operatorname{distBCDnmf}$ algorithm are a 2D distributed input matrix ($\mat[i,j]{X})$ and the rank value ($r$). The algorithm returns intermediate factors, $\mat{(W^i)^j}$ and $\mat{(H^j)^i}$, where the second matrix is used in the reshape operation of the TT algorithm for further decomposition until the last mode of the tensor is processed. In order to get optimal intermediate factors, we employed alternating nonnegative least squares (ANLS) strategy in BCD. This alternates updating each factor while fixing the other factor to be constant.

\begin{algorithm}[htp]
    \caption{$\mat{{MM^T}}$ = $\operatorname{distMM^T}(\mat{(M^i)^{j}})$ -- Distributed Gram calculation of $\mat{(M^i)^{j}}$}\label{alg:distMMt}
	\begin{algorithmic}[1]
    \Require  $\mat{(M^i)^{j}}$ or $\mat{((M^i)^{j}})^T$
    \Statex \textbf{/* To calculate $\mat{{HH^T}}$ or $\mat{{W^TW}}$ */}
    \State $\mat[i,j]{{U}}$ = $\mat{(M^i)^{j}}\mat{((M^i)^{j}})^T$  \Comment{$\mat[i,j]{{U}} \in \mathbb{R}_{+}^{r \times r}$}
    \State $\mat{{MM^T}}$ = $\sum \mat[i,j]{U}$ \Comment{$all\_reduce$ across all proc}

    \Ensure $\mat{{MM^T}} \in \mathbb{R}_{+}^{r \times r}$
    \end{algorithmic}
\end{algorithm}

\begin{algorithm}[htp]
    \caption{ $((\mat{{XH^T}})^i)^j$ = $\operatorname{distXH^T}(\mat[i,j]{X},\mat{(H^j)^{i}})$ -- Distributed matrix multiplication of $\mat[i,j]{X}$ and $\mat{(H^j)^{i}}$} \label{alg:distXHt}
	\begin{algorithmic}[1]
    \Require  $\mat[i,j]{X} \in \mathbb{R}_{+}^{\frac{m}{p_r} \times \frac{n}{p_c}}$ and  $\mat{(H^j)^{i}} \in \mathbb{R}_{+}^{r \times \frac{n}{p}}$

    \State $\mat[j]{H}$ = $\operatorname{all\_gather}(\mat{(H^j)^{i}})$ \Comment{across processor columns, $\mat[j]{H} \in \mathbb{R}_{+}^{r \times \frac{n}{p_c}}$}
    \State $\mat[i, j]{{V}}$ = $\mat[i, j]{X}\mat{H^{(j)T}}$  \Comment{$\mat[i, j]{{V}} \in \mathbb{R}_{+}^{\frac{m}{p_r} \times r}$}
    \State compute $(\mat{{XH^T}})^i$ = $\sum \mat[i, j]{V}$ \Comment{\textit{reduce-scatter} on processor rows for row-wise distribution}
     \State $(i,j)^{th}$ processor holds $((\mat{{XH^T}})^i)^j$ after \textit{reduce-scatter}
    \Ensure $((\mat{{XH^T}})^i)^j \in \mathbb{R}_{+}^{\frac{m}{p} \times r}$
    \end{algorithmic}
\end{algorithm}

\begin{algorithm}[htp]
    \caption{$((\mat{{W^TX}})^j)^i$ = $\operatorname{distWTX}(\mat[i,j]{X}, \mat{(W^i)^{j}})$ -- Distributed matrix multiplication of $\mat[i,j]{X}$ and $\mat{(W^i)^{j}}$} \label{alg:distWtX}
    \begin{algorithmic}[1]
        \Require $\mat[i,j]{X} \in \mathbb{R}_{+}^{\frac{m}{p_r} \times \frac{n}{p_c}}$ and $\mat{(W^i)^{j}} \in \mathbb{R}_{+}^{\frac{n}{p} \times r}$
    \State $\mat[i]{W}$ = $\operatorname{all\_gather}(\mat{(W^i)^{j}})$ \Comment{across processor rows, $\mat[i]{W} \in \mathbb{R}_{+}^{ \frac{n}{p_r}} \times r$}

        \State compute $\mat[i,j]{Y}$ = $\mat{W^{(i)T}} \mat[i,j]{X}$ \Comment{$\mat[i,j]{Y} \in \mathbb{R}_{+}^{r \times \frac{n}{p_c}}$}
        \State compute  $(\mat{{W^TX}})^j$ = $\sum \mat[i, j]{Y}$ \Comment{\textit{reduce-scatter} on processor columns for column-wise distribution}
          \State $(i,j)^{th}$ processor holds $((\mat{{W^TX}})^j)^i$ after \textit{reduce-scatter}
        \Ensure $ (\mat{{W^TX}})^j)^i \in \mathbb{R}_{+}^{r \times \frac{n}{p}}$
    \end{algorithmic}
\end{algorithm}

The main computational steps in the distributed BCD are the calculation of matrices $\mat{W^{T}}\mat{W}$ or $\mat{H}\mat{H^{T}}$ (Gram matrices), $\mat{X}\mat{H^{T}}$, and $\mat{W^{T}}\mat{X}$. Each of these computations are being performed multiple times in the BCD algorithm. Therefore, we present the distributed Gram calculation in Algorithm~\ref{alg:distMMt}, Algorithm~\ref{alg:distXHt} computes distributed $\mat{X}\mat{H^{T}}$, and Algorithm~\ref{alg:distWtX} computes distributed $\mat{W^{T}}\mat{X}$. We note the final conditions of BCD try to guarantee convergence to an optimal solution. In this step, the  $\mat{(W^i)^j}$ and $\mat{(H^j)^i}$ matrices are initialized to the initial values (lines~\ref{bcd:start_correct}--\ref{bcd:end_correct} in Algorithm~\ref{alg:distbcd}) when the objective value of optimization is worse than the previous iteration. Otherwise, the intermediate factors are updated (lines~\ref{bcd:start_extrapolate}--\ref{bcd:end_extrapolate} in Algorithm~\ref{alg:distbcd}) accordingly with the use of a user defined hyper parameter ($\delta$). We now evaluate our approach.

\section{Experiments and Results}
\label{sec:experiments}
We run the experiments on the HPC cluster {\em Grizzly}, located at Los Alamos National Laboratory (LANL). {\em Grizzly} has Intel Xeon Broadwell (E5-2695v4) processors with a total of $1490$ compute nodes, where each node has $18$-core dual socket Ivy Bridge processor. Each of the $36$ processors has a clock speed of $2.1$ GHz with a private $L_1$ and $L_2$ caches of sizes $64$KB and $256$KB. Both the sockets share an $L_3$ cache of size $45$MB, where each node contains $128$GB of memory. Grizzly uses Tri-Lab Operating System Stack (TOSS) version 3, while the interconnect is Intel {\em OmniPath} that uses a fat-tree topology.

Our source code is in Python, where the dependencies include Dask~\footnote{\href{https://docs.dask.org/en/latest/}{https://docs.dask.org/en/latest/}}, Numpy~\footnote{\href{https://numpy.org/doc/1.18/reference/index.html}{https://numpy.org/doc/1.18/reference/index.html}}, MPI4PY~\footnote{\href{https://mpi4py.readthedocs.io/en/stable/}{https://mpi4py.readthedocs.io/en/stable/}}, and Zarr~\footnote{\href{https://zarr.readthedocs.io/en/stable/}{https://zarr.readthedocs.io/en/stable/}}. Our framework supports dense tensors, and we leave the sparse implementations for future releases. We use Python (v3.7.0) compiler and the OpenMPI (v2.1.2) library available on Grizzly.

\subsection{Data Generation}
We generate a synthetic tensor, say $\ten{A} \in \mathbb{R}^{n_{1},...n_{d}}$ , with known tensor-train ranks $r_{1},\dots,r_{d-1}$, and selected dimensions $n_{1},....n_{d}$. Each of the tensor train factors, e.g., $\ten[][i]{G}$, gets elements sampled from a uniform distribution between 0 and 1. The tensor $\ten{A}$ is then generated as a product of the TT factors and is distributed among the processors if its size is too large. Specifically, if the tensor is too large, we first generate the TT cores as described above, reshape them into matrices, and distribute them along the 1D processor grid. We then perform a distributed matrix multiplication of the factors and repeat the process until we obtain the final unfolded tensor. Then, we apply distributed reshape operation on the unfolded tensor to obtain the desired $d$-dimensional tensor.

\subsection{Scalability}
We evaluate the scalability of the distributed tensor train algorithm looking at both strong and weak scaling. We also analyze the scaling performance with respect to tensor train ranks. We used 16, 32, 64, 128, and 256 processors/cores for these experiments. For the scaling experiments, we collect the total time taken for the nonnegative matrix factorization for all $d-1$ factors, having fixed the number of NMF iterations to 100. This is done ten times and the average times are reported. We also report the time spent on compute, communication, and I/O operations in the TT algorithm. The computation costs are: \emph{GR}---the  local  computations  of a Gram  matrix  ($\mat{W^{T}}\mat{W}$ or $\mat{H}\mat{H^{T}}$),  which  is  of  size  $r \times r$ ; \emph{MM}—matrix-matrix  multiplications  using  local  (MPI  rank  specific)  input  matrix  and  factor  matrices; \emph{MAD}---element-wise matrix multiplication and division operations; \emph{Norm}---l2 norm computation of local matrices; \emph{INIT}---initialization of factor matrices. The communication cost includes{:} all\textunderscore gather (\emph{AG})---time taken for global matrix-matrix multiplications while distributing the results across all processors; all\textunderscore reduce (\emph{AR})---the time required to compute global Gram matrices, and reduce-scatter (\emph{RSC})---the time to compute global matrix-matrix product using a reduce-scatter operation.
 We also report the breakdown of the scaling performance for distributed reshaping and I/O operations along with NMF operations.

\subsubsection{Strong Scaling}
We use data with fixed size to be $256\times256\times256\times256$ (i.e. 16GB) and vary the number of processors that we use to compute the distributed tensor train decomposition. The times taken for NMF of all factors, for data operations, and the overall TT time for two NMF algorithms, the BCD (block coordinate descent) and the MU (multiplicative update) times, across processor grid of sizes $2^{k} \times 2 \times 2  \times 2$, where $1 \leq k \leq 5$, are presented in Figure~\ref{fig:strong}.

\begin{figure*}
\centering
\begin{subfigure}{.33\textwidth}
  \centering
  \includegraphics[width=\linewidth]{images/Strong_data_overall.png}
  \caption{Strong scaling (Overall)}
\end{subfigure}%
\begin{subfigure}{.33\textwidth}
  \centering
  \includegraphics[width=\linewidth]{images/Strong_data_nmf.png}
  \caption{Strong scaling (NMF)}
\end{subfigure}
\begin{subfigure}{.33\textwidth}
  \centering
  \includegraphics[width=\linewidth]{images/Strong_data_op.png}
  \caption{Strong scaling (Data Operations)}
\end{subfigure}

\caption{Strong scaling experiments }
\label{fig:strong}
\end{figure*}

The strong scaling experiments are run with  tensor-train ranks (TT-ranks) set to be 1,10,10,10, and 1 in the respective dimensions. The scaling results show that the overall TT performance achieves better FLOPS with larger grid size as the overall running time decreases with larger processor size. The scaling saturates at larger core sizes due to the inter-processor communication and the smaller matrix operations within the local computation kernels.

\subsubsection{Weak Scaling}
To test the weak scaling performances of the proposed TT implementation, we report the performance per core for 100 iterations of each TT decomposition stage. Figure~\ref{fig:weak} illustrates the weak scaling performance for the proposed framework.
\begin{figure*}
\centering
\begin{subfigure}{.33\textwidth}
  \centering
  \includegraphics[width=\linewidth]{images/weak_data_overall.png}
  \caption{Weak scaling (Overall)}
\end{subfigure}%
\begin{subfigure}{.33\textwidth}
  \centering
  \includegraphics[width=\linewidth]{images/weak_data_nmf.png}
  \caption{Weak scaling (NMF)}
\end{subfigure}
\begin{subfigure}{.33\textwidth}
  \centering
  \includegraphics[width=\linewidth]{images/weak_data_op.png}
  \caption{Weak scaling (Data Operations)}
\end{subfigure}

\caption{Weak scaling experiments}
\label{fig:weak}
\end{figure*}

\begin{figure*}
\centering
\begin{subfigure}{.33\textwidth}
  \centering
  \includegraphics[width=\linewidth]{images/k_data_overall.png}
  \caption{Scaling with TT-ranks (Overall)}
\end{subfigure}%
\begin{subfigure}{.33\textwidth}
  \centering
  \includegraphics[width=\linewidth]{images/k_data_nmf.png}
  \caption{Scaling with TT-ranks (NMF)}
\end{subfigure}
\begin{subfigure}{.33\textwidth}
  \centering
  \includegraphics[width=\linewidth]{images/k_data_op.png}
  \caption{Scaling with TT-ranks (Data Operations)}
\end{subfigure}

\caption{Scaling experiments with respect to TT-ranks }
\label{fig:features}
\end{figure*}

For this experiment, the size of data is fixed per processor, while we scale up the processor  and the data sizes by the same factors.  Similar to the strong scaling, we use processor grid of sizes $2^{k} \times 2 \times 2  \times 2$, where $1 \leq k \leq 5$. In addition to that, we vary the data size as $256^{k} \times 256 \times 256  \times 256$, where $1 \leq k \leq 5$. The condition for $k=1$ matches with the $k=1$ setup in strong scaling. The size of data varies from 16GB to 256GB for a processor count ranging from 16 to 256. Again, the scaling performance degrades slightly for larger processors and data sizes due to the inter-nodal communications and the I/O involved.
\begin{figure*}%
    \centering
    \begin{subfigure}{.31\textwidth}
    \includegraphics[width=\linewidth]{images/compression_vs_error_vid_svd_nmf_final_loglog.png}
    \caption{ } \label{fig:yale}

    \end{subfigure}
    \qquad
    \begin{subfigure}{.31\textwidth}
    \includegraphics[width=\linewidth]{images/compression_vs_error_vid_svd_nmf_final_loglog.png}
    \caption{ } \label{fig:vid}

    \end{subfigure}
    \begin{subfigure}{.31\textwidth}
    \qquad
    \includegraphics[width=\linewidth]{images/compression_vs_error_512gb_nmf_gr_mu.png}
    \caption{ } \label{fig:synthetic}

    \end{subfigure}
    \caption{Compression ratio vs relative error with tensor train decomposition on datasets a)Yale face b) Video c) Synthetic data (500GB)}%
    \label{fig:TT}%
\end{figure*}
\subsubsection{Scaling with TT ranks}
To demonstrate the scalability of our framework for different tensor train ranks, we fix the total number of processors to be 256 for a tensor of size $256 \times 256 \times 256  \times 256$, and vary the inner tensor train ranks to analyze the effect of $r$ on scaling. Figure \ref{fig:features} shows the scaling with TT ranks results, where each TT rank $r$ is varied in $\{2,4,8,16\}$.

\subsection{Application to real-world dataset}

\subsubsection{Data Description}
\subsubsubsection{Extended Yale Face Dataset B} We first demonstrate the compressibility on Extended Yale Face Dataset B \cite{GeBeKr01,wang2019principal,wang2017efficient} that includes 38 people with 9 poses under 64 illumination conditions. Each image from the Yale Face dataset has size of $192 \times 168$, where each image is down-sampled to $48 \times 42$ for comparison with an existing method \cite{GeBeKr01}. The formatted 4D tensor dimensions are $48 \times 42 \times 64 \times 38$.  We also demonstrate denoising on the same dataset by adding Gaussian noise $N(0,900)$ to each voxel of the tensor.

\begin{figure}
\centering
\begin{subfigure}[b] {0.45\textwidth}
   \includegraphics[width=0.9\linewidth]{images/svd_denoise3.png}
   \caption{}
   \label{fig:Ng1}
\end{subfigure}

\begin{subfigure}[b]{0.45\textwidth}
   \includegraphics[width=0.9\linewidth]{images/nmf_denoise3.png}
   \caption{}
   \label{fig:Ng2}
\end{subfigure}
\caption{Demonstration of image denoising. (a) Denoising with SVD-TT (b) Denoising with NMF-TT } \label{fig:denoise}
\end{figure}

\subsubsubsection{Video}
 The video tensor ($100 \times 260 \times 3 \times 85$) obtained from a high-speed camera video for gun shooting \cite{youtube}, comprises 4 dimensions, where the first two dimensions correspond to a monochromatic image, the third dimension is the channel, and the fourth one is the frame count.

\subsubsection{Compression Ratio and Reconstruction Error}
For a $d$-way data tensor $\ten{A}$ of size $n_{1},n_{2},...,n_{d}$, the reconstructed tensor from the tensor train factors $\ten[i]{G}$ is given by $\ten{\tilde A}$. Then the relative error $\epsilon$ for reconstruction is given as
\begin{equation}
    \epsilon= \frac{\Vert \ten{A} - \ten{\tilde A}\Vert_F}{\Vert \ten{A} \Vert_F}.
\end{equation}

If the tensor train ranks of the decomposition are $r_{0}= r_d = 1,$ $r_{i} \geq 1$ for $1 \leq i < d$, then the compression ratio is measured as
\begin{equation}
    C = \prod_{i=1}^{d} n_i /(\sum_{i=1}^{d} n_{i}*r_{i-1}*r_{i})
\end{equation}

Each data point in the Figures~ \ref{fig:yale} and \ref{fig:vid} is the compression ratio $C$ at the desired relative error $\epsilon$ at each TT decomposition stage. The targeted relative errors at each TT stage for selection of $r_i$ was set to 0.5, 0.25, 0.125, 0.075, 0.01, 0.005 and 0.001. Lower TT rank produces higher compression and higher reconstruction error, whereas higher TT rank produces lower compression and better reconstruction. For the Yale face dataset, the compression ratio $C$ varies from 1.13 with $\epsilon= 0.04$ to $C$ of 2.55e4 for $\epsilon= 0.55$ with  nTT. Similarly with TT, the compression ratio $C$ varies from 1.32 with $\epsilon= 0.013$ to $C$ of 2.55e4 for $\epsilon= 0.55$. Also, for the video dataset, the compression ratio $C$ varies from 1.01 with $\epsilon= 0.015$ to $C$ of 1.47e4 for $\epsilon= 0.54$ with  nTT. Similarly with TT, the compression ratio $C$ varies from 1.007 with $\epsilon= 0.0004$ to $C$ of 1.47e4 for $\epsilon=  0.53$.

\subsubsection{Application to image denoising}

To demonstrate the efficacy of non-negative tensor train over the regular tensor train, we apply both techniques for the decomposition of a noisy Yale Face sample and report the denoising performance. The metric that we choose to evaluate the correctness of reconstruction compared to the noise-free sample is the structural similarity (SSIM) index~\cite{clifton1938osteopetrosis}. SSIM is a widely used metric for image similarity measures in computer vision applications. SSIM ranges on a scale of [0,1], where 1 is the best match.  In Figure~\ref{fig:denoise}, we report the SSIM values of the reconstructed images with respect to the noiseless ground truth image. Figure~\ref{fig:Ng1} corresponds to the regular tensor train based reconstruction, whereas Figure~\ref{fig:Ng2} corresponds to the non-negative TT based reconstruction. For both of these figures, the images from left to right correspond to the reconstructed images with decreasing TT-ranks and increasing compression rates. The value of the top of each image corresponds to the SSIM measure with respect to the ground truth and the value on the bottom correspond to the compression rate. We can observe that increasing the compression-rate-based decomposition eliminates noise significantly and results in the reconstructed image to resemble the original noise-free image. For regular TT/SVD-TT based reconstruction, the best SSIM reported for the reconstructed image is 0.85 whereas with the non-negative TT/nTT based reconstruction, the best SSIM reported is 0.88. For given TT ranks, the reconstructed image SSIM for nTT is better than that for the TT.

\subsubsection{Compression of large synthetic data}

We demonstrate the compression ratios with lower reconstruction error for a 500 Gigabyte (GB) matrix with dimensions $1024 \times 512 \times 512 \times 512$ and tensor train ranks =[1,20,30,40,1]. We synthetically generate the data in a distributed manner as discussed in the data generation section. Figure~\ref{fig:synthetic} shows compression ratios with two different NMF optimization methods, BCD vs multiplicative update, both of which are based on minimization of the Frobenius norm. For the BCD optimization-based NTF, the compression ratio $C$ varies from 163880 with $\epsilon=3e-6$ to $C$ of 6391320 for $\epsilon= 0.03$ with non-negative tensor train. Similarly, with multiplicative update algorithm based NTF, the relative error $\epsilon$  varies from 0.005 to 0.045 for the same compression range. This experiment demonstrates the tradeoff between two different NMF update algorithms BCD and MU for nTT. As per the scaling plots \ref{fig:strong},\ref{fig:weak} and \ref{fig:features}, MU algorithm demonstrates better timings. However, BCD achieves a better compression rate with lower reconstruction error compared to MU as per Figure \ref{fig:synthetic}.

\section{Conclusion and Future directions}
Here, we introduce a distributed non-negative tensor train, nTT, algorithm that is capable of computing tensor train to a prescribed relative error. We demonstrate the nTT scaling performance on synthetic data, and establish its ability to decompose a 500GB tensor. Finally, we apply the algorithm to various real and synthetic datasets to demonstrate the nTT data compression capabilities. In the future, we aim to apply our framework to large real-world datasets such as seismic datasets, satellite images, medical images, etc., in order to have an efficient highly compressed representation and to be able to do inference/classification from the nTT low-dimensional representations.

\section{Acknowledgements}
This research was funded by Laboratory Directed Research and Development  (20190020DR), and resources were provided by the Los Alamos National Laboratory Institutional Computing Program, supported by the U.S. Department of Energy National Nuclear Security Administration under Contract No. 89233218CNA000001.

\bibliographystyle{IEEEtran}