An Order Reduction Design Framework for Higher-Order Binary Markov Random Fields

The order reduction method is an important approach to optimize higher-order binary Markov random fields (HoMRFs), which are widely used in information theory, machine learning and image analysis. It transforms an HoMRF into an equivalent and easier reduced first-order binary Markov random field (RMRF) by elaborately setting the coefficients and auxiliary variables of RMRF. However, designing order reduction methods is difficult, and no previous study has investigated this design issue. In this paper, we propose an order reduction design framework to study this problem for the first time. Through study, we find that the design difficulty mainly lies in that the coefficients and variables of RMRF must be set simultaneously. Therefore, the proposed framework decomposes the design difficulty into two processes, and each process mainly considers the coefficients or auxiliary variables of RMRF. Some valuable properties are also proven. Based on our framework, a new family of 14 order reduction methods is provided. Experiments, such as synthetic data and image denoising, demonstrate the superiority of our method.


Introduction
The higher-order binary Markov random field (HoMRF) is a non-convex optimization model widely used in the fields of economy, information theory, quantum computing, machine learning and image analysis [1][2][3][4][5][6][7][8][9][10][11][12][13][14]. Recently, a new order reduction method has been developed to optimize HoMRF energies. The order reduction method transforms the HoMRF into a reduced quadratic binary Markov random field (RMRF) by elaborately setting the coefficients and auxiliary variables of RMRF to ensure the equivalence between HoMRF and RMRF. By equivalently transforming a complex HoMRF into an easier RMRF, the order reduction method achieves remarkable performance and thus has attracted increasing attention from researchers.
Since it is difficult to reduce an HoMRF into an RMRF in a single step, the mainstream approach is to transform a higher-order monomial into a sum of linear and quadratic monomials, and then iteratively perform this operation for all higher-order monomials of the HoMRF until a RMRF is obtained. According to the higher-order monomials to be reduced, previous order reduction methods can be divided into two categories: the methods of reducing higher-order negative monomials, and the methods of reducing higher-order positive monomials. In first category, Freedman et al. [15] developed an order reduction method for higher-order negative monomials for the first time. Then, Anthony et al. [16][17][18] proposed an order reduction method that is a more general form of the method [15]. Both two methods require only a binary auxiliary variable to produce submodular monomials, which can be easily minimized in polynomial time [19]. Anthony et al. [16][17][18] and Yip et al. [20] also developed two alternative methods to reduce higher-order negative monomials. Since these two methods produce nonsubmodular monomials that make minimizing RMRF to be an NP-hard problem, they are rarely used in practice.
Compared with reducing higher-order negative monomials, designing the order reduction method for higher-order positive monomials is much more difficult [1,2,[16][17][18]21,22]. Ishikawa [1,2] firstly developed a method of reducing higher-order positive monomials, by summarizing the rules in a large number of order reduction methods from the thirdorder to seventh-order positive monomials. Further, Anthony et al. [16][17][18] initiated a systematic study of lower and upper bounds on the number of binary auxiliary variables required by the order reduction methods. Nevertheless, Anthony et al. [16][17][18] only proved the existence of these upper and lower bounds, and failed to design any specific order reduction method that matches the corresponding upper and lower bounds. Until recently, Boros et al. [21,22] provided four order reduction methods, which utilize less auxiliary variables than Ishikawa [1,2]. Compared with the method of reducing higher-order negative monomial, the method of reducing higher-order positive monomials suffer from two critical problems: the nonsubmodular quadratic monomials obtained lead to the NP-hard problem of minimizing RMRF energies, and the introduced multiple binary auxiliary variables increase the computational complexity.
Therefore, designing more powerful order reduction methods for higher-order positive monomials is important and challenging. We observe that the main obstacle of designing order reduction methods for higher-order positive monomials is the interdependence between the coefficients and auxiliary variables of RMRF. Studies [21,22] also noticed the interdependence of auxiliary variables and coefficients. When designing order reduction methods that use logarithmic auxiliary variables, they proved that the coefficients increase exponentially. However, they did not further research the interdependence issue. Such interdependence forces us to consider the auxiliary variables and coefficients of RMRF simultaneously rather than separately, which greatly increases the design difficulty.
In this paper, we propose a unified design framework to reduce higher-order positive monomials. The proposed framework not only significantly decreases the difficulty of designing order reduction methods, but also provides 14 order reduction methods that allow applications in different fields to pick up their best method. Since the main difficulty of reducing HoMRF is due to the interdependence between the coefficients and auxiliary variables of RMRF, our core idea is to decompose this difficulty into two easier processes, where each process considers the coefficients or auxiliary variables separately. The flowchart of the proposed framework is shown in Figure 1. First, we generalize previous order reduction methods and propose a novel general reduction function (GRF) that requires only one auxiliary variable for any higher-order positive monomial. Different from the previous order reduction methods that utilize multiple binary auxiliary variables, our integral auxiliary variable is fixed, hence we can focus on setting the coefficients of RMRF. We rigorously prove the properties of the minimum value range of the integral auxiliary variable. Second, since the coefficients are already considered in the first process, we propose substitution and minimum transformations that convert the integral auxiliary variable of GRF into binary auxiliary variables. With the proposed transformations, we can concentrate on setting the binary auxiliary variables of RMRF. Based on the framework, a new family of 14 order reduction methods is proposed. Moreover, some state-of-the-art reduction methods can be easily derived from the proposed framework. We believe that our novel framework may take the design of order reduction methods to a new level. Comparison experiments show the superiority of our work.  The main contributions of our study are as follows: 1.
We propose an order reduction design framework that divides the complex design of order reduction methods into two simpler processes. Compared with existing works, our framework significantly decreases the design difficulty for order reduction methods.

2.
A novel GRF is developed to generalize previous order reduction methods. Unlike the previous methods with multiple binary auxiliary variables, GRF utilizes an integral auxiliary variable. Some valuable properties of GRF are also rigorously proved.

3.
Two sets of substitution and minimum transformations are developed to produce more order reduction methods. A variety of 14 order reduction methods are produced to enable applications in different fields to choose their most suitable method. Moreover, four state-of-the-art order reduction methods can be easily derived from our work.

Notation and Related Work
In this section, we provide some basic terminologies, as well as the notations and backgrounds necessary for understanding the subsequent sections.

Higher-Order Binary Markov Random Field
In this paper, we focus on the binary Markov random field, since a Markov random field with any multi-labels can be converted to a binary Markov random field via standard techniques [2,[23][24][25][26]. An HoMRF, also known as a pseudo-Boolean function [19,27], is a binary Markov random field whose order is greater than two. It is defined as E : {0, 1} n → R. Here, n is the number of variables and R is the set of real numbers. The energy function E can be uniquely expressed as multilinear polynomials [19]: where x 1 , . . . , x n are binary variables, [n] = {1, 2, . . . , n}, c S ∈ R is the coefficient of c S ∏ i∈S x i . In particular, when S = ∅, c S denotes the constant of E(x 1 , . . . , x n ). The size of S is denoted by |S|. If |S| > 2, the monomial c S ∏ i∈S x i is a higher-order monomial. For brevity, a general higher-order negative monomial is denoted as −x 1 x 2 · · · x |S| , and a general higher-order positive monomial is denoted as Although minimizing a general energy function is generally an NP-hard problem [19,28], the submodular energy function can be minimized globally in polynomial time. Proposition 24 in [19] is commonly used to judge whether E satisfies submodular.

Related Works of Order Reduction Methods
As shown in Definition 1, the goal of reducing the HoMRF E is to design an equivalent RMRF R for E. However, directly reducing E into R is extremely difficult. A mainstream order reduction method is to transform the higher-order monomials of E one by one until R is obtained. Note that an important property of the higher-order monomial c S x 1 x 2 · · · x |S| is symmetry, which means that it is invariant under any permutation of the coordinates {1, 2, . . . , m} of the variable x. Most of reduction methods [16,18,21,22] utilize the symmetry to reduce higher-order monomials based on the following concept: where x = (x 1 , x 2 , . . . , x |S| ) and y = (y 1 , y 2 , . . . , y m ) .
We introduce the order reduction methods of higher-order negative terms and higher-order positive terms, respectively. First, for the higher-order negative monomial, Freedman et al. [15] are the first to provide: This method produces a submodular [19] RF and requires only one auxiliary variable. Thus, it achieves excellent experimental performance [1][2][3][4]. A more general form of Equation (5) is proposed by paper [17]: where c 1 is the coefficient. To reduce the higher-order negative monomial, Anthony et al. [18] provide another method with one auxiliary variable, but their method do not satisfy submodular; alternatively, Yip et al. [20] propose a method that is an equivalent transformation of Equation (5). Second, to reduce the higher-order positive monomial, order reduction method generates a nonsubmodular RF and requires multiple auxiliary variables, which is the main reason for the difficulty in minimizing the RMRF. A lot of efforts have been devoted to designing the RF for a higher-order positive monomial. Ishikawa [1,2] observes the RF of x 1 x 2 x 3 , x 1 x 2 x 3 x 4 , . . . , x 1 x 2 x 3 x 4 x 5 x 6 x 7 and then summarizes the first order reduction method high order clique reduction (HOCR) for a general higher-order positive monomial: where m = |S| − 1 2 , c i,|S| = 1, |S| is odd and m = i, 2, otherwise.
The researches [21,22] systematically study the upper and lower bounds of the number of auxiliary variables for several classes of specially structured pseudo-Boolean functions, and provide two order reduction methods that require a logarithmic number of auxiliary variables. One is the logarithmic reduction type-1 (LogR-1) method: where m = log 2 (|S|) and c = 2 m − |S|. Another is the logarithmic reduction type-2 (LogR-2) method where m = log 2 (|S|) − 1 and c = 2 m+1 − |S|. However, they notice that Equations (8) and (9) require exponential coefficients, which have a negative impact on computational performance. Thus, the following linear reduction (LinR) method is developed: where |S| 4 m |S| 2 and c = |S| − 2m. Finally, they propose a square root reduction (SqrtR) method that nicely matches the lower bound in [17]: where m = 2 |S| + 1 is the number of auxiliary variables, c = m 2 , and λ is a large number. y i 0 and y i 0 must satisfy the condition c · i 0 · y i 0 + j 0 · y j 0 = |S|.

The First Process of the Proposed Framework: General Reduction Function (GRF)
In this section, we define a novel reduction function g(x, z) and demonstrate the convenience of designing it in Section 3.1. Then, we prove some critical properties of g(x, z) in Section 3.2.

General Reduction Function (GRF)
Based on Definition 2, setting the coefficients and auxiliary variables of RMRF is greatly simplified to designing f (x, y). However, designing f (x, y) is still very difficult, mostly owing to the interdependence between the coefficients and the auxiliary variables of f (x, y). Such interdependence has rarely been studied, yet it must be considered when designing f (x, y). Therefore, to avoid considering both coefficients and auxiliary variables, we propose a novel function g(x, z) with a fixed auxiliary variable z, which allows us to focus on designing the coefficient.
is the value space of z and |W| is the size of W.
Since the higher-order negative monomial has the excellent reduction method Equations (5) and (6), Definition 3 only considers the higher-order positive monomial. Furthermore, g(x, z) relaxes the restriction of f (x, y) in two ways: g(x, z) can be more than quadratic and z has not be binary. Since g(x, z) is symmetric about x 1 , . . . , x |S| , we denote an auxiliary function h : where g 1 : W → R and g 2 : is simple with only two variables. Thus, we can easily utilize it to formulate three specific forms of g(x, z).

Properties of GRF
The computational complexity of minimization algorithms, such as the sequential tree-reweighted algorithm (TRW-S) [44,45], depend heavily on the size of W. Therefore, it is natural to hope that the size of W is as small as possible. Before the discussion, we prove some important properties of g(x, z). Theorem 1. Let g(x, z) be the GRF of x 1 x 2 · · · x |S| and h(k, z) be the auxiliary function of g(x, z). Proof.

3.
Prove that g 2 (z) 0. We prove this case by contradiction. Suppose z = z 0 ∈ W such that g 2 (z 0 ) > 0 and min z∈W h(k, z) = h(k, z 0 ) = 1(k = |S|). We discuss it in two cases according to the value of k.
(a) If 1 k < |S|, then h(k, z 0 ) = 0. According to Equation (13), we have Directly from the definition of h(k, z), we can write Prove that g 1 (z) ≡ 0. We prove this case by contradiction. Suppose that g 1 (z) ≡ 0. Based on the definition of h(k, z) and Equation (13), there exists z 0 ∈ Z such that In other words, Theorem 2. Let g(x, z) be the GRF of x 1 x 2 · · · x |S| and h(k, z) be the auxiliary function of g(x, z). There exists z 0 ∈ W such that at most three integers 0 If k 3 = |S|, then c = 0 which contradicts Theorem 1; if k 3 = |S|, there exists a solution of Equation (17) Now we prove that z 0 ∈ W corresponds to at most three integers. We prove this by contradiction. Suppose that there exists z 0 ∈ W such that four integers 0 To hold the equivalence of Equation (18), k 4 must be less than |S|, i.e., k 4 < |S|. At this point, c must be zero, which contradicts Theorem 1. For the case that there are more than four integers, the proof is similar. Corollary 1. Let g(x, z) be the GRF of x 1 x 2 · · · x |S| and h(k, z) be the auxiliary function of g(x, z).
Proof. Based on Theorem 2, to make Equation (17) hold and not contradict the fact that c > 0, k 3 must be equal to |S| and the solution of Equation (17) is Theorem 3. Let g(x, z) be a GRF of x 1 x 2 · · · x |S| and W be the value space of z. The minimum size of W is |S| 2 .
Apparently, the size of W 3 in Equation (16) is |S| 2 that matches the minimum size of W. Although Equations (14) and (15) do not satisfy the minimum size of W, they can not only derive some state-of-the-art order reduction methods shown in the next section, but also help to better analyze our framework shown in the experiments.

The Second Process of the Proposed Framework: Transformation from GRF to RF
This section proposes two transformations that convert the GRF g(x, z) to the RF f (x, y) in Sections 4.1 and 4.2.

Minimum Transformation
Unlike Equation (19), the minimum transformation utilizes the minimum operation to obtain f (x, y): Compared with the proposed substitution transformation, the minimum transformation is capable of avoiding the exponential coefficients in Equation (22) and transforming Equation (14) into an RF. There are two specific forms. First, suppose a vector of auxiliary variables: y = (y 1 , . . . , y m ) ∈ {0, 1} m where m = |W| − 1. We define z(y) = y 1 + y 1 y 2 + · · · + y 1 y 2 · · · y m .

Experiments and Discussions
In this section, we compare the proposed method (GR-1, GR-2, GR-3, SR-1, SR-2, SR-3, SR-4, SR-5, SR-6, MR-1, MR-2, MR-3, MR-4 and MR-5) with state-of-the-art order reduction methods (HOCR, LogR-1, LogR-2, LinR and SqrtR) on synthetic data and image denoising. Table 1 summarizes the order reduction methods developed in this paper. To ensure the fairness of the comparison, we utilize the classical algorithm TRW-S [44,45] to minimize the RMRFs reduced by these order reduction methods, and conduct the comparison experiments with the same hardware and software environment (AMD Ryzen 5800U, RAM 16GB, Nvidia RTX 3050 and MATLAB 2022a). Table 1. Comparison of the proposed 14 order reduction methods for the higher-order positive monomial x 1 x 2 · · · x |S| .

GRF Transformation Number of Auxiliary Variables Type of Auxiliary Variables
GR-1 Equation (14

Synthetic Data Experiments
In the first experiment, we follow the approach in the paper [46] to generate a series of synthetic HoMRFs, which have the following form: where x 1 , x 2 , . . . , x n ∈ {0, 1}, c i ∈ R is the coefficient, and E S : {0, 1} |S| → R. We synthesize HoMRfs of order three to seven. For each order, we generate 500 HoMRFs with n = 50, 500 HoMRFs with n = 200, and 10 HoMRFs with n = 1000. The variables x 1 , x 2 , . . . , x n are randomly sampled from a uniform distribution, and the coefficient c i and the values of E S are randomly sampled from a standard Gaussian distribution. The energy results of each order achieved by different order reduction methods are shown in Table 2, Table 3 and Table 4, respectively. Note that for SqrtR, MR-4 and MR-5, λ = 5.  The order reduction methods below the solid line are designed by our framework. The order reduction methods below the solid line are designed by our framework.  The order reduction methods below the solid line are designed by our framework.
First, from Tables 2-4, it can be seen that for the third-order HoMRFs, the best energy performance is obtained by HOCR, LinR and MR-3; for the forth-order HoMRFs, the best energy performance is obtained by GR-3; for the fifth-order to seventh-order HoMRFs, the best energy performance is obtained by HOCR and MR-3. As shown in Equations (7), (16) and (33), these methods with the best performance all have fewer number of auxiliary variables and smaller coefficients, which is strongly supported by the research [18,21,22,46]. Although the performance of HOCR is identical to that of MR-3, designing HOCR is more difficult than MR-3. Since HOCR is designed by summarizing the RF f (x, y) of the third-order to seventh-order positive monomials [1,2], it is heuristic. Therefore, it must be fortunate enough to heuristically find the required RFs from a large number of RFs. As for LinR, it is designed by a complex mathematical theory in the references [21,22]. In contrast, our design framework explicitly gives two processes, each of which is much more rigorous than HOCR [1,2] and simpler than LinR [21,22].
Moreover, as we discussed in Section 4.2, the state-of-the-art order reduction methods LogR-1 and LogR-2 [21,22] can be derived from our design framework. We also present two methods SR-3 and SR-6 that are more concise forms of LogR-1 and LogR-2. For the energy performance, SR-3 and SR-6 outperform LogR-1 and LogR-2, respectively, shown in Tables 2-4. It shows that the proposed framework produces superior order reduction methods than state-of-the-art ones.
Third, to validate the importance of the proposed theorems in Section 3.2, we compare the energy performance of GR-1 Equation (14), GR-2 Equation (15) and GR-3 Equation (16). As seen in Tables 2-4, while GR-3 achieves the best energy performance on the fourthorder HoMRFs and the second best energy performance on the fifth-order to seventh-order HoMRFs, GR-1 and GR-2 perform poorly. This result shows that GR-3 is more preferable than GR-1 and GR-2. Furthermore, the RFs transformed from GR-3 are also more preferable than those from GR-1 and GR-2: SR-4, SR-5 and SR-6 are superior to SR-1, SR-2 and SR-3 in energy performance, respectively; MR-3 outperforms MR-1 and MR-2, and MR-5 outperforms SqrtR and MR-4. All these results indicate the importance of GR-3 and demonstrate that the size of value space |W| does significantly affect the efficiency of the reduction method.
Finally, we compare the energy performance of substitution transformations Equations (20)- (22) and minimum transformations Equations (30) and (34) in Section 3.2. As shown in Tables 2-4, SR-4, SR-5 and SR-6 have exactly the same energy performance on the third-order and fourth-order HoMRFs, since their formulas are the same in lower order HoMRFs. For the fifth-order, sixth-order and seventh-order HoMRFs, SR-4, SR-5 and SR-6 have different energy performance, and the performance of SR-5 is superior to that of SR-4 and SR-6. It demonstrates that for the GR-3, the substitution transformation Equation (21) is better than other transformations Equations (20) and (22). Similarly, it can be seen from SR-1, SR-2 and SR-3 that no substitution transformations Equations (20)- (22) has obvious advantages for GR-2. Furthermore, MR-1 and MR-3 show the better energy than SqrtR and MR-5 respectively in any case, and MR-2 shows the better energy than MR-4 in most cases. This result reveals that with the same GRF, the minimum transformation Equation (30) is more preferable than the minimum transformation Equation (34). We can conclude that no matter what transformation g(x, z) chooses, the fewer auxiliary variables and smaller coefficients lead to superior results.

Image Denoising
In this section, we conduct the second experiment on the benchmark [1,2], Fields of Expert (FoE) [47,48], which is a widely used HoMRF model for image denoising. FoE represents the prior probability of an image distribution as: where n is the number of pixels of the final image, a i is the value of the i-th pixel in the noise image, and A i * x i denotes the result of convolving the patch at pixel i with filter A i . The parameters A i and α i are learned from a database of natural images [47,48]. For more details, please refer to [1,2,47,48]. In the benchmark [1,2], 10 images from the widely known Berkeley Segmentation Database (BSDS500) [47][48][49] are selected. Gaussian noise with a standard deviation σ is added to each image. For efficiency, we shrink the image size by half. To make a fair comparison, we set the same starting point and maximum number of 100 iterations for order reduction methods. We utilize the energy and the peak signal to noise ratio (PSNR) and the calculation time to compare the performance of different reduction methods. The results are shown in Table 5. We also plot the curves of the energy and PSNR with iterations shown in Figure 2. A visual example of image denoising is illustrated in Figure 3. Note that for the FoE model, HOCR and LinR are the same. Thus, they have the same performance in terms of energy and PSNR and differ only slightly in terms of calculation time. The same is true for SR-4, SR-5 and SR-6. For SqrtR, MR-4 and MR-5, λ = 5. Table 5. Comparison of energy, PSNR and calculation time achieved by different order reduction methods, averaged over the benchmark [1,2]. We highlight the best performance in bold.

Methods
Energy   First, as shown in Table 5, the best performance of the energy, PSNR and calculation time is obtained by LinR, MR-5 and MR-3, respectively. For the denoised images with σ = 20, MR-3 has 1.68% higher energy than LinR and 0.89% lower PSNR than MR-5. However, the calculation time of MR-3 is 8.73% and 55.31% lower than that of LinR and MR-5, respectively. Thus, MR-3 balances energy, PSNR and calculation time and achieves the overall best performance. Similarly, for the denoised images with σ = {15, 25}, MR-3 also obtains the overall best performance. As discussed in the synthetic data experiment, MR-3 also has the advantage of lower design difficulty, which fully demonstrates the superiority of the proposed framework.
Moreover, as we discussed in Section 4.2, the proposed SR-3 and SR-6 are simpler and clearer than LogR-1 and LogR-2. It can be seen from the Figure 2a and Table 5 that SR-3 outperforms LogR-1 in energy, PSNR and calculation time, while SR-6 outperforms LogR-2 in calculation time. The results indicate that the proposed framework can improve the performance of some state-of-the-art order reduction methods.
Third, we compare different GR-1 Equation (14), GR-2 Equation (15) and GR-3 Equation (16) in Section 3.1. As shown in Figure 2b and Table 5, GR-3 has better performance than GR-1 and GR-2 in terms of energy, PSNR and calculation time. In Figure 3, GR-1 and GR-2 keep more noise, while GR3 preserves the important edge information while removing noise, which demonstrates the superiority of Equation (16). Similarly, MR-3 outperforms MR-1 and MR-2, and MR-5 outperforms SqrtR and MR-4, shown in Figures 2d and 3. This shows that the RF f (x, y) transformed from Equation (16) is superior to that from Equations (14) and (15). These results indicate that |W| significantly affect the efficiency of the reduction method and that proving the theorems in Section 2.2 is important.
Forth, we compare different transformations Equations (20)- (22), (30) and (34) in Section 3.2. As shown in Figure 2c and Table 5, SR-2 performs best among SR-1 and SR-3 in energy, PSNR and calculation time. The visual results in Figure 3 also show that SR-2 outperforms SR-1 and SR-3, which demonstrates that Equation (21) is more preferable than Equations (20) and (22). Similarly, MR-1 outperforms SqrtR in energy, PSNR and calculation time; MR-2 and MR-3 outperform MR-4 and MR-5, respectively, in energy and calculation time. This result demonstrates the significant advantage of Equation (30) over Equation (34) in energy and calculation time, while for PSNR, Equation (34) is slightly better than Equation (30).
Moreover, we analyze the stability and convergence of different order reduction methods. To numerically investigate the stability of order reduction methods, the input images in the benchmark [1,2] are added with Gaussian noise of different standard deviations σ = {15, 20, 25}. The results are shown in Table 5. As it can be seen, different order reduction methods show almost the same hierarchy of denoising performances for σ = {15, 20, 25}, which demonstrates that the order reduction methods exhibit a very robust response against the standard deviation changes. In Figure 4, to show the convergence of different order reduction methods, we plot the energy changes obtained by different order reduction methods for all images in the benchmark [1,2]. The experimental results show that all order reduction methods converge for all images in the benchmark [1,2]. Among these methods, the curves of LogR-1, SqrtR and SR-1 are not smooth enough, indicating that they are more susceptible to the influence of fusion move [1,2], which is a standard technique that converts Markov random field with any multi-labels into a binary Markov random field. Moreover, most of the proposed order reduction methods (GR-1, GR-2, GR-3, SR-2, SR-3, SR-4, SR-5, SR-6, MR-1, MR-2, MR-3, MR-4 and MR-5) converge smoothly to their respective energies, which illustrates the superiority of the proposed design framework.  Finally, we summarize two important conclusions that are helpful for designing effective reduction methods. First, the fewer the auxiliary variables y in f (x, y), the better the experimental performance. For example, as shown in Figure 2 and Table 5, LinR, LogR-2 and MR-3 with the least y in f (x, y) have excellent performance. Second, the smaller coefficients in f (x, y) are also beneficial to the experimental performance, especially the calculation time. For example, MR-3 has the fastest calculation time since it has the fewest coefficients y in f (x, y).

Conclusions
In this paper, we propose a novel framework that significantly reduces the design difficulty for order reduction methods. The framework generalizes previous order reduction methods and proposes a novel GRF, which is then transformed into more RFs. The experimental results validate the superiority of the proposed design framework. There are several interesting directions for future work. First, although the proposed order reduction methods are able to reduce HoMRF with any order, minimizing the RMRF reduced from HoMRF would fail when the order of HoMRF is particularly high. Therefore, the efficiency of order reduction methods is the priority for future research. Second, there is a few of applications for order reduction methods. We plan to introduce the order reduction method into more fields.