\chapter{Implementation and Complexity}\label{chapter4}

In this chapter, we will first describe how the input data should be stored and accessed along with practical improvements of the formal pseudocode that will improve its efficiency. Next, we will calculate the space complexity of the presented implementation and time complexity of one iteration. Finally, we present the actual interface of the implemented minimization algorithm in C++ and the structure of the input and output data files.

\section{Implementation Details}

Here, we give the implementation details of the previously presented pseudocode and used data structures.

\subsection{Storing Input Data}

Concerning the input data, it is straightforward how to store the vectors $\mathbf{m}$ and $\mathbf{b}$ but it is necessary to deal with efficiently storing the vectors $\mathbf{a}$. As it was presented in Algorithm \ref{transformationalg_summax}, the vectors are likely to be sparse. Therefore, we should store each one as a list of non-zero $a_{i,j,k}$ values along with the coordinates $k$ that these values correspond to. Observe that this kind of storage automatically allows efficient querying of $N(f_{i,j})$ elements and the corresponding coefficients $a$, which is frequently done.

Similarly, we could store for each coordinate $k$ the non-zero coefficients $a_{i,j,k}$ in a list along with the $(i,j)$ indices. This structure will be used for querying $N(x_k)$.

Later in this section, we will show that it is never necessary to access the values $a_{i,j,k}$ randomly, but always in the scope of iterating over $N(x_k)$ or $N(f_{i,j})$.

\subsection{Efficiency Improvements}

First of all, it is important to note that on each place in the algorithm where the value of a subfunction is queried in the current point, we should not calculate it every time but use the pre-defined vector $\mathbf{y}$ that satisfies $y_{i,j} = f_{i,j}(\mathbf{x})$ and is updated after each iteration. In addition, we keep a vector $\bm{\tau} \in \mathbb{R}^l$ such that
\begin{equation}
\tau_i = \max_{j = 1}^{m_i} f_{i,j}(\mathbf{x}),
\end{equation}
so that we can also replace the occurrences of cluster maxima and not evaluate them each time.

\subsubsection{Local Consistency Algorithm}

In the local consistency algorithm, we introduce some auxiliary variables for efficiency. One of these vectors, $\bm{\varrho} \in \mathbb{N}^l$, contains in each $\varrho_i$ the amount of remaining subfunctions that are alive in cluster $i$,
\begin{equation}
\varrho_i = \sum_{j = 1}^{m_i}\iverson{p(i,j) = \mathsf{ALIVE}}.
\end{equation}
The initial values of this vector are zeros but whenever a subfunction $f_{i,j}$ becomes alive, we need to increase the corresponding component and do the opposite when it becomes inactive or killed. The same condition that is checked on lines 2 and 13 in Algorithm~\ref{main_ac_alg} is equivalent to checking whether some $\varrho_i=0$, where $i$ are the clusters of the subfunctions that were just killed. Updating the vector does not create any overhead in the sense of asymptotic time complexity because in each assignment to the vector $\mathbf{p}$, we will update at most one value in $\bm{\varrho}$. 

Next, to see whether a coordinate is consistent, we do not need to iterate over all the subfunctions but instead, we can use vectors $\bm{\gamma}^+,\bm{\gamma}^- \in \mathbb{N}^n$, where $\gamma^+_k$ (resp. $\gamma^-_k$) is the amount of subfunctions $f_{i,j}$ that are alive and their corresponding coefficient $a_{i,j,k}$ is positive (resp. negative). Formally,
\begin{subequations}
\begin{align}
\gamma_k^+ & = \sum_{i=1}^l\sum_{j=1}^{m_i} \iverson{(a_{i,j,k} > 0) \wedge (p(i,j) = \mathsf{ALIVE})}\\
\gamma_k^- & = \sum_{i=1}^l\sum_{j=1}^{m_i} \iverson{(a_{i,j,k} < 0) \wedge (p(i,j) = \mathsf{ALIVE})}.
\end{align}
\end{subequations}
Testing for local consistency of a single variable $x_k$ is then equivalent to checking
\begin{equation}
(\gamma_k^+ = 0) \iff (\gamma_k^- = 0).
\end{equation}
Similarly as with $\bm{\varrho}$, the values of $\bm{\gamma}^+$ and $\bm{\gamma}^-$ need to be updated whenever a subfunction becomes alive, i.e. evaluate
\begin{subequations}
\begin{align}
\gamma^+_k \gets \gamma^+_k + 1, \quad \forall k \in N(f_{i,j}), a_{i,j,k} > 0 \\
\gamma^-_k \gets \gamma^+_k + 1, \quad \forall k \in N(f_{i,j}), a_{i,j,k} < 0 
\end{align}
\end{subequations}
and whenever a subfunction is killed or becomes inactive, it is necessary to subtract $1$ from the corresponding values.

It is also important to mention that the queue $Q$ that is used in the local consistency algorithm should be implemented so that it stores each value $k$ at most once -- i.e. if the function $Q$.push($k$) is called and the queue already contains $k$, then no action should be executed.

The final remark to the local consistency algorithm is that we should have a vector $\mathbf{\bar{a}}$ with values $\bar{a}_{i,j} \in \mathbb{R}\cup \{\mathsf{UNDEFINED}\}$ for each subfunction $f_{i,j}$. Initially, all $\bar{a}_{i,j} = \mathsf{UNDEFINED}$, but whenever some $p(i,j)$ value is set to a coordinate from $[n]$, then we set $\bar{a}_{i,j} \gets a_{i,j,p(i,j)}$. The benefit of this approach is that we can eliminate all randomly accessed values $a_{i,j,k}$ from the algorithm. In the pseudocode, all these values are used only while iterating through the corresponding sets $N(x_k)$ or $N(f_{i,j})$, except for some occurrences of $a_{i,j,p(i,j)}$, which can be replaced by $\bar{a}_{i,j}$. Therefore, we do not require anywhere to access a random $a_{i,j,k}$ value but only obtain these values while querying the sets $N$.

\subsubsection{Decreasing Direction, Line Search, and Update}

When Algorithm \ref{directionalg_dag} calculates the direction $\Delta \mathbf{x}$, the DAG or the augmenting DAG should not be explicitly created, except for storing which subfunctions are in the augmenting DAG (which is stored in $\mathbf{r}$) and their in-degrees (which are stored in $\mathbf{d}$). The existence of particular edges is then determined implicitly by their definition.

During the creation of $\Delta \mathbf{x}$ vector, whenever a value is assigned to some $\Delta x_k$, we push the coordinate $k$ into a queue $Q_x$, which is initially empty. And after that, when the set $N(x_k)$ is traversed, we put all the clusters $i'$ into a queue $Q_i$, which was also empty at the start of direction calculation algorithm. These queues should be implemented in the same manner as the queue in the local consistency algorithm, i.e. they will contain each value at most once. Eventually, the queue $Q_x$ will contain all the non-zero coordinates of $\Delta \mathbf{x}$ and $Q_i$ will contain all the clusters $i$ such that this cluster contains at least one subfunction $f_{i,j}$ that depends on a non-zero coordinate of $\Delta \mathbf{x}$.

Now, when the step size is calculated, we can replace the occurrences of $f_{i,j}(\mathbf{x})$ with $y_{i,j}$ and do the same with the cluster maxima, which can be replaced by $\bm{\tau}$. Additionally, we can read the cluster indices $i$ from the queue $Q_i$ to find all potentially increasing subfunctions, as it is done in the for-loops of the line search defined in Algorithm \ref{linesearch_dag}.

After that, in the update procedure, we iterate over the clusters stored in $Q_i$ again to update the corresponding values in $\mathbf{y}$ and $\bm{\tau}$ and during that, we empty the queue. Finally, the current point $\mathbf{x}$ is updated but we only need to recalculate the coordinates that are stored in the queue $Q_x$, which is emptied during this process.

\section{Time and Space Complexity}

In this section, we will discuss the asymptotic space complexity of the algorithm and the asymptotic time complexity of one iteration. We will assume that the vectors $\mathbf{a}_{i,j}$ are sparse, which formally means that the size of sets $N(f_{i,j})$, respectively $N(x_k)$ does not exceed a small constant $K$, i.e.
\begin{equation}
K \geq \max \left( \max_{i = 1}^l{\max_{j = 1}^{m_i}{\abs{N(f_{i,j})}}} , \, \max_{ k = 1}^n{\abs{N(x_k)}} \right),
\end{equation}
thus $\abs{N(f_{i,j})} \in O(1)$ for all subfunctions $f_{i,j}$ and $\abs{N(x_k)} \in O(1)$ for all $k \in [n]$. We will come to the result that the asymptotic time complexity of each iteration is $O(n+M)$, which is also the asymptotic space complexity of the algorithm. We denote the total number of subfunctions in the SMAF as $M$.

\subsection{Space Complexity}

Among all the variables that the algorithm uses, there are some whose size is proportional to the total number of subfunctions. These are for example $\mathbf{y}$, $\mathbf{c}$, $\mathbf{d}$, $\mathbf{r}$, $\mathbf{p}$ and $Q_f$. There are also structures whose size is proportional to the amount of variables $n$, e.g. $\mathbf{x}$, $\Delta \mathbf{x}$ and $Q$. There are no structures whose size would be asymptotically larger than $n$ or $M$.

The input is also stored as a sparse matrix, so the total space complexity of it is not higher than $O(nK+MK) = O(n+M)$. Thus, the overall asymptotic space complexity is $O(n+M)$.

\subsection{Time Complexity}

\subsubsection{Local Consistency Algorithm}

As mentioned previously in this chapter, the first condition in Algorithm \ref{main_ac_alg} can be checked in $O(l)$ using the precalculated values $\bm{\varrho}$.

Next, we would like to estimate an upper bound on the number of loops driven by the queue $Q$. Assume that the queue is initialized with a non-zero amount of variables, then it could happen that even before all these variables are processed, the termination condition is satisfied in less than $n$ loops. However, if all the initial elements of the queue are processed and the algorithm did not end yet, we can say that we have killed $v_1$ subfunctions until now. That means that there is at most some $K\cdot v_1$ variables that were added while the subfunctions were being killed. Then, we can use the same reasoning and see whether the algorithm ends before all these variables are processed. If it does not, we can denote $v_2$ as the number of subfunctions that were killed due to the previously added variables.

\begin{figure}
\begin{tikzpicture}
\node at (-0.5,0.35) {$Q$:};
\draw (5,0) -- (0,0) -- (0,0.7) -- (5,0.7) -- (5,0);
\node at (2.5,0.35) {\small initial content};
\draw [->] (2.5,-0.2) -- (2.5,-1.5); 
\node at (2.5,-2) {\small killed $v_1$ subfcns};
\draw [->] (3.7,-1.5) -- (5.7,-0.2); 
\draw (8.1,0) -- (5.1,0) -- (5.1,0.7) -- (8.1,0.7) -- (8.1,0);
\node at (6.55,0.30) {\small added by $v_1$};
\draw [->] (6.5,-0.2) -- (6.5,-1.5); 
\node at (6.5,-2) {\small killed $v_2$ subfcns};
\draw [->] (7.2,-1.5) -- (8.7,-0.2); 
\draw (11.2,0) -- (8.2,0) -- (8.2,0.7) -- (11.2,0.7) -- (11.2,0);
\node at (9.6,0.30) {\small added by $v_2$};
\draw [->] (9.6,-0.2) -- (9.6,-1.5); 
\node at (9.6,-2) {\small killed $v_3$ subfcns};
\draw [->] (10.3,-1.5) -- (11.8,-0.2); 
\draw (12.3,0) -- (11.3,0) -- (11.3,0.7) -- (12.3,0.7);
\node at (12.3,0.35) {\small ...};
\node at (12.3,-2) {\small ...};
%\draw[ ultra thick, -{Latex[length=6mm,width=4mm]}] (11,-6) -- (13,-6); - tady se nastavuje jak délka, tak šířka toho hrotu

\end{tikzpicture}
\caption{Processing coordinates in queue.}
\label{processing_coordinates_fig}
\end{figure}


This idea is shown visually in Figure \ref{processing_coordinates_fig}. In the figure, there is some initial content of the queue during which $v_1$ subfunctions are killed. Because these subfunctions are killed, they add at most $K \cdot v_1$ variables to the queue that kill $v_2$ subfunctions etc. But it is important that the algorithm ends at some point. Necessarily, the total amount of killed subfunctions until this point is lower than the total amount of subfunctions, i.e.
\begin{equation}
\sum_i v_i \leq M.
\end{equation}
Because each killed subfunction added at most constant amount of variables into the queue, the total amount of added variables during the run of the algorithm is $O(KM)$. Assuming that there could have initially been at most $n$ variables in the queue, there are at most $O(KM+n) = O(M+n)$ iterations of the while loop.

Each iteration of the while loop takes constant time because checking consistency is easy after introducing vectors $\bm{\gamma}^+$ and $\bm{\gamma}^-$, there is only constant amount of subfunctions that depend on the inconsistent variable and each of these subfunctions depends on constant amount of variables. The checking of the terminating condition can be done for the clusters of the subfunctions that are just being killed.

Therefore, the local consistency algorithm has asymptotic time complexity $O(n+M)$.


\subsubsection{Decreasing Direction Calculation}

In the calculation of $\Delta \mathbf{x}$ direction, the augmenting DAG is traversed in order to find the degrees of its nodes. Obviously, there can be at most $M$ nodes in the DAG, therefore the maximum number of edges is $O(KM) = O(M)$, which means that the first traversal is at most linear w.r.t. the number of subfunctions.

Then, we go through the augmenting DAG again to calculate the direction. Because the size of augmenting DAG is at most $O(M)$ and each subfunction is processed in $O(K) = O(1)$ time, we can conclude that the whole process of calculating $\Delta \mathbf{x}$ has time complexity $O(M)$.

\subsubsection{Line Search and Update}

During the line search, we iterate over the subfunctions in the cluster $i^*$ at most twice and then over all clusters that may contain a subfunction with non-zero slope. This procedure can be performed in $O(M)$ time because in the worst case, we would iterate over all the subfunctions.

The same can be said about the update -- we might have to update all the subfunctions (and also all the values of the maxima) and all components of $\mathbf{x}$ and then iterate over the updated values to find out which subfunctions are active for the next iteration. This can be done in $O(2M+n) = O(M+n)$ time. It is important to note that the deactivate function runs in constant time $O(K)$ and the accumulated runtime of activate function will not exceed $O(M)$ because each subfunction's status could be set to $\mathsf{ALIVE}$ at most once and recursive resurrecting will not run more than once for each subfunction.

To conclude this section, we have estimated the time complexity of one iteration of the algorithm to be $O(n+M)$, where $n$ is the number of variables and $M$ is the total number of subfunctions in the SMAF.

\subsubsection{Complexity of the Whole Algorithm}

In Section \ref{finite_num_of_iter_section}, we introduced an upper bound that showed that the maximum number of the iterations could theoretically be exponential in the size of the input. It was out of scope of this thesis to find a better bound on its convergence rate and we do not know whether a better bound exists. However, as said above, the number of iterations was experimentally shown to be much lower than the exponential bound.


\section{Interface}

In this section, we describe the input format of the file that defines a SMAF and also the output format of the file that returns the result of the algorithm. We also briefly comment on our implementation in C++, namely its division into classes and description of their functionalities.

\subsection{Input and Output File Formats} \label{io_section}

This section describes the way the input and output files are structured to allow the user to interpret them in own custom applications.

\subsubsection{Input Format}

The input format was designed so that it would be easy and efficient to load the data to an internal representation. The input file is a text file with content according to Figure \ref{inputfileformat}.
\begin{figure}
\begin{tcolorbox}
$l$ $\,$ $n$ $\,$ $K$

$m_1$ $\,$  $m_2$  $\,$ ... $\,$  $m_l$ 
\iffalse
\begin{subequations}
\begin{align*}
e_{1,1}  & \quad k_1^{1,1} & \quad  a_{1,1,k_1^{1,1}} & \quad  k_2^{1,1} & \quad  a_{1,1,k_2^{1,1}}  & \quad \cdots  & \quad k_{e_{1,1}}^{1,1} & \quad  a_{1,1,k_{e_{1,1}}^{1,1}} & \quad  b_{1,1}\\
e_{1,2}  & \quad k_1^{1,2} & \quad  a_{1,2,k_1^{1,2}}  & \quad k_2^{1,2}  & \quad a_{1,2,k_2^{1,2}}  & \quad \cdots  & \quad k_{e_{1,2}}^{1,2}  & \quad a_{1,2,k_{e_{1,2}}^{1,2}}  & \quad b_{1,2}\\
\vdots \\
e_{1,m_1}  & \quad k_1^{1,m_1} & \quad  a_{1,m_1,k_1^{1,m_1}}  & \quad k_2^{1,m_1}  & \quad a_{1,m_1,k_2^{1,m_1}}  & \quad \cdots  & \quad k_{e_{1,m_1}}^{1,m_1}  & \quad a_{1,m_1,k_{e_{1,m_1}}^{1,m_1}}  & \quad b_{1,m_1} \\
e_{2,1}  & \quad k_1^{2,1} & \quad  a_{2,1,k_1^{2,1}}  & \quad k_2^{2,1}  & \quad a_{2,1,k_2^{2,1}}  & \quad \cdots & \quad k_{e_{2,1}}^{2,1}  & \quad a_{2,1,k_{e_{2,1}}^{2,1}}  & \quad b_{2,1} \\
\vdots \\
e_{l,m_l}  & \quad k_1^{l,m_l} & \quad  a_{l,m_l,k_1^{l,m_l}}  & \quad k_2^{l,m_l}  & \quad a_{l,m_l,k_2^{l,m_l}}  & \quad \cdots  & \quad k_{e_{l,m_l}}^{l,m_l}  & \quad a_{l,m_l,k_{e_{l,m_l}}^{l,m_l}}  & \quad b_{l,m_l}
\end{align*}
\end{subequations}
\fi

$e_{1,1}$  $\,$ $k_1^{1,1}$ $\,$  $a_{1,1,k_1^{1,1}}$ $\,$  $k_2^{1,1}$ $\,$  $a_{1,1,k_2^{1,1}}$  $\,$ ...  $\,$ $k_{e_{1,1}}^{1,1}$ $\,$  $a_{1,1,k_{e_{1,1}}^{1,1}}$ $\,$  $b_{1,1}$

$e_{1,2}$  $\,$ $k_1^{1,2}$ $\,$  $a_{1,2,k_1^{1,2}}$  $\,$ $k_2^{1,2}$  $\,$ $a_{1,2,k_2^{1,2}}$  $\,$ ...  $\,$ $k_{e_{1,2}}^{1,2}$  $\,$ $a_{1,2,k_{e_{1,2}}^{1,2}}$  $\,$ $b_{1,2}$

$\vdots$

$e_{1,m_1}$  $\,$ $k_1^{1,m_1}$ $\,$  $a_{1,m_1,k_1^{1,m_1}}$  $\,$ $k_2^{1,m_1}$  $\,$ $a_{1,m_1,k_2^{1,m_1}}$  $\,$ ...  $\,$ $k_{e_{1,m_1}}^{1,m_1}$  $\,$ $a_{1,m_1,k_{e_{1,m_1}}^{1,m_1}}$  $\,$ $b_{1,m_1}$

$e_{2,1}$  $\,$ $k_1^{2,1}$ $\,$  $a_{2,1,k_1^{2,1}}$  $\,$ $k_2^{2,1}$  $\,$ $a_{2,1,k_2^{2,1}}$  $\,$ ...  $\,$ $k_{e_{2,1}}^{2,1}$  $\,$ $a_{2,1,k_{e_{2,1}}^{2,1}}$  $\,$ $b_{2,1}$

$\vdots$

$e_{l,m_l}$  $\,$ $k_1^{l,m_l}$ $\,$  $a_{l,m_l,k_1^{l,m_l}}$  $\,$ $k_2^{l,m_l}$  $\,$ $a_{l,m_l,k_2^{l,m_l}}$  $\,$ ...  $\,$ $k_{e_{l,m_l}}^{l,m_l}$  $\,$ $a_{l,m_l,k_{e_{l,m_l}}^{l,m_l}}$  $\,$ $b_{l,m_l}$

\end{tcolorbox}
\caption{Input file format.}
\label{inputfileformat}
\end{figure}

On the first line, there should be the number of clusters $l$, number of variables $n$ and a constant $K$ that is an upper bound on the number of non-zero elements in any $\mathbf{a}_{i,j}$ and also the number of subfunctions that depend on a variable $x_k$.

The second line contains the sizes of clusters, i.e. how many subfunctions there are in each cluster.

Next, there are $M = \sum_{i \in [l]}m_i$ lines which define the individual subfunctions. The order of the lines is straightforward -- we first define the $m_1$ subfunctions of the first cluster, then the $m_2$ subfunctions of the second one etc. 

When defining a subfunction $f_{i,j}$, the first number on the corresponding line is $e_{i,j}$, which equals to $\abs{N(f_{i,j})}$, i.e. the number of variables that this subfunction depends on. This number is followed by exactly $e_{i,j}$ pairs of values, where the first number of the pair is a coordinate $k$ and the second number is the corresponding $a_{i,j,k}$ value. Obviously, only coordinates with non-zero  $a_{i,j,k}$ values are listed and after all of them are indexed, the line is ended with the value of $b_{i,j}$. Notice that in the input files, coordinates are not numbered from 1 to $n$ but instead from 0 to $n-1$.

\subsubsection{Output Format}

The output format is much simpler when compared to the input format. On the first line, it contains the input number of clusters $l$, number of variables $n$, and $\epsilon$ such that the returned point $\mathbf{x}$ is locally $\epsilon$-consistent. On the second line, the values of $\mathbf{x}$ are listed. On the third line, there are $l$ numbers $h_1$, ..., $h_l$ that contain the index $j$ of the subfunction in the corresponding cluster which has status $\mathsf{ALIVE}$. The values of $\mathbf{p}$ are taken from the result of the local consistency algorithm which was run on $\epsilon$-active subfunctions. If there are multiple alive subfunctions in a cluster, $h_i = -1$. Formally, we define
\begin{equation}
h_i = \left\{
        \begin{array}{ll}
            \sum\limits_{j \in [m_i]} \, j \cdot \iverson{p(f_{i,j}) = \mathsf{ALIVE}} & \text{if there is only one $j$ with $p(f_{i,j}) = \mathsf{ALIVE}$} \\
            -1 & \text{otherwise}
        \end{array}.
    \right.
\end{equation}

Similarly as with the variable indices in the input, the subfunctions are numbered from 0 to $m_i-1$ instead of 1 to $m_i$ in the actual input.

\begin{figure}

\begin{tcolorbox}
$l$ $n$ $\epsilon$

$x_1$ $x_2$ $\cdots$ $x_n$

$h_1$ $h_2$ $\cdots$ $h_l$

\end{tcolorbox}
\caption{Output file format.}
\label{outputfileformat}
\end{figure}

\subsection{The Code and its Usage}

The implementation is divided into 3 classes. Class $\mathtt{Instance}$ stores the input instance of SMAF and provides the data to the minimization algorithm, this class is initialized using a path to a source file in the format shown in Figure \ref{inputfileformat}.

The minimization algorithm is implemented in class $\mathtt{Solver}$, which is initialized by an object of class $\mathtt{Instance}$. The minimization process then can be started by invoking the $\mathtt{solve()}$ function, which accepts the initial $\epsilon$ as its optional parameter -- if it is not specified, then it uses the initialization shown on line 3 in Algorithm~\ref{decreasing_epsilon_integer_alg}. By default, the program writes the information concerning every iteration on standard output but this can be turned off by setting the verbose flag to false.

After the $\mathtt{solve()}$ function terminates, the user can call the $\mathtt{result()}$ function that reports the reached function value and the minimum $\epsilon$ value for which the current point is locally $\epsilon$-consistent. It also uses the theorems presented in Section~\ref{opt_criteria_section} to decide the optimality of the result (if it is possible). The function $\mathtt{result()}$ also accepts a path to a text file as its optional parameter. If the path is specified, then the output file in the form shown in Figure \ref{outputfileformat} is created and the reached point $\mathbf{x}$ is stored there. It is also possible to change the values of vector $\mathbf{b}$ by invoking the function $\mathtt{changeB(i, j, d)}$ that sets $b_{i,j} \gets b_{i,j} + d$. After this change is done, the result can be re-optimized by the $\mathtt{solve()}$ function that would start from the previously found point.

The third class is $\mathtt{Queue}$, which is an implementation of the queue that contains only unique elements and is used in multiple places throughout the algorithm. Except for the classes, we also attach an example code $\mathtt{main.cpp}$ that shows the usage of our implementation.