1 Introduction↩︎

Combinatorial optimization focuses on finding the best solution from a finite but often exponentially large set of possibilities. Within this domain, binary optimization, where decision variables take values in \(\{0,1\}\), plays a central role. Applications span scheduling, network design, facility location, portfolio selection, feature selection, and machine learning [1]–[4]. Since these problems are typically NP-hard, research has explored both exact approaches and heuristic algorithms to achieve a balance between computational tractability and solution quality.

Two standard formulations are Quadratic Unconstrained Binary Optimization (QUBO) and its generalization, Higher-Order Unconstrained Binary Optimization (HUBO). QUBO restricts the objective to linear and quadratic terms, benefiting from an extensive ecosystem of solvers, both classical and quantum [5]–[9]. HUBO, in contrast, allows higher-order interactions among three or more variables, offering richer modeling fidelity without introducing auxiliary variables. However, higher-order formulations may increase complexity and hinder scalability. Thus, the choice between HUBO and QUBO reflects a fundamental trade-off between expressiveness and efficiency.

A promising direction for such problems is optimization with tensor train sampling. Representing a probability tensor as a tensor train (TT) enables efficient probabilistic sampling in high-dimensional binary spaces without exhaustive enumeration. Because tensors naturally encode multi-way dependencies, they align well with HUBO’s higher-order interactions [10], [11]. The PROTES framework [12] demonstrates how tensor-based sampling achieves competitive performance in black-box optimization tasks [13]–[15]. This makes tensor sampling a useful setting in which to ask whether the algebraic simplicity of QUBO compensates for the dimension growth caused by reducing higher-order terms.

These observations lead to a central investigation:

What are the trade-offs between HUBO and QUBO formulations in probabilistic tensor sampling optimization, and how do these trade-offs evolve when the number of variables is fixed but the order of interactions increases?

1.1 Contributions↩︎

We formalize probabilistic binary optimization under HUBO and QUBO, highlighting structural trade-offs and reformulation costs. We conduct a systematic study of fixed-variable problems with varying interaction orders to isolate higher-order effects. Then we evaluate tensor sampling methods under both native and transformed formulations, comparing efficiency, accuracy, and robustness. Building on these results, we identify when HUBO’s expressiveness reduces the effective search dimension and when QUBO’s simpler form remains a more natural representation.

1.2 Paper Structure↩︎

The paper is organized as follows: Section 2 reviews the necessary background on HUBO, QUBO, and tensor-based optimization methods. Section 3 compares HUBO and QUBO formulations theoretically. Section 4 describes the experimental setup and presents the results. Section 5 discusses our results, and Section 6 concludes with future directions.

2 Background↩︎

2.1 HUBO and QUBO Formulations↩︎

A Higher-Order Unconstrained Binary Optimization (HUBO) objective of order \(k\) is \[E_{\text{HUBO}}(\mathbf{x}) = \sum_{p=1}^k \sum_{i_1<\dots<i_p} w_{i_1\dots i_p} \prod_{t=1}^p x_{i_t},\] where \(w_{i_1\dots i_p}\) are coefficients describing interactions among \(p\) binary variables. The order \(k\) of the HUBO corresponds to the polynomial degree of the objective.

The Quadratic Unconstrained Binary Optimization (QUBO) problem is the special case \(k=2\): \[E_{\text{QUBO}}(\mathbf{x}) = \sum_i a_i x_i + \sum_{i<j} b_{ij} x_i x_j,\] where \(a_{i}\) are linear coefficients and \(b_{ij}\) capture pairwise interactions. In practice, QUBO constraints are absorbed into the objective via penalty terms, and the objective can be expressed in matrix form: \[E_{\text{QUBO}}(\mathbf{x}) = \mathbf{x}^T Q \mathbf{x}.\]

Both HUBO and QUBO objectives define discrete polynomial optimization problems over \(\mathbf{x} \in \{0,1\}^d\). Directly enumerating all \(2^d\) configurations is infeasible for large \(d\), motivating specialized optimization strategies.

2.2 Tensor Representations of Probability Tensors↩︎

A binary optimization problem can be associated with a tensor representation. For HUBO/QUBO themselves, one could in principle construct an order-\(d\) tensor \[\mathcal{E} \in \mathbb{R}^{\underbrace{2 \times \dots \times 2}_{d}},\] where \(\mathcal{E}(i_1,\dots,i_d)\) stores the objective value at \(\mathbf{x}=(i_1,\dots,i_d)\). However, storing and manipulating this tensor is intractable for large \(d\) due to exponential growth.

Instead, tensor decompositions are applied to the probability tensor, which encodes a distribution over binary solutions.

2.3 Tensor Train Format↩︎

In the TT format, a high-dimensional tensor is represented as a chain of 3D cores: \[\mathcal{T}[n_1, n_2, \dots, n_d] = \sum_{r_0=1}^{R_0} \cdots \sum_{r_d=1}^{R_d} \mathcal{G}_1[r_0, n_1, r_1] \mathcal{G}_2[r_1, n_2, r_2] \dots \mathcal{G}_d[r_{d-1}, n_d, r_d],\] where \(\mathcal{G}_k \in \mathbb{R}^{R_{k-1} \times N_k \times R_k}\) are TT-cores, and \(R_k\) are TT-ranks. This representation reduces memory requirements from \(O(2^d)\) to \(O(d \cdot N \cdot R^2)\), where \(N = \max_k N_k\) and \(R = \max_k R_k\). Basic linear algebra operations (norms, contractions, etc.) can be performed in polynomial time in \(d\).

2.4 Tensor Train in Optimization↩︎

The key idea is that the probability distribution over solutions \(\mathcal{P}\) can be stored in TT-format. This makes it possible to efficiently sample solutions and update the distribution in high dimensions. Two leading tensor-based optimization methods are:

1. PROTES [12]: A probabilistic method that parameterizes \(\mathcal{P}\) in TT-format, samples candidate solutions, evaluates them, and updates \(\mathcal{P}\) to bias toward high-quality solutions. See Appendix 8 for complexity analysis.

2. TTOpt [16]: A deterministic method based on the maximum-volume principle. It identifies representative submatrices during tensor unfolding and uses them to approximate the maximum entry of \(\mathcal{P}\). See Appendix 8.

These methods are not applied to the HUBO/QUBO tensor directly, but rather to the probability tensor that guides the search over the exponentially large solution space.

A practical comparison of scalability is shown in Table 1.

3 Theoretical Comparison↩︎

In different combinatorial optimization problems, the choice between HUBO and QUBO can be more natural depending on the structure of the objective. In this section, we compare the two formulations from the point of view of expressiveness and reduction overhead. The main issue is not only whether a high-order polynomial can be quadratized, but also how many additional variables are introduced and how this changes the tensor dimension handled by the sampler.

3.1 Expressiveness↩︎

In the context of combinatorial optimization, expressiveness refers to the ability of a problem formulation to accurately represent interactions between variables. QUBO models are limited to pairwise interactions, i.e., terms of the form \(x_i x_j\), whereas HUBO formulations can include higher-order terms \(x_{i_1} x_{i_2} \dots x_{i_k}\) for \(k > 2\) and therefore represent dependencies among several objects directly. We provide examples to illustrate when a QUBO formulation arises more naturally, and when a HUBO formulation is more appropriate.

3.1.0.1 Max-Cut Problem.

Given an undirected graph \(G(V, E)\) with a set of vertices \(V\) and a set of edges \(E\), the Max-Cut problem asks for a partition of \(V\) into two sets such that the number of edges crossing the partition is maximized.

We can formulate this problem as a QUBO by introducing binary variables such that \(x_j=1\) if vertex \(j\) is in the first set and \(x_j=0\) otherwise. The expression \(x_i+x_j-2x_ix_j\) is equal to 1 only if exactly one of the two variables is equal to 1, and it is equal to 0 in all other cases [5]. Thus, this expression indicates whether the endpoints of an edge belong to different parts of the cut.

So, the Max-Cut problem can be formulated as:

\[\label{eq:max-cut32first} y = \max \sum_{(i,j) \in E} \left( x_i + x_j - 2x_i x_j \right)\tag{1}\]

We see that this is a quadratic form, so we can reformulate this problem as a QUBO:

\[y = \max_{\mathbf{x} \in \{0, 1\}^n} \mathbf{x}^TQ\mathbf{x}\]

The condition of a graph’s vertex set, such as the existence of an edge, is naturally formulated for two vertices; therefore, the QUBO formulation is the most interpretable in this case.

3.1.0.2 RSA Problem.

Let a large composite number \(N\) be given. The RSA problem consists of representing this number as the product of two prime numbers \[N = p \times q\]

We can easily reduce this to an optimization problem [18]:

\[\mathop{\mathrm{argmin}}_{p,q \in \mathbf{N}} \|pq - N\|_2 \label{eq:rsa}\tag{2}\]

\[\|pq - N\|^2_2 = p^2q^2 - 2pqN + N^2\]

Next, we will quantize the numbers \(p\) and \(q\) and substitute them into the objective functional:

\[p = \sum_{l=0}^{n-1} 2^l q_l, \quad q = \sum_{l=0}^{n-1} 2^l q_{n+l}.\]

\[\begin{align} p^2 q^2 &= \left( \sum_{l=0}^{n-1} 2^l q_l \right)^2 \left( \sum_{l=0}^{n-1} 2^l q_{n+l} \right)^2 \nonumber \\ &= \left( \sum_{l=0}^{n-1} 2^{2l} q_l + \sum_{l_1 < l_2} 2^{l_1 + l_2 + 1} q_{l_1} q_{l_2} \right) \left( \sum_{l=0}^{n-1} 2^{2l} q_{n+l} + \sum_{l_1 < l_2} 2^{l_1 + l_2 + 1} q_{n+l_1} q_{n+l_2} \right) \nonumber \\ &= \sum_{l_1=0}^{n-1} \sum_{l_2=0}^{n-1} 2^{2(l_1 + l_2)} q_{l_1} q_{n+l_2} + \sum_{l_1 < l_2} \sum_{l_3 < l_4} 2^{2l_1 + l_2 + l_3 + 1} (q_{l_1} q_{n+l_2} q_{n+l_3} + q_{l_2} q_{l_3} q_{n+l_4}) \nonumber \\ &\quad + \sum_{l_1 < l_2} \sum_{l_3 < l_4} 2^{l_1 + l_2 + l_3 + l_4 + 2} q_{l_1} q_{l_2} q_{n+l_3} q_{n+l_4}. \end{align}\]

\[-2pqN = -2N \left( \sum_{l=0}^{n-1} 2^l q_l \right) \left( \sum_{l=0}^{n-1} 2^l q_{n+l} \right) = \sum_{l_1=0}^{n-1} \sum_{l_2=0}^{n-1} (-2^{l_1 + l_2 + 1} N q_{l_1} q_{n+l_2}).\]

Minimising the function in this form provides an equivalent reformulation of the RSA problem as a HUBO problem.

In addition to their natural interpretability, HUBO and QUBO formulations often differ in the order of the tensor to which they can be reduced. Recall that HUBO and QUBO can be expressed as polynomial objectives in \(d\) binary variables, with degree two for QUBO. As discussed in Section 2, the probability tensor used by tensor sampling methods has one mode per binary variable. Therefore, a reduction that increases the number of variables directly increases the tensor dimension, even if the polynomial degree becomes lower. This makes it critical to determine how many variables \(x_i\) are required by each formulation and what practical tensor-size limits are imposed by the chosen sampling method (Table 1).

3.2 Reduction of HUBO to QUBO↩︎

A typical reduction introduces auxiliary variables \(y\) and penalty terms to enforce \(y = \prod x_{i_t}\). The overhead grows with the number and order of monomials, and this overhead is particularly important for tensor-based methods because each auxiliary variable adds a new tensor mode.

As was shown in [18], we can reduce the order-3 monomial \(cxyz\), where \(c \in \mathbb{R}\) in the following way:

\[\label{eq:3to2} c x y z = \begin{cases} c w(x + y + z - 2), & c < 0 \\ c \left\{ w(x + y + z - 1) + (xy + yz + zx) - (x + y + z) + 1 \right\}, & c > 0 \end{cases}\tag{3}\]

Continuing this idea, we can consider higher-degree monomials by temporarily grouping variables into a single one to obtain a cubic term. For example, we can reduce the monomial \(a_1 a_2 b_1 b_2\) with a positive coefficient as follows: \[\begin{align} a_1 a_2 b_1 b_2 &= a_1 \left( x_1 a_2 + x_1 b_1 + a_2 b_1 + a_2 b_2 + b_1 b_2 - a_2 - b_2 - x_1 + 1 \right) \\ &= x_2 x_1 + x_2 a_1 + x_1 a_2 + x_1 a_1 + x_1 a_2 - x_1 - a_1 - a_2 - x_2 + 1 \\ &\quad + x_3 x_1 + x_3 a_1 + x_1 b_1 + x_1 b_1 + a_1 b_1 - x_1 - a_1 - b_1 - x_3 + 1 \\ &\quad + x_4 x_1 + x_4 a_1 + x_1 b_2 + x_1 b_2 + a_1 b_2 - x_1 - a_1 - b_2 - x_4 + 1 \\ &\quad + x_5 x_1 + x_5 a_1 + a_1 a_2 + a_2 b_1 - a_1 - a_2 - b_1 - x_5 + 1 \\ &\quad + x_6 x_1 + x_6 a_1 + a_1 a_2 + a_2 b_2 - a_1 - a_2 - b_2 - x_6 + 1 \\ &\quad + x_7 x_1 + x_7 a_1 + a_1 b_1 + a_1 b_2 + b_1 b_2 - a_1 - b_1 - b_2 - x_7 + 1 \\ &\quad - a_1 a_2 - a_1 b_1 - a_1 b_2 - a_1 x_1 + a_1 \end{align}\] We observe that reducing a cubic monomial requires the introduction of one auxiliary variable, whereas a quartic monomial may require as many as seven. As the degree of the monomial increases, the number of auxiliary variables introduced by this reduction method grows quickly. Given the limitations of tensor sampling methods, this is critical for our purposes: a HUBO-to-QUBO reduction may produce a formally quadratic model that is nevertheless harder for tensor sampling because of the resulting increase in dimension.

3.3 Max-Cut formulation as a HUBO problem↩︎

Recall that, as a consequence of expression 3 , for \(x_{0}, x_{1}, x_{2}, x_{3} \in \{0,1\}\) and \(c \in \mathbb{R}\) with \(c<0\), the following equality holds: \[\label{eq:max-cut1} c x_{1} x_{2} x_{3} \;=\; c x_{0}x_{1} + c x_{0}x_{2} + c x_{0}x_{3} - 2c x_{0}\tag{4}\]

Thus, when \(c=-2\), the following equality holds: \[\label{eq:max-cut2} -2 x_{1} x_{2} x_{3} - 4x_{0} \;=\; -2 x_{0}x_{1} -2 x_{0}x_{2} -2 x_{0}x_{3}\tag{5}\]

This means that if we can identify terms of the form \(-2x_{0}x_{1} - 2x_{0}x_{2} - 2x_{0}x_{3}\) in our QUBO formulation 1 , then we can equivalently replace them with cubic and linear terms by increasing the order of the problem.

To this end, we modify the formulation of the Max-Cut problem. Let our graph \(G = (V, E)\) be defined on \(n\) vertices. Then \(E_{\text{full}}\) denotes the set of all possible edges \((i, j)\) with \(0 \leq i < j < n\) in this graph. Let \(M = E_{\text{full}} \setminus E\) be the set of edges that are missing from our graph.

First, we compute the objective function 1 for the case of a complete graph. Then, we successively take quadruples of coordinates of the vector \(x\), namely \(x_{i}, x_{i+1}, x_{i+2}, x_{i+3}\), and replace their occurrence in the objective function in the form \(-2x_{i}x_{i+1} - 2x_{i}x_{i+2} - 2x_{i}x_{i+3}\) with cubic and linear terms using expression 5 .

After that, we subtract from the obtained result the value of the objective function computed on the set \(M\). The resulting polynomial is equivalent to the original Max-Cut objective, but it contains cubic terms introduced by the replacement step. We illustrate this procedure in Algorithm [alg:Max-Cut95hubo].

Theorem 1: Algorithm [alg:Max-Cut95hubo] is correct.

Proof Indeed, since we only performed equivalent substitutions, even after increasing the order of the problem, we can write \(Q_{\text{full}}\) in the form \[Q_{\text{full}} = \sum_{(i,j) \in E_{\text{full}}} \bigl(x_{i} + x_{j} - 2x_{i}x_{j}\bigr) = \sum_{(i,j) \in M} \bigl(x_{i} + x_{j} - 2x_{i}x_{j}\bigr) + \sum_{(i,j) \in E} \bigl(x_{i} + x_{j} - 2x_{i}x_{j}\bigr).\]

Therefore, in the final step of the algorithm we have \[Q_{\text{sparse}} = \sum_{(i,j) \in M} \bigl(x_{i} + x_{j} - 2x_{i}x_{j}\bigr) + \sum_{(i,j) \in E} \bigl(x_{i} + x_{j} - 2x_{i}x_{j}\bigr) - \sum_{(i,j) \in M} \bigl(x_{i} + x_{j} - 2x_{i}x_{j}\bigr).\]

This simplifies to \[Q_{\text{sparse}} = \sum_{(i,j) \in E} \bigl(x_{i} + x_{j} - 2x_{i}x_{j}\bigr),\] which is exactly the original Max-Cut objective. \(\square\)

4 Experimental Setup and Results↩︎

We conducted experiments to evaluate tensor sampling performance on HUBO, QUBO, and transformed HUBO-to-QUBO or QUBO-to-HUBO formulations. Test cases included arbitrary polynomial functions, RSA factorization problems, and Max-Cut instances on random graphs. For each problem, we measured solution quality and computation time or evaluation budget across different problem sizes (\(n\)) and polynomial orders (\(k\)).

Both PROTES and TTOpt solvers were applied where appropriate to examine the trade-off between solution quality and computational budget, and to verify whether the observed results generalize across different optimization methods. For transformed QUBO formulations, auxiliary variables introduced by the reduction were included in the search space; therefore, the comparison reflects the effective dimension seen by the solver rather than only the polynomial degree of the objective. Convergence and best-found solutions were tracked over iterations, and the Max-Cut results are presented in Tables 6 and 7, as well as in Figures 1 and 2, for comparison.¹

4.1 Arbitrary Functions↩︎

We first compare HUBO and HUBO-to-QUBO formulations on synthetic polynomial objectives. These functions are small enough to make the effect of the reduction transparent, but they already contain cubic, quartic, and fifth-order interactions: \[f_1(\mathbf{x})= -3x_1x_3x_6 + 6x_2x_1x_3 - x_3x_4x_5x_6 + x_3x_6 + 2x_2x_5 - x_5 - x_1\] \[f_2(\mathbf{x})= 5x_1x_2x_3 - 4x_2x_3x_4 - 2x_1x_3x_4 + x_1x_3 + x_2 - x_4\] \[f_3(\mathbf{x})= -2x_1x_2x_3x_4x_5 - x_1x_3x_4x_5x_6 + 2x_1x_2x_3x_6 + x_2x_3x_4x_5 - 3x_1x_5 + x_2x_5\] \[\begin{gather} f_4(\mathbf{x}) = -2x_1x_2x_3x_4x_5 + x_2x_3x_5x_6x_7 - 3x_2x_3x_4x_5 + x_1x_2x_3x_4 \\ - x_3x_5x_6x_7 + 2x_1x_2x_3 - 2x_2x_3x_6 + x_1x_6 - x_6x_7 - x_7 \end{gather}\]

For the first four functions, both formulations return the same objective value, although the reduced QUBO is sometimes faster because the problem size is still small. However, this behavior does not persist for all high-order objectives. For example, for the third-degree function \[f_5(\mathbf{x}) = \sum_{i=0}^{n-1} \sum_{j=i+1}^{n-1} \sum_{\substack{k=j+1 \\ \text{k is odd}}}^{n-1} (-1)^j (j - i) x_i x_j x_k\] the output and time become different for some values of \(n\).

For fourth-degree functions, the situation is more complex because degree reduction with a positive coefficient adds seven auxiliary variables in the construction considered above, and a negative coefficient also increases the number of variables. Consider the function \[f_6(\mathbf{x})= \sum_{i=0}^{n-1} \sum_{j=i+1}^{n-1} \sum_{\substack{k=j+1 \\ k\equiv j+1 \pmod{2}}}^{n-1} \sum_{\substack{p=k+1 \\ p\equiv k+1 \pmod{2}}}^{n-1} (-1)^p (j - i + k) x_i x_j x_k x_p\] After the HUBO-to-QUBO transformation, PROTES does not recover the correct answer for the tested cases, while the native HUBO formulation does. This indicates that the additional variables introduced by the reduction can dominate the apparent advantage of a quadratic objective.

4.2 RSA problem↩︎

As shown in 2 , the RSA factorization problem can be reduced to a HUBO optimization problem: \[\mathop{\mathrm{argmin}}_{p,q \in \mathbf{N}} \|pq - N\|_2\]

For this problem, an objective value of zero corresponds to exact factorization. After reducing the HUBO formulation to QUBO, PROTES does not reach the global minimum in the tested cases, while the native HUBO formulation does. Results for different values of \(n\) are shown in Table 5.

As shown, even for these small examples PROTES fails to obtain the exact answer after the QUBO reduction. In addition, the HUBO formulation is faster for the larger tested instances, which is consistent with the dimension-growth argument.

4.3 Max-Cut problem↩︎

In Tables 6 and 7 we compare the results of QUBO and HUBO (obtained using Algorithm [alg:Max-Cut95hubo]) for the Max-Cut problem. The optimization was carried out using PROTES and TTOpt. In addition to the solution found by each method, we also report the budget required to obtain it. For all experiments, random graphs were used, and their characteristics are provided in the tables. Figures 1 and 2 illustrate the comparison of QUBO and HUBO Max-Cut for \(n=64\) and \(n=128\).

The Max-Cut results differ from the RSA and synthetic high-order cases. Since the standard Max-Cut objective is naturally quadratic, the original QUBO formulation does not require auxiliary variables and is often competitive. The HUBO version obtained by Algorithm [alg:Max-Cut95hubo] can still match QUBO quality in many instances, but it does not provide the same structural advantage as in problems whose native model is already high-order. This supports the main interpretation: the preferred formulation depends on whether a transformation reduces or increases the effective dimension faced by the tensor sampler.

5 Discussion↩︎

Our theoretical and empirical investigation reveals a clear and consequential trade-off between HUBO and QUBO formulations within the context of tensor sampling optimization. The central finding is that for problems with inherent high-order interactions, the HUBO formulation is not just more expressive but often significantly more efficient.

The Computational Advantage of HUBO: The prevailing assumption that QUBO is universally more tractable is challenged by our results. The HUBO-to-QUBO reduction process introduces a large number of auxiliary variables, drastically increasing the effective dimensionality of the optimization problem. Since the cost of tensor sampling methods scales sharply with dimensionality (see Appendix 8), this expansion often nullifies any theoretical advantage of the quadratic form. As our experiments show, for problems like RSA factorization and high-order synthetic functions, the native HUBO formulation consistently found better solutions (Tables 3, 4, and 5) and, in several cases, did so with a lower computational budget (shorter runtimes, fewer evaluations). This demonstrates that preserving the native high-order structure creates a simpler, lower-dimensional search space for tensor methods to explore, leading to superior performance.

When QUBO Remains Relevant: The QUBO formulation retains its value for problems that are naturally quadratic (like the standard Max-Cut problem) or where the number of variables is low. In these cases, as shown in Tables 6 and 7, both formulations can find solutions of similar quality. The choice may then hinge on the maturity and wider availability of QUBO-specific solvers outside the tensor sampling paradigm.

The Critical Role of Dimensionality: The effectiveness of any formulation is dictated by the capabilities of the underlying optimizer. Methods like PROTES and TTOpt are powerful but have strict practical limits regarding the number of variables they can handle (Table 1). Therefore, the decision between HUBO and QUBO is, in practice, a decision about which formulation produces a tensor of manageable dimensionality. Our work shows that for high-order problems, the native HUBO formulation is the only viable path to a solution, as the QUBO-reduced form often exceeds the feasible dimension for tensor samplers.

Limitations: The experiments are intended to isolate the effect of formulation and dimension growth rather than to exhaustively benchmark all available solvers. The conclusions are therefore most directly applicable to tensor-sampling optimizers and to reductions that introduce many auxiliary binary variables. Other QUBO reductions, penalty choices, or specialized QUBO solvers may change the practical balance for particular applications. Nevertheless, the dimension-growth mechanism identified here remains relevant whenever the optimizer scales strongly with the number of binary variables.

6 Conclusion↩︎

This work presented a comparative analysis of Higher-Order (HUBO) and Quadratic (QUBO) Unconstrained Binary Optimization models for probability optimization using tensor sampling techniques. We demonstrated that the traditional pathway of reducing complex problems to QUBO can be counterproductive when using modern tensor methods. The introduction of auxiliary variables inflates the problem dimension, often pushing it beyond the feasible scope of samplers like PROTES and TTOpt and resulting in poorer solutions. Our central conclusion is that for problems with inherent high-order interactions, the HUBO formulation is often superior: it provides a more accurate representation of the problem structure and, crucially, can lead to more computationally efficient optimization. This refines the conventional wisdom that QUBO is always the more tractable choice and provides a dimensionality-aware framework for selecting an optimization formulation. Future work should extend the comparison to additional quadratization schemes, larger benchmark families, and hybrid solvers that can exploit both native high-order structure and mature QUBO optimization tools.

The study was supported by the Ministry of Economic Development of the Russian Federation (agreement No. 139-10-2025-034 dd. 19.06.2025, IGK 000000C313925P4D0002).

7 Tensor Train for optimization↩︎

7.0.0.1 What is Tensor Sampling?

Probabilistic distributions are often high-dimensional, and storing them explicitly is computationally infeasible due to the curse of dimensionality. To address this problem, we can approximate the distribution as a high-dimensional tensor and represent it in the Tensor Train (TT) format. In this way, we avoid explicitly storing the full distribution while still being able to efficiently generate random samples (vectors) from it.

7.0.0.2 TTOpt.

The TTOpt [16] algorithm is designed to efficiently locate the maximum or minimum elements of a high-dimensional tensor \(\mathcal{T} \in \mathbb{R}^{N_1 \times \cdots \times N_d}\) that is defined implicitly through a multivariate function \(T(\boldsymbol{x})\). Unlike gradient-based methods, TTOpt requires no derivatives and is well-suited for black-box, discrete, or non-smooth optimization problems in high dimensions.

This method is based on the following max-volume observation [19]: if \(\hat{\boldsymbol{T}}\) is an \(R \times R\) submatrix of maximal volume (in selected rows and columns), then the maximal (by modulus) element \(\hat{T}_{\text{max}} \in \hat{\boldsymbol{T}}\) bounds the absolute maximal element \(T_{\text{max}}\) in the full matrix \(\boldsymbol{T}\): \[\hat{T}_{\text{max}} \cdot R^2 \geq T_{\text{max}}.\] Numerical experiments show that even this bound is pessimistic, that is why TTOpt and other maxvol-based methods aim to find the submatrix with a large volume instead of directly searching for the largest element.

We begin by considering the first unfolding of the tensor, \(\boldsymbol{T}_1 \in \mathbb{R}^{N_1 \times (N_2 \cdots N_d)}\), which reshapes \(\mathcal{T}\) into a matrix by grouping all but the first mode. Then, we sample \(R_1\) random column indices \(I_1^{(C)}\), each corresponding to a multi-index \((n_2, \dots, n_d)\), and construct a submatrix \(\boldsymbol{T}_1^{(C)} \in \mathbb{R}^{N_1 \times R_1}\) by sweeping over all values of \(n_1 = 1, \dots, N_1\).

Next, TTOpt applies the maxvol algorithm to \(\boldsymbol{T}_1^{(C)}\) to select a subset of \(R_1\) rows that form a submatrix \(\hat{\boldsymbol{T}}_1 \in \mathbb{R}^{R_1 \times R_1}\) of maximal volume. The corresponding row indices \(I_1^{(R)}\) are stored and used in the next step.

The algorithm reinterprets the selected \(R_1\) rows as a new implicit matrix \(\boldsymbol{J}_2 \in \mathbb{R}^{(R_1 N_2) \times (N_3 \cdots N_d)}\), which corresponds to reshaping the partially selected tensor slice. We then sample \(R_2\) random columns \(I_2^{(C)}\) in this new matrix and compute the submatrix \(\boldsymbol{J}_2^{(C)} \in \mathbb{R}^{(R_1 N_2) \times R_2}\) explicitly.

Again, the algorithm applies maxvol to extract \(R_2\) row indices \(I_2^{(R)}\), which correspond to combinations of the previously selected indices and new values along the second mode. This process is repeated sequentially across all modes \(k = 3, \dots, d\), alternating between random column sampling and maxvol-based row selection.

After reaching the last mode, the algorithm proceeds backward from mode \(d\) to mode \(1\), now treating previously fixed row indices as variables and selecting new column indices. This alternating forward and backward procedure (referred to as sweeps) continues until convergence or until a predefined budget of function evaluations is exhausted.

After \(T\) sweeps, the algorithm returns the multi-index \((n_1^*, \dots, n_d^*)\) corresponding to the largest observed value of \(T(\mathbf{n})\), along with the value itself.

Remark: The convergence of this method, in the case when it is applied to tensors, has not been proven. We can only be sure that the result will improve monotonically.

7.0.0.3 Optima-TT

This method [17] finds the maximum or minimum of a tensor given in the TT-format.

The main idea of this algorithm is to associate each tensor element with the values of the probability density function of a random vector \(\boldsymbol{\xi} = \{\xi_1, \ldots, \xi_d\}\), where \(\xi_i \in \{1,2,\dots,N_i\}\) and \(N_i\) is the size of the corresponding mode (\(N_i = 2\) in our case).

To avoid negative values, we need to square the tensor element-wise: \(p(n_1, \dots, n_d) = C (\mathcal{T}[n_1, \dots, n_d])^2\) for all \(n_k = 1, 2, \dots, N_k\) (\(k = 1, 2, \dots, d\)), where \(p\) is the probability distribution of \(\boldsymbol{\xi}\) and \(C\) is a normalization constant.

In this formulation, finding the maximum modulo element in the tensor \(\mathcal{T}\) is equivalent to finding the most probable value of \(\boldsymbol{\xi}\).

The algorithm builds the probability distribution, follows the [20], and sequentially samples the vector \(\boldsymbol{\xi}\) from the corresponding marginal distribution at each step.

\[\xi_i \sim \boldsymbol{p_i}(\xi_i|\xi_1, \dots, \xi_{i-1})\]

At each step, we keep the \(K\) most probable candidates. If in the \(l\)-th step of the algorithm we have \(\{\tilde{\xi}_1, \dots, \tilde{\xi}_{l-1}\}\), then the marginal distribution function \(p_l\) is given by the following expression:

\[p_t(\xi_l \mid \bar{\xi}_1, \dots, \bar{\xi}_{l-1}) = \sum_{n_{l+1}=1}^{N_{l+1}} \cdots \sum_{n_d=1}^{N_d} p(\bar{\xi}_1, \dots, \bar{\xi}_{l-1}, \xi_l, n_{l+1}, n_{l+2}, \dots, n_d).\]

In this case, efficient computation is ensured by the following theorem.

Theorem 2:

Let \(\mathcal{T}\) be a tensor with the TT-representation. Let the TT-cores of this representation be orthogonalized. Then the result of the convolution of this tensor with itself by the last \(l\) indices (\(0 < l < d\)) is given by the following explicit expression: \[\sum_{n_{d-l+1}=1}^{N_{d-l+1}} \sum_{n_{d-l+2}=1}^{N_{d-l+2}} \cdots \sum_{n_d=1}^{N_d} \mathcal{T}[n_1, \dots, n_d] \mathcal{T}[n_1, \dots, n_d] = \left\| \mathcal{G}_1[1, n_1, :] \mathcal{G}_2[:, n_2, :] \cdots \mathcal{G}_{d-l}[:, n_{d-l}, :] \right\|_2^2.\]

7.0.0.4 PROTES.

The PROTES [12] (Probabilistic Optimization with Tensor Sampling) method for black-box optimization is based on probabilistic sampling from a probability density function given in the low-parametric TT format.

This method, following [11], represents the distribution \(p(x)\) as a tensor \(\mathcal{P}_{\boldsymbol{\theta}} \in \mathbb{R}^{N_1 \times N_2 \times \dots \times N_n}\) (in our case \(N_1=\dots=N_n=2\)), where high-quality solutions should receive high probability. \(\mathcal{P}_{\boldsymbol{\theta}}\) is generated in a low-parametric representation with a set of parameters \(\boldsymbol{\theta}\), having the same shape as the target tensor \(\mathcal{T}\).

This method starts from the random tensor \(\mathcal{P}_0\) and, at each iteration of the algorithm, generates \(K\) trial points \(x_i\) from the current distribution \(\mathcal{P}_\theta\), represented in the TT format.

Among them, \(k\) points with the smallest values of the objective function \(f(x_i)\) are selected. These candidates have indices \(\mathcal{S} = \{s_1,\dots,s_k\}\) in the set \(\mathcal{X}_K\) of sampled candidates.

The parameters \(\theta\) are updated using gradient ascent, maximizing the logarithm of the likelihood of the selected candidates:

\[\widehat{L}_{\boldsymbol{\theta}}(\{\boldsymbol{x}_{s_1}, \boldsymbol{x}_{s_2}, \dots, \boldsymbol{x}_{s_k}\}) = \sum_{i=1}^k \log \left( \mathcal{P}_{\boldsymbol{\theta}}[\boldsymbol{x}_{s_i}] \right).\]

8 Complexity Analysis↩︎

8.0.0.1 TTopt

method has the following asymptotic complexity: \[\mathcal{O}\left(T \cdot d \cdot \max_{1 \leq k \leq d} (N_k R_k^2)\right),\] where \(T\) is the number of sweeps, \(d\) is the dimension, \(N_k\) is the size of the \(k\)-th mode (in our case \(N_k=2\) for every \(k\)), and \(R_k\) is the TT rank.

8.0.0.2 Optima-TT

method has the following asymptotic complexity: \[\mathcal{O}\left(d \cdot K \cdot N \cdot R^2\right),\] where \(d\) is the tensor dimension, \(K\) is the number of selected candidates for each TT-core, and \(N\) and \(R\) are the typical mode size and TT-rank, respectively.

8.0.0.3 PROTES

has the following asymptotic complexity: \[\mathcal{O}\left( M \cdot d \cdot \left( \frac{k}{K} \cdot k_{gd} \cdot \overline{R}^2 + (\overline{N} + \overline{R}) \cdot \overline{R} + \alpha(N) \right) \right),\] where \(M\) is the number of requests to the target function, \(d\) is the tensor dimension, \(k\) is the number of best candidates with minimal objective value, \(K\) is the number of sampled candidates, \(k_{gd}\) is the number of gradient ascent steps, \(\overline{R}\) and \(\overline{N}\) are the average TT-rank and mode size, respectively, and \(\alpha(N)\) is the complexity of sampling from the generalized Bernoulli distribution with \(N\) outcomes.

References↩︎

---

1. We removed Optima-TT because it can only solve small problems; see Table 1.↩︎