A Benchmark for Maximum-A-Posteriori Inference Algorithms in Discrete Sum-Product Networks

Heitor Reis Ribeiro. Algorithms for Maximum-a-Posteriori inference in discrete SumProduct Networks. 2021. 56 f. Tese (Mestrado) Instituto de Matemática e Estatística, Universidade de São Paulo, São Paulo, 2021. The solution to Maximum-a-Posteriori Inference problems in Sum-Product Networks provides the most probable configuration of the Random Variables encoded in its structure; a key step in Probabilistic reasoning that can be used for many applications, such as image auto-completion. It has been proven that this problem is NP-Hard (even to approximate) in Sum-Product Networks. Multiple algorithms have been developed to reach either approximate or exact solutions to this problem, but the experiments have been limited. In this Dissertation we provide descriptions, analysis, and a benchmark for experimental testing for algorithms that solve this problem. We conclude that, given limited time, a Local Search algorithm starting with a solution found by the Argmax-Product algorithm reaches, on average, better results on the tested datasets.

What are Sum-Product Networks?
Sum-Product Networks (SPNs) are probabilistic models introduced in [PD11] that represent a joint probability distribution over a set of random variables as a rooted, directed graph of arithmetic operations (Sums and Products). Sum operations are represented with sum nodes ( + ), and they are connected to other nodes with a weighted arc. Product operations are represented with product nodes ( × ) and they are connected to other nodes with unweighted arcs. The leaves in these graphs are indicator nodes ( • ) and they represent one particular assignment to a random variable. SPNs have been used in problems such as activity recognition [AT15], language modeling [CKP + 14], and image classification [SC16]. Figure 1.1 shows an example of a small SPN over two variables, X and Y , both binary, where x represents an indicator of X = 1, and x, X = 0. The same for variable Y .
Similar to the more known Probabilistic Graphical Models (PGMs), SPNs are used as a way of reasoning under uncertainty, that is, using data with an unknown probabilistic distribution that can be inferred with a large enough data sample, such as a collected or generated dataset, in order to make future decisions (or predictions) over it. For example, given a dataset of census data (age, education, occupation, marital status, income, etc.), a probabilistic graphical model can predict, with measurable uncertainty, the value of (at least) one of the obtained variables (such as income) given the value of some of the other variables -a process called conditional inference. The main advantage of SPNs is the calculation of exact, marginal and conditional inferences in linear time in relation to their size -operations that are NP-Hard on the other PGMs.

The Maximum-A-Posteriori Problem
A key component in artificial intelligence problems that use graphical models to reason under uncertainty is finding the maximum probability configuration, that is, an attribution (or assignment) to a set of variables that yields, upon inference, the highest probability value. One can use such an operation to answer questions such as "What is the most probable income for a person older than 35 that lives in a specific city?" or "What is the most probably entire image given half of its pixels?", given enough relevant data to be analyzed. It has been used previously in image segmentation [GG84], 3D image reconstruction [BVZ98], natural language processing [KRC + 10], and sentiment analysis [ZNSS11].
The operation of finding this configuration on a PGM is named Maximum-a-Posteriori (MAP) Inference. In SPNs with minimum height of 3, this type of inference is NP-Hard [CMdC17] even for finding an approximation of the value obtained by the configuration. That is, under the standard assumption of complexity theory that P = N P , one cannot efficiently find such a configuration for an arbitrary SPN. 1 There have been few studies that examine this problem. Besides the previously mentioned article, where the problem's complexity was proven, the MAP Inference problem has been described in more detail in [Peh15]. Up until then, we had one algorithm for this problem, the Max-Product (first described in [PD11]), which is fast (it takes linear time on the size of the SPN), but has no guarantee about the solution it finds. In [CMdC17] they introduced the ArgMax-Product algorithm, an extension on Max-Product that takes quadratic time on the size of the SPN, and its results are never worse than the results found by Max-Product.
Then in [MJT18] the authors proposed and evaluated two new approximate algorithms: Beam Search, an adaptation of hill-climbing to SPNs that uses derivatives to find new solutions, and K-Best Tree, an extension of Max-Product that finds K solutions for each node and tries to evaluate the best one. They also introduced the first (and, for now, only) exact algorithm: Max Search, which is slow but can use different heuristics to achieve a faster running time, and can be stopped at anytime to obtain a solution that improves with more time given. However, their evaluations had limitations: They were restricted to binary SPNs, and they only checked running time and in how many runs the algorithms reached the highest value within the allocated time. They did not test the quality of the found solutions, and had only a couple of heuristics for Max Search. 2 Recently, [MRKA20] proposed a new approximate algorithm, the Hybrid-Product Belief Propagation, and a new exact algorithm based on a MILP reformulation of the MAP problem, but they performed similar or worse than previous algorithms. An up-to-date survey about SPNs was also recently given in [PSCD20], with some discussion about the MAP problem, but no new algorithms were developed.

Proposed Setup and Outcomes
One approach of dealing with the MAP problem is using the approximate algorithms, since they find a result that is "good enough": not a measurable approximation, since such a thing would also be NP-Hard, but that is usually 3 not far from the optimal solution. These algorithms may even find the optimal solution in some cases, but are not reliable enough to do so consistently. Still, their running-time is considerably lower than the algorithms that find the optimal solution (both the trivial greedy-search and Max Search), and thus they may have utility.
In this dissertation we examine, experiment, and provide empirical results for algorithms designed to find a solution to the Maximum-a-Posteriori inference problem in Sum-Product Networks (SPNs) over discrete random variables. We created a benchmark of MAP Inference problems for SPNs with a sizeable collection of representative instances, and evaluated the performance of existing algorithms in terms of speed and quality of the found solution, under different metrics. In particular, we measure the solutions found by approximate algorithms and compare them against the fastest algorithm available, Max-Product, both in time and in the values found using many datasets, including non-binary ones.
We also test different applications of these algorithms, such as applying the found solutions to datasets with missing instances, testing classification models on auto-completed instances, and testing different heuristic choices for the Max Search algorithm that were not included in the algorithm's original article.
The main questions we seek to answer are: 1. Given a limited time-frame, which algorithms perform best in comparison with Max-Product?
2. Which heuristics in Max Search find good solutions faster?
3. Does a larger value for a MAP solution provide any gain in a downstream task (such as classification) that depends on that solution?
Our experiments show that, given a limited time, the approximate algorithms most often outperform the value found by Max Search, and thus they may be preferred for some applications. Also, the Max-Product algorithm is the quickest way of finding an initial solution that can be improved upon given more time allocated; and the ArgMax-Product algorithm introduced in [CMdC17] is a slower way of finding solutions that are as good as or better than the Max-Product's solution.
For the algorithms that need an initial solution to improve upon, such as the Max Search (exact solution) or Local Search (approximate solution), the experiments show that Local Search outperforms Max Search on larger SPNs if given the Max-Product solution to improve upon, and on most other SPNs if given the ArgMax-Product solution.
Our main contributions are: • An SPN library containing implemented algorithms for learning SPNs from data, and executing the MAP algorithms. It is available at https://gitlab.com/marcheing/spn-rs.
• A comparative study of the running-time and MAP values obtained by the existing algorithms.
• An examination of the heuristics for variable and Upper Bound selection on Max Search.

1.3
Chapter 2 Background: introduction to SPNs and MAP inference

Sum-Product Networks
We assume the reader is familiar with probability distributions and (discrete) random variables. We represent random variables by capitalized letters such as X and Y , a set of random variables by bold capitalzed letters such as X = {X 1 , X 2 , . . . , X n }, the support of a variable X as val(X) = {x 1 , x 2 , . . . , x n }, and the set of supports of a set of random variables X as val(X) = × n j=1 val(X j ) (the cartesian product of all the val(X j ) sets). The size of a set X is represented by |X|≥ 0. A partial assignment to a single variable X is defined as X ⊆ val(X), and if there is only one value for the variable X in X , we call it an instantiation of X. A partial assignment to a set of variables X is defined by X = × n j=1 X j , that is, X ⊆ val(X). Given the variables X, if a partial assignment to this set of variables X contains a single assignment to every variable X ∈ X (that is, one value for each variable), we call X a complete assignment and represent it with a lowercase bold letter x. We define the operation of obtaining the values assigned to a variable X on a partial assignment X with X [X], the result being a partial assignment X to a single variable X. We represent the operation of removing the values assigned to a variable X from a partial assignment to a set of variables X with X \ {X} = X , where X is now a partial assignment to a set of variables excluding X. We can combine partial assignments to disjunct sets of variables X and Y with the operation X × Y.
With these notations, we can define Sum-Product Networks properly.
Definition: Sum-Product Networks (SPNs) 1 An SPN over a set of random variables X is a rooted, directed acyclic graph (DAG) containing three types of nodes: • Sum Nodes + .
• Leaf/Indicator Nodes • , each one associated with a variable X ∈ X and a value x ∈ val(X).
• Whenever there is an arc from node A to node B, we say that A is a parent of B, and B is a child of A.
• The set of all child nodes of a node A will be represented by ch(A). The set of all parent nodes of a node B will be represented by pa(B). With an associated scope function: Definition: Scope of a Node • The scope of a leaf/indicator node associated with a Variable X is {X}.
• The scope of a sum or a product node is the union of the scopes of their child nodes.
And following these restrictions: • Each sum node is connected to its children by weighted arcs with non-negative weights.
• Each child of a sum node has the same scope as the other children of the same sum node.
• Each product node is connected to its children by unweighted arcs.
• Each child of a product node has a different scope from the other children of the same product node.
This definition does not include any restriction enforcing normalization of the weights on a sum node, but the rest of the SPNs shown in this dissertation will. That is, given a sum node S and its weights w 1 , w 2 , . . . , w n , we have n i=1 w i = 1. Even if the weights of an SPN are not normalized, they can be normalized in linear time [PTPD15]. An SPN with all weights normalized is called a normalized SPN.

Value of an SPN
We represent an SPN with the symbol S, and use it to compute the value of an SPN given a partial assignment X as follows: • The value of an indicator associated with a variable X and value x is 1 if x ∈ X [X] and 0 otherwise.
• The value of a sum node is the weighted sum of the value of its child nodes (using the weight parameters on the arc between the sum node and the child node).
• The value of a product node is the product of the value of its child nodes.
• The value of the SPN S is the value of the node at the root.
We represent the operation of obtaining the value of an SPN S given a partial assignment X as S(X ). This value inference takes linear time in the size of the SPN [PD11]. Figure 2.1 contains an example of this inference on an SPN over three variables X, Y, Z with the partial assignment X = {0, 1, 2}, showing the original SPN at the top and the value of each node during inference at the bottom.
Given that obtaining the value of an SPN is a function of partial assignments to the variables of the SPN, and, on a normalized SPN, it is always a real number between 0 and 1, SPNs represents a joint probability distributions.

A note on marginalization
The term "marginalization" is used many times in this dissertation with a particular meaning that may confuse the reader given knowledge of its original definition. This section aims at explaining this term's differences of meaning when applied to SPNs and the MAP algorithms.
Suppose we have an SPN S over three discrete random variables X, Y , and Z. Using S, we can obtain the value of queries over the probability distribution of these three variables, such as P (X = x 1 , Y = y 2 , Z = z 2 ), which will be represented by S({x 1 } × {y 2 } × {z 2 }) or, more succintly, with S(x 1 , y 2 , z 2 ). However, there are also queries over the probability distribution without assigning specific values to some variables. For example, P (X = x 2 , Z = z 3 ) ignores the assignment of the variable Y , which means a probability over the event {X = x 2 , Z = z 3 } and considering all possible instantiations of Y. Given the definitions of the probability distribution, this is the same as y∈val(Y ) P (X = x 2 , Y = y, Z = z 3 ). This is often referred to as a query with a marginalized variable, and their value is obtained by the sum of other queries.
The previous sum, in SPNs, can be obtained in a single evaluation by setting every indicator for each of the marginalized values to 1. The notation for its value, y∈val(Y ) S(x 2 , y, z 3 ), therefore, can quickly become cumbersome when dealing with more than just one variable marginalized. For example: z∈val(Z) S(x 2 , y, z). Because of this we will follow the convention of using S(x 2 , z 3 ) to mean that, in this query to the SPN S, the value of the variable Y is marginalized.
Additionally, depending on the context, there will be situations in which the "marginalization" will not include all values from a given variable. Consider the partial assignments X such that there is a value y 2 ∈ val(Y ) and {y 2 } / ∈ X . If we marginalize in Y given the set X and the query {X = x 1 , Z = z 3 }, what this actually means is the sum y∈{y 1 ,y 3 } S(x 1 , y, z 3 ).
Therefore, whenever we use the term "marginalization" on a variable X, the context will make it clear whether we are marginalizing on val(X), or on X [X].

Learning SPNs from data
Simple SPNs, such as the examples in the figures, can be built by hand, but when dealing with datasets of different domains and applications, this will not be the case. In order to obtain the structure and weights of an SPN, we need a learning algorithm.
[GP13] introduced a simple, recursive algorithm for learning the structure of an SPN and the weights of its sum nodes. It is called LearnSPN and works as follows: The input is a set T of instantiations of variables V 1 , V 2 , . . . , V n ∈ V . The output is an SPN representing a distribution over V learned from T . The most important steps are: 1. Handle the case of a univariate distribution, where V = {X i }. In this scenario, the data instances will be used to generate a probability distribution estimation for the variable X i . • If they can be partitioned, run LearnSPN on each of the partitions V 1 , V 2 , . . . , V n where V = n i=0 V i and create a product node as a parent of the n SPNs returned. • If they cannot be partitioned into approximately independent subsets, then partition the set T into subsets of similar instances T 1 , T 2 , . . . , T n where T = n i=0 T i and create a sum node as a parent of the n SPNs returned, having weight equal to |T i | |T | for each T i . There are two major decisions related to the partitioning of variables and data instances to be made during the implementation of this algorithm: Approximate independence of variables, and instance similarity. The implementation used for the experiments in this dissertation used the Two-Way Likelihood Ratio 2 for the former, and the K-Medoids 3 algorithm for the latter.
For the Univariate distribution estimation (when |V |= 1), and since this dissertation will only run these experiments against datasets with discrete variables, we will generate a multinomial distribution with event probabilities obtained from the frequency of such events found in the data instances.

Maximum-a-Posteriori problem in SPNs
Completing the missing pixels of an image, finding the most probable classification, or simply finding the values that are more probable as assignments to a subset of the variables are all applications of the Maximum-a-Posteriori (MAP) Inference with SPNs, and simply describing them is enough to understand their usefulness, and thus, the usefulness of this type of inference.
Given a set of variables X = {X 1 , X 2 , . . . , X n }, a Query set of variables Q ⊆ X, an Evidence set of variables E ⊆ X, and a Marginalized set of variables M ⊆ X, with Q, E, M disjunct, and Q ∪ E ∪ M = X, a complete evidence for the variables in E represented by e, the MAP problem in an SPN S is defined as: The interpretation of this formula is the following: • There is a set of variables Q and we want to obtain assignments for all of them such that the value of the SPN is maximal. For example, they can be the missing pixels of an image, or the classification variable.
• There is a set of variables E with assignments e, and we cannot change their assigned values. They are interpreted as, for example, the given pixels of an image that are not missing, or the entire image for a classification task.
• There is a set of variables M , and we cannot change their values, but they are marginalized (either with all the values the variables can assume, or some subset of these values). They can be interpreted, in the previous examples, as the pixels or variables that are not important or can safely be ignored.
With these interpretations, we have that the MAP inference obtains the maximum probability configuration given a certain query, evidence, and marginalized sets. The search space for this problem is the product of the size of the value sets for every variable in The obtained configuration is, therefore, the complete assignment for the variables in Q.
The MAP inference is an NP-Hard problem [CMdC17] in the SPNs described so far 4 , meaning there are no polynomial-time algorithms that can obtain the exact solution to this problem in valid SPNs if P = N P .
However, we can still use algorithms that may reach a solution "close enough" to the optimal solution that their found results can still be used on the same applications 5 . This "close enough" measure is named performance guarantee and will be denoted with the symbol α.
Consider an SPN S and its MAP configuration c, where S(c) = v is the optimal MAP value. Following the definition in [WS11], an approximate algorithm that reaches an α-approximation of the optimal value v will reach a value v such that αv ≤ v ≤ α. According to [CMdC17], a provably achievable performance guarantee is 2 |S| for SPNs with height larger than or equal to 3, where |S| is the number of nodes in the SPN, and 0 ≤ < 1. 4 There are different kinds of SPNs, with specific structural restrictions that are learned with different algorithms. In particular, with Selective SPNs, whose learning algorithm is introduced in [PGD14], the MAP problem is not NP-hard since it can be obtained with Max-Product, an algorithm that runs in linear time on the size of the SPN and will be described on chapter 3 5 It is important to differentiate applications in this step since, for some of them, quickly obtaining a value that is close to the optimal is enough, and in those cases, how much time can be expended in return of a closer approximation becomes the critical trait we will measure.

Algorithms for MAP Inference in SPNs
There are two main classes of algorithms for the MAP problem in SPNs.
• Exact Algorithms will find the exact solution for the MAP problem.
• Approximate Algorithms will find a complete assignment that may not be the solution to the problem.
The approximate algorithms are generally faster than the exact algorithms and their solution may have use, depending on the problem. For example, image completion with an approximate algorithm can generate an output that looks good enough for the application, and does not need to spend more time in order to find a perfect solution that differs by a small amount of pixels.
Some algorithms are characterized as Anytime, meaning that they can be stopped during execution and still return a meaningful solution (in our case, a complete assignment). Giving more time to the algorithm would provide a better solution, even if it is still restricted to an approximate solution.
The given examples in this chapter are run on a small SPN (see figure 3.1) that illustrates the difficulties of reaching the MAP configuration with the simpler algorithms. 1 It is easy to see that the MAP configuration for this SPN is {X = x 1 , Y = y 1 } with the value 0.34.

Max-Product or Best Tree
Max-Product was introduced in the same article that presented Sum-Product Networks [PD11]. As noted in [Peh15], this algorithm is not guaranteed to find the MAP configuration on SPNs that are not either selective [PGD14] or augmented [PGPD16].
This algorithm runs in linear time on the size of the SPN (the same complexity as the evaluation of a partial assignment), making it especially useful to efficiently find lower bounds to the MAP value.
Max-Product(S) replaces all sum nodes in S with max nodes in order to perform a bottom-up search that greedily selects the highest output on each of them, returning a configuration. The major steps of the algorithm are: 1. Change all sum nodes in S to Max Nodes 2. Perform a bottom-up evaluation on the SPN, propagating the maximum-valued configurations upwards

Return the configuration at the root
We now examine these steps in detail.

Change all sum nodes in S to Max Nodes
A Max Node outputs the highest value of its child nodes, multiplied by the weight of that child node.
This step is purely symbolical, since the functionality of a sum node can be emulated on the algorithm implementation -that is, without creating a new node type.

Perform a bottom-up evaluation on the SPN, propagating the maximum-valued configurations upwards
A recursive step that makes use of the max nodes and evaluates the SPN in a way of searching for the maximum values of each node. Leaf nodes (indicators) have only one value, and thus it is returned. Product nodes return the product of their children. Max nodes return the maximum value of their children as calculated with the rules described on the previous step.
When this calculation reaches the root, each node will have an associated value it can output. This value is the Max-Product of each node.

Return the configuration at the root
The Max-Product value at the root is the Max-Product value of the SPN. Since the SPN is valid, this value is the result of the independent selection of indicators, one for each variable, and thus results in a valid complete assignment, 2 which will be the returned configuration.
The path used to retrieve the variable assignments is one particular induced SPN from the SPN S [ZPG16], also called a Parse Tree [MJT18]. Since this is the parse tree with the largest value, the algorithm is sometimes called Best Tree.

Max-Product Pseudocode
Algorithm 2: Max-Product Data: An SPN S Result: A complete assignment x 1 if S is a sum node with weights w and child nodes c then 2 return The partial assignment X = argmax c i ∈c w i S(Max-Product(c i )); 3 else if S is a product node with child nodes c then 4 return The union of partial assignments of Max-Product(c i ) for c i ∈ c; 5 else 6 S is an indicator for variable X i with value x j ; 7 return The partial assignment to variable X i , {x j }; 8 end

ArgMax-Product
Argmax-Product, first presented in [CMdC17] computes a complete assignment for the variables of an SPN, and the value obtained by evaluating the SPN on this assignment, guaranteeing an equal or better result than the one found by Max-Product. The difference lies in the checking step that slows the algorithm by calculating the value of a sub-SPN. An example execution of this algorithm is shown in figure 3.4, where the two values that need to be recalculated on sub-SPNs are shown as v 1 vs v 2 , where the blue one is chosen instead of the red one.
ArgMax-Product(S) Performs a bottom-up search on an SPN S, selecting for the maximum evaluated values of sum nodes, and returning a complete assignment, and the value obtained when this assignment is evaluated. The major steps of the algorithm are: 1. Perform a bottom-up search of maximum assignments on each node, measured by the value of each node 2. Return the complete assignment and value found at the root 3.1 We now examine these steps in detail.

Perform a bottom-up search of Maximum configurations on each node, measured by the value of each node
If the node is a leaf node (indicator) over variable X and value x, we return its maximum value: 1 and the partial assignment is represents: {x}.
If the node is a product node, we return the union of the maximizations found on its children, and the product of the values they found.
If the node is a sum node we test each partial assignment obtained from running Argmax-Product on its children, evaluating them as if the current child node was the root of an SPN. The assignment that, when multiplied by the weight associated with its sub-SPN root, returns the highest value, is stored (along with its value). As an example, a sum node with 3 children will have to test 3 partial assignments on each of the 3 child nodes (9 evaluations of sub-SPNs) in order to check which one reaches the largest value. However, each child node that arrived at an Argmax-Product solution has already been evaluated in order to have that configuration, and thus we have at least one configuration less to check for each child node.
This step is costly for sum nodes due to the amount of potential assignments to test. Consider the sum node S and its child nodes ch(S). The number of sub-SPNs to evaluate is, at most, n = |ch(S)|, each one taking linear time on the size of the sub-SPN rooted at the child node. If child nodes reach equal partial assignments, then this number can be lowered. 3 2. Return the complete assignment and value found at the root As in the Max-Product algorithm, this one will find a parse tree used to extract a complete assignment. However, unlike the previous algorithm, the value of the assignment has already been calculated and will be returned as well.

ArgMax-Product Pseudocode
Algorithm 3: ArgMax-Product Data: An SPN S Result: A complete assignment x and the value S(x) 1 if S is a sum node with child nodes c then 2 return The configuration X and S(X ) such that S(X ) = max c i S(ArgMax-Product(c i )); 3 else if S is a product node with child nodes c then

Local Search
Local Search(S, x) searches for local optimal complete assignments, starting with an initial complete assignment x, by modifying the associated value of one variable at a time, evaluating it in S.
It is also a simple algorithm to run after running one of the other approximate algorithms to improve on the assignment they found. When the task of MAP inference has some pre-allocated time and another algorithm finished before this limit, it could be useful to use the remaining time on an algorithm such as this one to attempt an improvement. The major steps of the algorithm are: 1. Evaluate the complete assignment x in S and store its value.
2. Choose a variable X i whose values val(X i ) have not been evaluated yet. If no such variables exist, jump to step 4.

For each
Store the complete assignment that reaches the highest value.
4. Return the complete assignment of highest value.
We now examine these steps in detail.
1. Evaluate the complete assignment x in S and store its value.
The complete assignment x can be obtained in many ways, such as a random selection, or another approximate algorithm. The value S(x) is the starting point of the search and, as expected from an algorithm whose goal is to find the local maximum, will heavily influence the last solution found.
2. Choose a variable X i whose values val(X i ) have not been evaluated yet. If no such variables exist, jump to step 4.
There are many heuristic choices to use when selecting the variable X i . However, due to their similarity to the heuristics used in the faster implementations of the Max-Search algorithm, and how much more detailed they are in their implementation when first introduced [MJT18], we will postpone further discussions on these heuristics to their own section.

For each
Store the complete assignment that reaches the highest value.
Each variable X i chosen will have the possible values on the set of val(X i ). The current best complete assignment (stored) will have the assignment to the variable X i changed to one of the values in val(X i ), their values recomputed, and compared with the previous one.
Each time this step is executed, S will be evaluated (at most) |val(X i )| times, which results, if no optimizations are considered, in the total running time of T (S) n i=1 |val(X i )|, where T (S) is the time to evaluate the SPN S and n is the number of variables.
4. Return the complete assignment of highest value.
The complete assignment with largest value stored will be returned, and if there is an use for it, the value obtained by evaluating this configuration has already been obtained on a previous step of the algorithm and can be returned as well.

Local Search Pseudocode Algorithm 4: Local Search
Data: An SPN S, a set of Query Variables Q, an initial complete assignment of the variables in sc(S), x Result: A complete assignment e 1 best ← S(x); 2 e ← x;

Beam Search
Beam search, as an approximate algorithm for MAP, was first presented in [MJT18] and described as an extension of the Hill Climbing algorithm for MAP inference in Bayesian Networks from [Par02].
[ Dar03] introduced arithmetic circuits that encoded the network polynomial of a Bayesian Network, and proposed differential methods of obtaining posterior marginal probabilities, a result that was used in [MJT18] to develop an algorithm to obtain the derivatives of an SPN. Since SPNs were built on the same ideas [PD11], the extension of this algorithm to valid SPNs was straightforward.
This algorithm allowed the development of the Beam Search method, and also enabled the development of some heuristics for the Max-Search algorithm, presented in the same article.
Given a partial assignment X over the variables in X, and an SPN S, we can choose one variable X ∈ X and a value x ∈ val(X) such that {x} ∈ X and calculate the following derivative: That is, the partial derivative of variable X on a value x ∈ val(X) over an SPN S with a current set of partial evidences X will calculate the value of the SPN on X with the assignment of X changed to {x}. The advantage of this method is that this choice of {x} can be made for all variables X i ∈ X, and all assignments to X i , at the same time, with the time complexity of an evaluation on the structure of the SPN [MJT18].
The beam search algorithm needs a couple extra parameters: b, called the beam size that will control the size of the solution space to investigate, and X , an initial partial assignment. The major steps of the algorithm are: 1. Obtain the derivatives for each value x ij ∈ val(X i ) for each X i in sc(S) given the current set of partial assignments X .
2. Select the b best assignments of variables to create b new sets of partial assignments.
3. Repeat steps 1 and 2 until there are no more marginalized variables to choose from.
4. Return the complete assignment with the highest value.
We now examine these steps in detail.
1. Obtain the derivatives for each value x ij ∈ val(X i ) for each X i in sc(S) given the current set of partial assignments X .
This step runs in linear time on the size of the SPN, following the referenced algorithm for finding the derivatives described in [Dar03]. We start with a set of partial evidences X , having either single assignments for some variables, such as {x i }, or having variables marginalized, such as {x 0 , x 1 , . . . , x n }.
Each value x ij for each variable X i will generate a new partial assignment having X i = x ij with the other variables having the same value (or values) as in X . The derivatives will result in the values for each of these generated partial assignments.
For example, given an SPN over the variables X and Y where val(X) = {x 0 , x 1 } and val(Y ) = {y 0 , y 1 , y 2 }, and starting with the partial assignments representing all possible values X = {val(X) × val(Y )}, the algorithm for finding the derivatives will find 6 different values, one for each assignment of each variable, which represent the probabilities of that assignment with the other variables marginalized, that is: S(x 0 ), S(x 1 ), S(y 0 ), S(y 1 ), and S(y 2 ). 4 2. Select the b best assignments of variables to create b new sets of partial assignments.
Continuing from the last example, the algorithm has obtained 6 values for derivatives. In this step, we select the b variable assignments with the largest values in order to continue the algorithm.
3. Repeat steps 1 and 2 until there are no more marginalized variables to choose from.
From this point the process will have b sets of partial assignments, each of them with at least one variable with a single value (instead of being marginalized). On the next steps, when referencing the use of the set of partial evidences, X , the algorithm should then use each of the b sets 5 when finding the derivatives.
This kind of search will then find more partial assignments to be tested, but after finding the values for the derivatives of the b partial evidence sets, it will select only b of them, stopping a potential exponential grow of the number of partial assignments to evaluate.
With each iteration, each set of partial evidences will have fixed one value for one additional variable, and therefore the maximum number of iterations this algorithm can have is |sc(S)|, the number of variables in this SPN.

Return the complete assignment with the highest value
Once the best b partial assignments have all variables set to a single value, we return the one whose evaluation results in the largest value.

K-Best Tree
The K-Best Tree algorithm was defined in [MJT18] and created after the empirical observation that the MAP configurations on some SPNs correspond to a complete assignment that is represented by one of its Parse Trees. However, not necessarily the Parse Tree with the largest value (these cases are already handled by the Max-Product algorithm).
In order to understand the association between complete assignments and Parse Trees, consider the example from the Max-Product algorithm shown in Figure 3.3. It shows only the nodes selected at the end of the algorithm run, and its indicators are limited to only one value per variable, representing a complete assignment for the problem. Removing the arc from the sum node to the indicator of X = x 0 and adding the arc from this same sum node to the Indicator of X = x 1 will generate a different Parse Tree, with a different complete assignment ({x 1 , y 1 }), and a different evaluated value at the root.
There are multiple Parse Trees that represent the same variable configurations 6 , and the K-Best Tree algorithm ranks them by their value, selecting only the K largest ones and returning the one with the largest value when evaluated on the entire SPN.
The major steps of the algorithm are: 1. Change all sum nodes in S to max nodes.
2. Do a bottom-up evaluation of the SPN propagating the K best assignments upwards.
3. Select the best of the K propagated assignments at the root.
We now examine these steps in detail.
1. Change all sum nodes in S to max nodes. This is the same step done in Max-Product, and with the same desired effect.
2. Do a bottom-up evaluation of the SPN propagating the K best assignments upwards.
The evaluation will now use max nodes instead of sum nodes to compute the results. More than that, given the type of the node, the results may vary: • Indicator nodes can only return value 1 indicating the variable on their scope assuming the value they represent. This is not changed on this algorithm.
• Max nodes will now return a list of the K maximum values obtained by multiplying their children values with their respect weights.
• Product nodes have disjoint scopes on each child, therefore they will receive a list of K separate assignments for each child. In order to obtain the value for this node, we need to use one assignment from each child node in order to create a complete assignment for the product node. For example, take a product node with 2 children, one returning {x 0 , x 1 } (supposing K = 2) and the other {y 0 , y 1 }. There are 4 combinations of the two possible values on each variable, and so 4 values will be computed on this node, but returning only the K = 2 best ones.
For more examples, and to consider the costs dependent on both the structure of the SPN and K: -Considering still K = 2 and 3 child nodes, it is easy to see that there will be 8 assignments to test. For 4 child nodes, there will be 16 assignments to test, or 2 4 . -Considering K = 3, and 2 child nodes, and the child nodes returning the assignments {x 0 , x 1 , x 2 } and {y 0 , y 1 , y 2 }, there will be 9 assignments to test. With 3 child nodes, there will be 27 assignments to test.
It is now easy to see that the number of assignments to test is at most K |ch(P )| , 7 although they are limited by the upper bound of the total amount of possible assignments on the scope of the node, that is X∈sc(P ) |val(X)|. Still, implementation-wise, we can limit the number of assignments to check even further by considering our options with the number K, since we are only interested in the K best assignments returned by this node. The values will always be the product of what was obtained from the child nodes, and it is easy to see that, for example, we will not need the assignment obtained by merging all the lowest-valued assignments from the children (It will never be among the K best values). We can achieve a faster running-time for this algorithm by only examining the largestvalued assignments from each child. Consider a product node with 2 children, A and B, the first one returning its K best assignments as A 1 , A 2 , . . . , A k and the second one as B 1 , B 2 , . . . , B k , in order of the largest-valued to the lowest. We will always return the assignment resulting from the merging A 1 × B 1 . After that, we examine other combinations in order of largest to lowest such as A 1 × B 2 and A 2 × B 1 , and so on. For each merging of configurations examined, A i × B j , we examine in the next steps A i+1 × B j and A i × B j+1 , on the ascending order of values. 8 3. Select the best of the K propagated assignments at the root.
At the root, the selected K complete assignments can have their pre-computed values (the ones used to sort the list of values) stored, but in order to obtain the correct value of each 3.2 configuration on the SPN, each of them will have to be evaluated on the proper structure, that is, using sum nodes instead of max nodes to calculate the value of each of the K complete assignments. This can be done in time proportional to K|S|.

K-Best Tree Pseudocode
Algorithm 6: K-Best Tree Data: An SPN S, an integer K Result: A complete evidence x 1 foreach Node N in S in reverse topological order do 2 if N is a sum node with weights w and child nodes c then

Branch-and-Bound algorithms: Max Search
Given the lack of algorithms for an exact solution of this problem, a solution that is capable of surpassing the trivial solution is a good improvement. One approach of doing so is a Branch-and-Bound algorithm, which has yet to be bested after the experiments in [MJT18].

The Branch-and-Bound paradigm
Branch-and-Bound is an algorithm design paradigm used to develop algorithms for solving integer and combinatorial problems. Currently, the only non-trivial, exact algorithm for solving the MAP problem in SPNs uses this paradigm, therefore it will be useful to know its general characteristics before delving into one of its implementations.
Branch-and-Bound has 3 main steps: 1. Setup the problem with a graph enumerating all the possible assignments to the variables, and an initial solution.
2. Choice of a node to expand and obtain a solution to the problem having additional restrictions based on the selected variable.
3. Prune the selected node if the solution is worse. Change the solution if a better one has been found.
4. Repeat steps 2 and 3 until no nodes can be expanded, and return the best found solution.
We now explain each step in detail.
• Setup the problem with a graph enumerating all the possible assignments to the variables, and an initial solution.
A simple example for this graph on two binary variables in figure 3.5. It is easy to see that, with more variables or variables with more than 2 values, this graph can grow exponentially. 9 Each node represents a set of assignments made to a subset of the variables of the problem. On the figure, the root has no assignments (it is an empty set), the nodes on height 1 have an assignment on variable X, and the nodes on height 2 have the assignments of their parent node and an assignment of the variable Y . Each leaf contains a complete assignment to all the variables in the problem, and each of them is a potential solution.
The initial solution can be any assignment easily obtainable (such as a random assignment to all variables), that we can represent with X , and the value obtained with this configuration, represented by x * . The initial solution will also be called the first bound to the problem, and since the MAP problem is a maximization one, this Bound is a lower bound -we know the solution to the problem will never be lower than this bound. We set the current best solution as X .
The algorithm depends on a division between nodes: active or inactive. Inactive nodes will not be used on subsequent steps. The initial division is to set the root node as active and the rest as inactive.
• Choice of a node to expand and obtain a solution to the problem having additional restrictions based on the selected variable.
By choosing a node to expand, we will have a restriction to the problem in the form of the set of assignments represented by that node. We can only choose active nodes to expand and, when the choice is made, the node becomes inactive for the rest of the execution.
The solution found for the problem with the restriction encoded in this node may be feasible to the original problem, such as a complete assignment obtained from one of the leaves. So, we have three choices for what to do on the next step.
• Prune the selected node if the solution is worse. Change the solution if a better one has been found.
-If the solution is feasible and is better (in our case, larger) than the bound, change the bound to this new value, and set the current best solution to this newly-found one.
-If the solution is worse than the current bound, we prune this node, ignoring any node that may come as a child to this one (and they will never be active for the rest of the algorithm).
-If the solution is not feasible, but the value found for this relaxation is larger than our bound, we activate the child nodes of the chosen node. For MAP problems, this value found for this node with a solution that is not feasible is called an upper bound.
• Repeat steps 2 and 3 until no nodes can be expanded, and return the best solution found.
The algorithm will then repeat until no active nodes can be found on the graph. At this point, the current best solution will be the best solution, and thus returned, and the value associated with that solution is the best value found for the problem.

Max Search and Branch-and-Bound applications
MaxSearch(S, X ) performs a recursive, depth-first, Branch-and-Bound search on the partial assignments X in order to find the complete assignment of maximum probability x * using the SPN S to find the value of each assignment. In other words, this algorithm finds x * such that x * = argmax x∈X S(x). The major steps of the algorithm are: 1. Obtain an initial assignment solution 2. Check for termination or select a variable to branch 3. Pruning of inferior search trees

Repeat steps 2 and 3 until termination is reached
We now examine these steps in detail.
1. Obtain an initial assignment solution.
The initial solution will be considered the best solution found until this moment of the execution. If this is not a recursive call -meaning it can only be the first call to the algorithmit is possible to set this solution to any x ∈ X .
On subsequent calls, the initial solution is passed on from the previous call. This solution is only changed when the algorithm reaches a leaf of the search tree.
It is useful to store this solution as an assignment x * and its value S(x * ) for future comparisons since this value needs time O(n) to be computed.

Check for termination or select a variable to branch
Termination occurs when there is only one assignment for each variable in X , that is, a complete assignment. In this case, the algorithm does not have to compare its initial solution value because, given the restriction described above, the best solution should be inside the set X . Therefore, it just returns the assignment x ∈ X .
If there is more than one assignment in X , then there is at least one variable X ∈ scope(S) being assigned different values in X . There are two main ways of making this choice: (a) Choose any variable that respects the precondition at random (a heuristic we call First Variable) (b) Choose the variable whose number of different values being assigned in X is the lowest (The heuristic described in [MJT18] as Ordering).

Pruning of inferior search trees
For each value x the variable X was assigned, we can obtain upper bounds to the remaining assignments in X that follow the restriction X = x. By comparing these upper bounds with the best value we found so far (from the best solution x * ) we can shrink the set X . The article [MJT18] gives two algorithms for doing this: (a) Marginal Checking: Compare the current best value S(x * ) with S(X ). If the current best value is higher, we can remove all assignments from X .
It is easy to see that this algorithm performs the evaluation of the SPN marginalizing the variables that have not already been added to the solution on the lower bound, and it fits more easily with the definition of the Branch-and-Bound step shown previously.
(b) Forward Checking: Compare the current best value S(x * ) with the derivatives of each individual variable assignment {x } in X . As described before with the Beam Search algorithm, all these values can be obtained at once with an algorithm that runs in linear time.
The derivatives in an SPN S give the following information: for each variable X and set of assignments {x 1 , x 2 , . . . , x n } in X , we have the other variables X \ X and their assignments marginalized, 10 and obtain the values of S((X \ X ) × {x i }) for each x i ∈ {x 1 , x 2 , . . . , x n }. We do not marginalize all the values that a variable can be assigned; only the variables for which there is a corresponding assignment in X . After each of these comparisons, we may obtain many individual assignments {x i } that we can remove from the set of assignments X . If we remove any of these, we must recheck the derivatives, since their values will have changed. If we did not obtain any assignment we can remove from X , that is, all values are equal to or greater than the current best value, we continue the algorithm with the remaining partial assignments. Unlike the Marginal Checking algorithm, Forward Checking is capable of pruning many nodes of the search tree in one pass, which speeds up the checking for upper bounds 11 .

Repeat steps 2 and 3 until termination is reached
If the modified partial assignments X obtained from step 3 is not empty, we can continue the search with the additional restriction X = x of the variable obtained in step 2. We, then, call the algorithm again, passing forward the best solution (and its best value), and X . This new set of assignments will have only one assignment for the variable X, and the next search will continue with a different variable, which ensures that, at some time, the algorithm will terminate.
At any time the algorithm can be stopped and return its current best solution, and therefore it is an Anytime Algorithm.

Max-Search Pseudocode
Algorithm 7: Max-Search Data: An SPN S, an initial complete assignment x, a set of assignments to search X Result: The MAP solution of S: x * 1 lower ← S(x); 2 if X is a complete assignment then 3 return x, the only assignment in X 4 end 5 X ← one variable from sc(S) with more than one assignment in X ; foreach Variable X ∈ X do end 10 end 11 end 12 until X is no longer changed ; 13 return X

Heuristics for variable selection
On line 7 of the Max Search algorithm, we select one assignment of {x i } for some variable X that still has more than one value to be chosen. Which value should we choose? This question would, in this algorithm, have the same answer as "Which variable and assignment would make the algorithm finish faster?". With Max Search, convergence of the result is related to finding the lower bound, and expecting the next found upper bounds to not surpass it. The original article for this algorithm [MJT18] described one heuristic for this choice, and left another one implicit: • Ordering (the described heuristic): Choose the variable that has the least amount of values to choose. So, in X there are the sets of values each variable X i can assume, represented by • Random Selection (the implicit one): Either choose the variable that can be selected the fastest, or a random variable (if random selection is faster). This would be the heuristic to compare as a baseline, since it is supposed to take the least time, in comparison to any other. 12 And in our work we have also developed and tested additional heuristics for comparison: • Lowest (or Largest) Marginal: With the derivative calculation on the forward checking algorithm, it is possible to obtain the marginal probabilities of every value still to be checked. For example, if we still have not set the value of a variable X, 13 and we are still searching on the set of assignments X , and X can have the value of either x 1 or x 2 , we are able to obtain both P ((X \ X) × {x 1 }) and P ((X \ X) × {x 2 }), which are the marginal values of these two assignments. 14 These heuristics will choose, from {x 1 , x 2 }, the assignment that has either the lowest or the largest marginal value.
• Lowest (or Largest) Entropy: We use the definition and, especially, the alternative view of entropy described in [KF09] in order to create a heuristic to select the variable and value. Considering a variable X with values {x 1 , x 2 , . . . , x n } still to be set, the entropy of X = x is H P (X = x)log 1 P (X=x) . 15 . This value can be interpreted as a measure of uncertainty about the assignment of {X = x}. These heuristics will select either the variable and value with the lowest or the largest entropy.
We have described the algorithms and the MAP problem we want to examine. Now we can begin describing the experiments we did in order to answer the questions at the introduction.

Experiment Setup
The SPNs created for these experiments, with one exception, have been learned with repeated executions of the Learn-SPN algorithm, each one testing a different combination of parameters, being K ∈ {2, 3, 4, 5}, the number of clusters, and p ∈ {0.0001, 0.001, 0.01, 0.1, 0.0002, 0.002, 0.02, 0.2, 0.0005, 0.005, 0.05, 0.5, 0.0007, 0.007, 0.07, 0.7, 0.0009, 0.009, 0.09, 0.9} the p-value to be used on the test of independence. We performed a greedy-search on the learned SPNs by their log-likelihood, on each dataset. For the Optdigits dataset, we divided the training dataset into 10 parts, one for each digit, and learned 10 SPNs in the same way the other SPNs were learned. Afterwards, the main SPN was obtained by adding a Sum Node as a root having each of the learned SPNs as a child, with weight 0.1.
The table 4.1 shows relevant statistics of the SPNs used in the experiments. In order to understand the problem size, we need to point out that some SPNs are binary, and thus all their variables can be assigned either 0 or 1, and others are not, meaning that at least one variable can have 3 or more values. The size of the MAP problem, on an algorithm that lists all possible configurations for the X variables in an SPN will list S = X∈X |val(X)| configurations. This number gets exponentially large in some SPNs, and thus we decided to show ln(S) instead, simplifying the comparisons between these measures.
It is important to note that, given the scenarios we have chosen, the problem size will change with the test. For example, in an SPN with 30 variables, we may choose 10 for the Query (Q), 10 for the Marginal, and 10 for the evidence. The problem size for this specific scenario will be S = X∈Q |val(X)|. We have chosen 30 scenarios for each SPN such that each scenario would guarantee that the solution found by the Max-Product algorithm would not be optimal, and thus allowing for the other algorithms, given the allocated time, to find a better solution. For each scenario, we choose 1/3 of the variables for each set: Query, Marginalized, and Evidence. Where applicable, we have chosen these scenarios from the list of scenarios already chosen in [MJT18], but since we have additional SPNs in our experiments, we have generated our own scenarios for them. 1 The running time limit varies with the experiment, being 10 minutes for the experiments comparing the results and running time of each algorithm, and 5 minutes for the others. When the time runs out, for the algorithms that can be stopped and still return a meaningful value (such as Local Search and Max Search), the best obtained result will be used, while for the others, the result is 0.
For the algorithms that depend on parameters, such as Beam Search and K-Best Tree, we use b = 3 for beam size and K = 3 for the number of evaluated parse trees, both for the algorithm's efficiency (and thus being able to finish before the time limit) and in order to reduce memory usage.  For the Max Search algorithm, we follow the chosen heuristics at [MJT18] in order to reach somewhat comparable results to theirs. That is, we use ordering as the heuristic to choose the next variable and values to expand, and we use forward checking for pruning.
The experiments ran on a machine with the following configurations: 2 • CPU: Intel(R) Core(TM) i7-8700K CPU @ 3.70GHz • Memory: 64 GiB, L1 Cache: 384KiB, L2 Cache: 1536KiB, L3 Cache: 12MiB We obtained the datasets from the publicly available repository of Density Estimation Benchmark Datasets (DEBD) [LD10] [VHD12] [BDC + 15] [LM11] 3 and from the UCI Machine Learning Repository [DG17]. 4 With this setup and the following benchmark, we look to answer the first question from the introduction: given a limited time-frame, which algorithms perform best in comparison with Max-Product?  found by other algorithms in comparison with Max-Product, that is, given the Max-Product value is v, the other columns show the value found by the other algorithms, v as v v . It is important to note that, in time computations, we do not sum the time needed to obtain the initial configuration on the algorithms that need one. Local Search and Max Search both can use an initial solution to start their search. In order to obtain the initial solution, we use either Max-Product or Argmax-Product for Local Search, depending on the column, and Max-Product for Max Search. This may change the conclusions to be taken from the data, since the Local Search algorithm could be faster than Max-Product (see the MSNBC line on the time table) when the total time taken to obtain that solution would be closer to the sum of both these algorithms' running times.

Comparing Running time and found value
The rank is given by the relative "win" of the algorithm on each SPN, that is, having the largest value in proportion to the value found by Max-Product, or having the lowest running time. The algorithm that found the best value gets rank 1, and the algorithms that found the same value also get rank 1. The algorithm that finds the next best values get higher ranks, therefore the lower the rank value, the better. On the time table, the fastest algorithms get lower ranks accordingly.
The running time results were as expected, with Max-Product being overall faster than any other algorithm, and Max Search being the slowest. Argmax-Product being slightly slower than Max-Product conforms with the necessary reevaluations of smaller sub-SPNs.
The differences in both Local Searches are not accounted by the algorithm, and the difference in their averages is of about 1.3 seconds. In comparison to the other algorithms, they are close enough to dismiss this difference.
We conclude that, given the allocated time of ten minutes, running the Argmax-Product algorithm should be the best, and if more time can be allocated, using the found solution to run the Local Search algorithm looks, on average, to be the best way of obtaining a good solution. If the SPN is small enough for the Max Search algorithm to run until completion, then there is no reason to not use it, but few of our instances met this requirement.

Testing Heuristics for Max Search
From the test results on the last section, we know that Max Search is the slowest algorithm on average. However, Max Search is an Anytime Algorithm, and therefore we can stop it at anytime and obtain a usable solution. The question then becomes the second questions listed at the introduction: Which heuristics in Max Search find good solutions faster?
The ways we have to test this is by changing the heuristics for choosing the next variable to be expanded, and for the nodes to be pruned. Table 4.4 shows the results of our experiments with the following heuristics: First Variable (FV), Lowest and Largest Marginals (LoM and LaM), Lowest and Largest Entropy (LoE and LaE), Ordering (Ord), Marginal Checking (MC), and Forward Checking (FC). Each "win" marked on the table signifies an algorithm that reached, in the time limit of 5 minutes, the best value in the least amount of time, that is, given their final value, we give the win to the combination of heuristics that reached such value the soonest. If two or more heuristics reached the same value at approximately the same time (with a margin of error of ten milliseconds).
The results vary depending on the SPN. For a small SPN such as NLTCS, the algorithm runs so fast that most heuristics will receive wins, but the heuristics that involve less calculations (such as getting the first variable, or simply calculating the marginals) reach the best solutions first.
One surprising result was the heuristics that we chose for the experiments in the previous section, and the heuristics used in [MJT18] perform, on average and in total, worse than most other combinations of heuristics. The variable-choice heuristics that require less calculations (First Variable, Lowest Marginal, Largest Marginal) tend to perform better than the others. For the pruning heuristic, Marginal Checking still receives more wins in both total and average.
Evidence set as given, fixed variables. This is a different way of evaluating the algorithms since both the obtained value and the complete assignment are important in order to compare the solution's quality.
We measured three classifiers: A Support-Vector Machine, a Random Forest, and a Decision Tree, all from the scikit-learn library [PVG + 11] with their default parameters. We have chosen the Digits dataset, a small dataset (229 test instances) of binary images of digits for classification, and the Optdigits dataset, a larger dataset (1797 test instances) with the same purpose but with values ranging from 0 to 16. 6 We also learn an SPN from the same training dataset, and use it in order to auto-complete selected test instances with values removed.
In order to formulate a MAP problem with these instances, we remove the first 1 3 values and add the variables corresponding to them to the query set. We then remove the next 1 3 values and add the variables to the marginalized set. The rest of the variables, including the classification label, are given as evidence. Table 4.5 shows the base accuracy of each tested classifier, that is, the accuracy of classification with no missing instances and no auto-completion. We can consider this an upper limit to the classification accuracy on the auto-completed instances.
Besides the Max-Product and Argmax-Product algorithms we described before, we also have compared: • Random Completion: we complete the missing values with a random value, chosen uniformly from the set of values from each variable.
• Mode Completion: we complete the missing values with the value that occurs the most on the training dataset for that variable.
On tables 4.6 and 4.7 we have the accuracy as the percentage of correctly-classified images for each classifier, using each algorithm for auto-completion.
We can see that using the SPN and one of the MAP algorithms is, in the majority of scenarios, better than using random values or the mode for completing instances with very few exceptions. A smaller but inconsistent difference is in the accuracy provided by Argmax-Product. It shows that, sometimes on the Digits dataset, we may get a lower accuracy rate. On the Optdigits dataset, it is consistently higher accuracy. We know from the experiments on the previous sections that the values obtained by Argmax-Product are, on average, higher than the values obtained by Max-Product, but it is hard to see if there is any relation between this higher value and the task of classification. Because of this problem, we have compared the log-values 7 of each solution generated by Max-Product and Argmax-Product against the mean square error of the auto-completed solution (which is a rough measure of how far from the original image the auto-completion is). The results are in figure 4.1. We use blue points on the chart to show instances where the difference in log-values is zero, and red points to show where the Argmax-Product reached a higher value.
On the digits dataset we cannot see a clear relation between the value being better and a lower error (in fact, sometimes the solution is much better but the error is also higher). However, we can see a relation between the log-value and the lower error rates for optdigits, which could partially explain the increase in accuracy we have seen before on table 4.7.

Final Considerations
We have developed a benchmark to test MAP algorithms for SPNs and executed it on a collection of them, obtaining running time and average value comparisons. Additionally, we have provided and tested additional heuristics for variable selection in order to improve the speed of Max Search algorithms in terms of how quick they reach their solutions.
At the introduction, we have listed three questions, and now we can expand on them: 1. Given a limited time-frame, which algorithms perform best in comparison with Max-Product?
Depending on the SPN, and on the specific MAP problem scenario being tested, some algorithms will invariably perform better. The application depending on this calculation will dictate which algorithm should be used. Max Search is optimal when the SPN is small or its heuristics are able to reach the final values faster, and in these situations it would be preferred.
However, either the Argmax-Product or a combination of a solution from Argmax-Product with Local Search performed better on average, for small or large SPNs, and thus it should be the ideal choice when Max Search cannot be used.

Which heuristics in Max Search find good solutions faster?
On average, we have seen that choosing either the first variable, or the one with lowest marginal, and using the marginal checking algorithm yielded the best results.
However, there may not have been enough test instances, and they have not reached the end of the algorithms execution, so the other heuristics could, at the end, obtain the final value (the solution to the MAP problem) faster than the ones we have listed.
3. Does a larger value for a MAP solution provide any gain in a downstream task (such as classification) that depends on that solution?
We do not have a conclusive answer for this, and thus more experiments are needed. On the digits dataset, the accuracy did not increase, but the opposite happened on the optdigits.

Suggestions for future work
There is still a space for discovering and implementing an algorithm capable of surpassing Max Search in terms of speed while still obtaining the exact solution to the problem. SPNs still lack a fast algorithm for solving this problem, but in order to develop something that surpasses Max Search, these are the things to be considered: • For an algorithm that works with bounds, such as the Branch-and-Bound, it needs a better way of obtaining these bounds. Forward Checking on Max Search can obtain as many upper bounds as there are variables and values still to be tested (and is able to prune them at the same time they are checked). Obtaining better bounds, lower than the marginalized values, would be enough for this. There is still space for experimenting with an algorithm that finds less bounds, but with better quality (higher lower bounds, or lower upper bounds).
• The algorithm was not developed to take advantage of parallelism, and this may be an advantage that could be explored. Splitting the partial assignments set X into two, and restricting the search on each of these splits to be conducted on different threads could speed up the algorithm, but it clashes with the effect of Forward Checking (that is capable of pruning from the entire X ) and therefore may not be ideal.
The effect of dataset differences must be given some attention. On [MJT18] they have shown that developing an algorithm more specialized for the binary SPNs they used works to reach faster executions of the Max Search algorithm. Further, our experiments show that, depending on the SPN, the MAP values will vary much more in proportion to the Max-Product value. SPNs such as CWebkb have achieved larger values on other algorithms to the order of 10 7 , while SPNs such as Ad, 10 −2 . It is fair to conclude that a study on specific datasets, and the way each algorithm behaves on it, may provide more conclusive information for the MAP problem on that specific application than a generalized study using unrelated SPNs.
The technique of changing (removing nodes from) the SPN during an algorithm execution presented on that same article should also be investigated, even if their results were not consistently good, and their conditions for applying not fully described.
Added to the differences of datasets, the differences in learning algorithms that generates these SPNs completes the full scenario of what can change the relevant attributes for these algorithms: Structure and parameters. More information about the effect of these two variables should be useful for the next studies.


Smooth nonlinear fitting scheme for analog multiplierless implementation of Hindmarsh - Rose neuron model

The Hindmarsh-Rose (HR) neuron model is built to describe the neuron electrical activities. Due to the polynomial nonlinearities, multipliers are required to implement the HR neuron model in analog. In order to avoid the multipliers, this brief presents a novel smooth nonlinear fitting scheme. We first construct two nonlinear fitting functions using the composite hyperbolic tangent functions and then implement an analog multiplierless circuit for the two-dimensional (2D) or three- dimensional (3D) HR neuron model. To exhibit the nonlinear fitting effects, numerical simulations and hardware experiments for the fitted HR neuron model are provided successively. The results show that the fitted HR neuron model with analog multiplierless circuit can display different operation patterns of resting, periodic spiking, and periodic/chaotic bursting, entirely behaving like the original HR neuron model. The analog multiplierless circuit has the advantage of low implementation cost and thereby it might be suitable for the hardware implementation of large-scale neural networks. ________________________________________

Circuit implementations of mathematical neuron models and large-scale neural networks on a hardware level have received much attention recently [21]. Based on the off-the-shelf analog components and on-chip digital devices, numerous analog and digital neuromorphic circuits were developed to simulate standalone neurons and coupled neural networks [22,23]. Due to the power-efficiency and compactness of analog circuits, mathematical neuron models were widely implemented in analog [24][25][26].
Based on numerous off-the-shelf discrete components, Bao et al optimized an analog circuit design to implement an adapting synapse-based neuron model [24], Behdad et al presented an experimental electronic neuron to achieve a complete Morris-Lecar model [25], and Arthur and Boahen developed a silicon-based integrate-and-fire neuron using a dynamical system approach [26].
However, when complex nonlinearities are involved, hardware circuit implementations for mathematical neuron models become extremely difficult. To solve this issue, a large number of digitally FPGA-based piecewise liner approximation approaches were proposed to implement neuromorphic circuits without multipliers [27][28][29][30][31][32][33]. Using the digital multiplierless implementation method, Jokar and Soleimani aimed at a calcium-based plasticity model [27], Gomar and Ahmadi focused on a biological adaptive-exponential neuron model [28], and Hayati et al sighted a two-coupled biological HR neuron model [29]. Meanwhile, Rahimian et al implemented a two-compartmental Pinsky-Rinzel pyramidal neuron model in digital [30], Imani et al investigated the multiplierless realization of a coupled Wilson neuron model in digital [32], and Haghiri et al employed a low-cost digital design to implement the noisy Izhikevich neuron model without multipliers [33].
These electronic neurons with analog and digital circuit implementations are capable to reproduce lots of different neuron dynamics and spiking/bursting behaviors that might appear in biological neurons [25,31]. Encouraged by the above piecewise liner approximation approaches, this paper presents a novel smooth nonlinear fitting scheme to implement an analog multiplierless circuit for the 2D/3D HR neuron model. Using this smooth nonlinear fitting scheme, the polynomial nonlinearities in the 2D/3D HR neuron model can be perfectly fitted by the composite hyperbolic tangent functions. Therefore, the presented smooth nonlinear fitting scheme is suitable for implementing electronic neurons in analog and the designed nonlinear fitting functions are smooth and continuously differentiable. By contrast, the piecewise liner approximation approaches reported in [27][28][29][30][31][32][33] are applied by the digitally FPGA-based electronic neurons, and the utilized piecewise linear functions are non-smooth and continuously non-differentiable.
The contributions of this paper are summarized as follows. (1) A novel smooth nonlinear fitting scheme is presented to multiplierlessly implement the 2D/3D HR neuron model and different neuron dynamics and spiking/bursting behaviors are reproduced. (2) An analog multiplierless circuit for the 2D/3D HR neuron model is designed and fabricated on a printed circuit board (PCB). Hardware experiments show that the analog multiplierless circuit can perfectly generate different operation patterns, such as the resting pattern, periodic/chaotic bursting patterns, and periodic spiking patterns with different oscillating frequencies.
The rest of this paper is organized as follows. Section 2 presents a smooth nonlinear fitting scheme and constructs two nonlinear fitting functions for the HR neuron model. Section 3 designs the circuit modules of the smooth nonlinear fitting functions and develops a hardware device for the analog multiplierless circuit. Finally, this paper is briefly concluded in Section 4.

Nonlinear fitting scheme for the HR neuron model
This section presents a smooth nonlinear fitting scheme to implement two polynomial nonlinearities in the HR neuron model.

Brief descriptions for the HR neuron model
Simplified from the Hodgkin-Huxley neuron model [4], the 2D HR neuron model was firstly introduced in 1982 [7], and it is described as 3 2 2 , , x y ax bx I y c dx y where x, y, and I represent the membrane potential, spiking variable, and steady current, respectively.
To exhibit the rich electrical activities of neuronal membrane potential, the 3D HR neuron model was then introduced [8] and it is expressed as , x y ax bx I z y c dx y z r s x x z where z stands for the bursting variable and x 1 represents the resting potential. The extra variable z is coupled into the 2D HR neuron model to regulate the current I. Generally, the model parameters in (1) and (2) are selected as a = 1, b = 3, c = 1, d = 5, x 1 = 1.6, s = 4, and r = 0.01 along with an adjustable current I [9]. Note that the parameter r is a small value and it is amplified to facilitate the selection of circuit parameters.
For the sake of simplicity, denote 3 2 2 ( ) 3 , After substituting the aforementioned parameters and the above equations into (1) and (2), the original 2D HR neuron model in (1) can be rewritten as ( ) , ( ) , and the original 3D HR neuron model in (2) can be rewritten as ( ) , ( ) , 0.04 0.064 0.01 , Then, the 2D/3D HR neuron model involves two polynomial nonlinearities: the cubic polynomial F(x) and the quadratic polynomial G(x).

Novel smooth nonlinear fitting scheme
To physically implement the 2D/3D HR neuron model in analog, the multiplier is a basic component for constructing the quadratic and cubic polynomials. However, the analog circuit with multiplier has higher implementation cost than that without multiplier. To reduce the implementation cost of the 2D/3D HR neuron model, we present a smooth nonlinear fitting scheme without any multiplier to implement the polynomial nonlinearities of the 2D/3D HR neuron model.
The hyperbolic tangent function tanh(·) in mathematics is smooth differentiable and has the upper and lower bounds. It can be easily implemented using some off-the-shelf analog electronic components. In the presented smooth nonlinear fitting scheme, two composite hyperbolic tangent functions are utilized to fit the two polynomial nonlinearities of the 2D/3D HR neuron model.
Corresponding to the cubic and quadratic polynomials F(x) and G(x) shown in (3), two smooth nonlinear fitting functions H 1 (x) and H 2 (x) are realized by the composite hyperbolic tangent functions and described as 1 1 where m j , κ j , and δ j with j = 1, 2, 3, 4, 5, as well as δ o1 and δ o2 are the parameters of the fitting functions H 1 (x) and H 2 (x), as listed in Table 1. Thus, H 1 (x) consists of three hyperbolic tangent functions while H 2 (x) contains two.   It should be stressed that the nonlinear fitting functions H 1 (x) and H 2 (x) using the presented nonlinear fitting scheme are smooth and continuously differentiable, whereas the piecewise linear approximation functions appeared in [27][28][29][30][31][32][33] are non-smooth and continuously non-differentiable.

Spiking and bursting in the fitted HR neuron model
The classical HR neuron model can show two operation patterns, namely the spiking and bursting. In the spiking pattern, the slow variable z is set as zero and the corresponding system is reduced to a 2D HR neuron model. In this way, the oscillating frequency is closely related to the steady current I. By contrast, in the bursting pattern, the slow variable z is what generates the bursting patterns. Correspondingly, the system is described by a 3D HR neuron model and its operation patterns have the mode transitions from resting to periodic bursting, to chaotic bursting, and finally to spiking patterns as the steady current I increases successively.
To demonstrate the feasibility of the presented smooth nonlinear fitting scheme, we introduce the two fitting functions H 1 (x) and H 2 (x) into the 2D/3D HR neuron model to replace the original nonlinear functions F(x) and G(x). Thus, the fitted 2D/3D HR neuron model is coined. To exhibit the nonlinear fitting effects, several representative values of the steady current I are selected for the original and fitted 2D/3D HR neuron models.
When the steady current is set as I = 0.5 and 2, respectively, the original and fitted 2D HR neuron models operate in the spiking patterns with different oscillating frequencies, and their membrane potentials are shown in Figs. 2a, b. As can be observed, the oscillating frequency increases as the steady current increases. Meanwhile, the two spiking patterns generated by the original and fitted 2D HR neuron models are basically consistent.  For different values of the steady current I, the 3D HR neuron model shows different operation patterns. The four representative values of the steady current are considered as I = 0.1, 2, 3.3, and 5, respectively. When imposing these steady currents on the original and fitted 3D HR neuron models, the membrane potentials can be induced and illustrated in Fig. 3. As can be seen, with the increasement of the steady current, the operation patterns undergo the resting, periodic bursting, chaotic bursting, and periodic spiking patterns sequentially. In the bursting patterns, the busters of the original and fitted 3D HR neuron models have slight difference but their dynamical behaviors are the same. As a result, the fitted 3D HR neuron model can be used to simulate the rich electrical activities of the original 3D HR neuron model, indicating the feasibility of the presented smooth nonlinear fitting scheme.
Since the original and fitted 2D HR neuron models only have the simple periodic spiking dynamics, we take the original and

Analog multiplierless circuit design and hardware experiment
This section designs an analog multiplierless circuit of the fitted 2D/3D HR neuron model using the circuit modules of composite hyperbolic tangent functions. Thereafter, a hardware device is developed for the analog multiplierless circuit.

Analog Multiplierless circuit design
The circuit modules for positive and negative hyperbolic tangent functions were designed in [24,35], as shown in Figs. 5a, b. They are designed using some off-the-shelf analog components, including operational amplifiers, voltage sources, transistors, and resistors. Note that the current source I 0 is implemented using two transistors, three resistors, and one voltage source [35], and its value is achieved by adjusting a resistance.
respectively, where E is a fixed offset voltage, μE = δ j , κ = R F /(2RV T ), and μ = R F /(2R E V T ) = κR/R E . In our design, the offset voltage is connected to the inverting input terminal of the operational amplifier via an adjustable resistor. The intent is to easily achieve different offset values of the hyperbolic tangent function by adjusting the resistance R E . However, when realizing the offset values of the hyperbolic tangent function, the offset voltage connection method reported in the literature [24,35] should adjust the offset voltage. This requires the power supply to provide multiple and high precision DC voltage outputs, and thus greatly increases the difficulty of hardware circuit experiment. In short, the circuit modules of positive and negative hyperbolic tangent functions in Figs. 5a, b are designed optimally in this paper.  With these circuit modules of hyperbolic tangent functions shown in Figs. 5a, b, two nonlinear fitting functions H 1 (x) and H 2 (x) given in (6) and (7) can be directly designed and their circuit modules are shown in Figs. 5c, d, respectively. When applying an input v i , the outputs of the two circuit modules in Figs. 5c, d can be represented as where m j = R/R mj with j = 1, 2, 3, 4, 5, V 1 = δ o1 = RE/R o1 , and V 2 = δ o2 = RE/R o2 .
The circuit parameters R = 10 kΩ, R C = 1 kΩ, E = 15 V, I 0 = 1.3 mA, and V T = 26 mV are determined. Using the parameter settings of the nonlinear fitting functions H 1 (x) and H 2 (x) shown in Table 1, the desired resistances of all circuit modules in Fig. 5 can be calculated and listed in Table 2. Thus, the cubic and quadratic polynomials F(x) and G(x) given in (3) can be approximately achieved using relatively cheap commercial electronic components.
Using the two circuit modules shown in Figs. 5c, d, the analog multiplierless implementation circuit of the fitted 2D/3D HR neuron model can be designed and drawn in Fig. 6. The upper part is two integral channels with the circuit modules of H 1 (v x ) and H 2 (v x ), and the lower part is an extra integral channel for implementing the slow variable v z . Based on three capacitor voltages v x , v y , and v z in Fig. 6, the circuit state equations under the physical time t = RCτ are established as where R I is a resistance to adjust the steady current.
2D HR neuron model 3D HR neuron model Compared (11) with (5), the resistance R I for adjusting the steady current I can be expressed as Thus, different spiking and bursting patterns can be displayed from the designed analog multiplierless circuit by adjusting the resistance R I achieved by a precise potentiometer.

PCB-based hardware experiments
Based on the circuit schematics shown in Figs. 5 and 6, a hardware device for the analog multiplierless circuit of the fitted 2D/3D HR neuron model is welded on a PCB level. Fig. 7 shows a picture of the hardware device. Here, the resistors, capacitors, precise potentiometers, MPS2222 transistors, and TL082CP operational amplifiers are employed. The off-the-peg ±15 V DC power module is used to provide the supply voltages of the operational amplifiers and offset voltages of the circuit modules. The experimental outputs are captured by the WaveSurfer 510 oscilloscope. Due to the limitation of saturation output level, it is necessary to ensure that the operational amplifier supplied by ±15 V DC power module operates within a linear operating range about ±13 V. In addition, I 0 can be adjusted by a specific potentiometer and one can refer to [35] for more details. Fig. 8 shows the screenshots of the hardware experimental prototype with the captured transfer characteristic curve for the circuit module H 1 (v i ) and periodic bursting pattern for the fitted 3D HR neuron model.  with the numerical simulations given in Fig. 1, manifesting the feasibility of the presented smooth nonlinear fitting scheme in the analog circuit implementation. Following the relation in (12), the resistance R I is employed for adjusting the values of the steady current I. First, disconnect the jumper and select the implementation circuit of the fitted 2D HR neuron model. Corresponding to Fig. 2, the spiking patterns for two representative values of R I are measured and shown in Fig. 10. Next, connect the jumper and select the analog implementation circuit of the fitted 3D HR neuron model. Corresponding to Fig. 3, the resting, periodic/chaotic bursting, and periodic spiking patterns for the four representative values of I are measured, as shown in Fig. 11. Consequently, various electrical activities can also be acquired from the analog multiplierless circuit of the fitted 2D/3D HR neuron model as well. The hardware experiments show that the presented smooth nonlinear fitting scheme can implement the 2D/3D HR neuron model without any multiplier.  In the hardware experiments, only the resistances in the circuit modules H 1 (v i ) and H 2 (v i ) need to be fine-tuned. Corresponding to the desired resistances listed in Table 2, the measured resistances are determined by Precision LCR Meter and listed in Table 3. By comparing the values in Tables 2 and 3, one can see that the measured and desired resistances have some slight deviations. These slight deviations are caused by the model idealizations, parasitic parameters, measurement errors, and so on.
The Mean Absolute Percentage Error (MAPE) introduced in [36] can be used to evaluate the deviations between the measured and desired resistances. Denote R desiredk as the desired value in Table 2 and R measuredk as the measured value in  (13) where N represents the number of resistance samples. According to the values listed in Tables 2 and 3, we can obtain that N = 17 and MAPE = 2.4% and this evaluation result is an acceptable precision for the PCB-based hardware experiments.

Conclusions
In this paper, we presented a smooth nonlinear fitting scheme for the 2D/3D HR neuron model. When implementing an analog circuit, the presented smooth nonlinear fitting scheme can avoid multipliers by employing two nonlinear fitting functions.
Therefore, without any multiplier, the analog circuit has low implementation cost, and thus is suitable for the hardware implementation of large-scale neural networks. Numerical simulations of the fitted 2D/3D HR neuron model and hardware experiments of the analog multiplierless circuit can reproduce various neuronal electrical activities with an acceptable error.
Certainly, the presented smooth nonlinear fitting scheme can also be applied in other neuron models or neural networks, such as the Morris-Lecar neuron model [25] and coupled Hindmarsh-Rose neuron model [29], which deserve our future study.