PROTES: Probabilistic Optimization with Tensor Sampling
Authors: Anastasia Batsheva
Andrei Chertkov
Gleb Ryzhakov
Ivan Oseledets
Venue: NeurIPS 2023 poster (2023)
Keywords: Tensor Train
Black Box Optimization
Sampling
Optimal Control
TL;DR: We propose a new black-box optimizer based on sampling from tensor train probability distribution and outperform other popular algorithms on 20 optimization problems.

Abstract: We developed a new method PROTES for black-box optimization, which is based on the probabilistic sampling from a probability density function given in the low-parametric tensor train format. We tested it on complex multidimensional arrays and discretized multivariable functions taken, among others, from real-world applications, including unconstrained binary optimization and optimal control problems, for which the possible number of elements is up to $2^{1000}$. In numerical experiments, both on analytic model functions and on complex problems, PROTES outperforms popular discrete optimization methods (Particle Swarm Optimization, Covariance Matrix Adaptation, Differential Evolution, and others).

PROTES: Probabilistic Optimization with Tensor Sampling
Anastasia Batsheva"
Andrei Chertkov"
Skolkovo Institute of Science and Technology
Skolkovo Institute of Science and Technology and AIRI, Moscow, Russia
Moscow, Russia
a. batsheva@skolleehh.rr
chertkovv@kklleechhho
Gleb Ryzhakor
Ivan Oseledets
Skolkovo Institute of Science and Technology
Skolkovo Institute of Science and Technology and AIRI, Moscow, Russia i. osseledets@skoltechh..
Moscow, Russia
g @.ryyhaakov@kkottech..t
Abstract
We developed new method PROTES for black-box optimization, which based the probabilistic sampling from probability density function given in the low-parametric tensor train format. We tested it on complex multidimensional arrays and discretized multivariable functions taken, among others, from real-worke applications, including unconstrained binary optimization and optimal control problems, for which the possible number of elements up 2 In numerical experiments, both on analytic model functions and on complex problems, PROTES outperforms popular discrete optimization methods (Particle Swarm Optimization, Covariance Matrix Adaptation, Differential Evolution, and others).
1 Introduction
The multidimensional optimization problem is one of the most common in machine learning.
includes the important case of gradient free discrete optimization [30 20
41)
a Pmmin = min f(x:) S.f. x
[M1, n,. na).
(1. 2, N.).
(1)
where d is the dimensionality of the problem, and N1, N₂, NJ are the numbers of items for each dimension. Such settings arise when searching for the minimum or maximum element in an implicitly given multidimensional array (tensor") including when considering the discretization of functions from continuous argument. Multidimensional discrete optimization problems are still computationally difficult in the case of complex target functions or large dimensions, and efficient direct gradient-fres optimization procedures are highly needed.
The development of methods for low-rank tensor approximations has made it possible to introduce fundamentally new approaches for the approximation, storage, and operation with multidimensional
'Equal contribution.
tensor is multidimensional arry with several dimensions (d 2 1). A two-dimensional tensor (d = 2) is a matrix, one-dimensional (d = 1) vector. We use lower case and upper case (where this docs not lead to confusion) letters for scalars, bold letters (a, b, ... for vectors. upper case letters (A, B, C....) for matrices, and calligraphic upper case letters (sl, B, B, C....) for tensors with
The (n1, 722 ... na))th entry of d-dimensional tensor y E IR NoNNxX,, NgN11NN2x... X......x Na is denoted by = y (1), 72, ...... nd), where 1, 2, N; (i 1,2, d), and N; is size of the i-th mode. The mode-i slice of such tensor is denoted by / [[m....,,n-1,.11 -
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

tensors [3
6 6
40 One of the common methods for low-rank approximation is the tensor train (TT)
decomposition 26 which allows bypassing the curse of dimensionality. Many useful algorithms (e. element- wise addition, multiplication, solution of linear systems, convolution, integration, etc.)
have effective implementations for tensors given in the TT-format. The complexity of these algorithms turns out to be polynomial in the dimension d and the mode size N. It makes the TT-decomposition extremely popular in a wide range of applications, including computational mathematics [29 [ and
machine learning
31 19
In the last few years, new discrete optimization algorithms based on the TT-format have been proposed: TTOpt [37] Optima- TT (5), and several modifications [34, 35 24 (6). However, the development of new, more accurate, and robust TT-based methods for mu multidmmnnioonll discrete optimization is possible. In this work, we extend recent approaches for working with probability
distributions and sampling in the TT-format [TO, 25 to the optimization area. That is, we specify a multidimensional discrete probability distribution in the TT-format, followed by efficient sampling from it and u its parameters to approximate the optimum in a better way. This makes possible to build an effective optimization method, and the contributions of our work are as follows:
We develop new method PROTES for optimization (finding the minimum or maximum value) of multidimensional data arrays and discretized multivariable functions based on a sampling method from a probability distribution in the TT-format;
We D for various analytic model functions and for several multidimensional QUBO and optimal control problems demonstrate its performance and compare it with popular discrete optimization algorithms (Particle Swarm Optimization, Covariance Ma- trix Adaptation, Differential Evolution, and NoisyBandit) as well as TT-based methods (TTOpt and Optima-TT). We used the same set of hyperparameters of our algorithm for all expeeieeenee experiments and obtained the best result for 19 the 20 problems considered.
2 Optimization with probabilistic sampling
Our problem is to minimize the given multivariable discrete black-box function f [ It can be formulated in terms of the multi-index search in an implicitly given d-dimensional tensor
NER.NN......... XNN, y [n1, n22 ...,, nal f(z)) x x= [m], 12. ..... nal.
for all n. 1,2, NI ...., d). The essence of our idea is to use a probabilistic method to find the minimum Emin- We propose establishing a discrete distribution p(x) from which the minimum could be sampled with high probability. This distribution can be specified tensor Po E IR Ni X X N2 in some low-perametric representation with set of parameters A, having the same shape as the target tensor 3/.
We start from a random non non-negative tensor Py and iteratively perform the following steps until the budget is exhausted or until convergence (see graphic illustration in Figure D):
1. Sample candidates of 2min from the current distribution Py: Ig = {x1 Z2..... o&
2. Compute the corresponding function values: 1/1 f(xx), 1/2 f(x2). .... WK f(xx);
3. Select & best candidates with indices & = {51 .,..... } from Ix, with the minimal objective value, i.e., Vj S yu for all j 8 and J E {1, 2, .... K} \ 8;
4. Update the probability distribution Py - ((new)) increase the likelihood of selected candidates S. We make several (kiga) gradient ascent steps with the learning rate À for the tensor Pg. using the following loss function
Lo (f Lzz. Do = Σ log (Folxx,1))
(2)
1=1
"As will be shown in the work, we use several gradient ascent iterations for the likelihood to update the parameters of the IT-decomposition. However, it is important to note that the proposed optimization method gradient-free, that is, it does not use the internal structure and gradient of the objective function. Further, for concreteness, we will consider the minimization problem in this paper. while the proposed method can be applied the discrete maximization problem without any significant modifications.
"The program code with numerical examples is available in the repository https: https://gtt //gith/ com/ anabatsh/PROTES

---

Sample
Compute
Probability
K multi-indices
black box values
tensor
y/1
X2
Black Box
N
XK
VK
XK
yK
Update
Select k top items
for selected multi-indices
Figure 1: Schematic representation of the proposed optimization method PROTES.
After a sufficient number of iterations, we expect the tensor Pg to represent an almost Kronecker delta-funce function with pronounced peak in the value of the minimum of the target function (or several peaks if the minimum is not unique). Therefore, this value will be sampled during the steps of our algorithm since the probability of sampling other values will be sufficiently small.
From low-parameter representation Pg expect efficient sampling algorithm and efficient calculation procedure for the logarithms in (2) with automatic differentiation capability to enable gradient ascent methods. As will be shown below, the TT-representation of tensors satisfies these requirements. Further, we will omit the index 0, assuming that the parameterized representation of the tensor Se corresponds to the TT-format. Note that the values K, kgd, À and the number of parameters in (i. e., rank of the ITT-decomposition) are the hyperparameters of our algorithm.
3 Basic properties of the tensor train format
Let us dwell on the concept of the TT-format. A d-dimensional tensor 9P E NN.. is said be in the TT-format 26
if its elements are represented by the following formula
R
R
Ra-1
F0 n1, n2, nal
[1,,n,, Be(ri, n2, r2]
Ba nd, 1), (3)
T(-1=1
where (ni, T2,,.... na) is a multi-index (n: = 1,2,.... N₂ for = 1, .)) integers Ro, .I...., Ra (with convention Ro Ra 1) named TT-ranks, and three-dimensional tensors '3 E RR Rt-1 xN.xit, (i (i=10 1, d) are named TT-cores. The TT-decomposition (3) (see also an illustration in Figure 2 allows to represent tensor or discretized multivariable function in a compact and descriptive cow-parameter form, which is linear in dimension d, i.e., it has less than d N. R² parameters, where N and R are effective ("average"" mode size and TT-rank respectively.
ma.:=
Linear algebra operations g., element-wise addition, solution of linear systems, convolution, integration, etc.) on tensors the TT-format, respectively, also have complexity linear in dimension if the TT-ranks are bounded. The TT-approximation for tensor discretized multivariable function may be built by efficient numerical methods, e. g., TT-SVD (27). TT-ALS (4), and TT-cross (88
detailed description of the TT-format and related algorithms are given in the works [26]
Below, we discuss only three operations in the TT-format, which will be used later in the work.
Construction of the random TT-tensor. given size (N1, N2,..., S1. Go,..., ., Sa (3-dimensional tensors) with random elements from the uniform distribution on the interval (0. 1). We will refer to this method as tt raraaadnn,,,, (R, [N1, No.... Nal).
order to build random non-negative TT-tensor of Iag...,N, Na) with a constant TT-rank R. it is enough to generate d TT-cores

---

No.
7d
Ra
Ra.
71g
Ns
R1
R
R2
R3
R4.
Ra-1
N2
N1
Na
N2
11
nd-1
N
R. xeMxxxx.... x
92
9s
G
Ga
N₂
N₂
R
Rg
Rg
Ra
Ra
y
=
9
G2
Gs
Ga-1
Ga
Na
N
N
N2
N3
Na-1
N&
Figure 2: Schematic representation of the IT-decomposition. The top picture demonstrates the calculation of the specific tensor element æ [n1, . na] from its TT-representation, and the bottom picture presents the related tensor network diagram.
Computation of the log-likelihood in the TT-format. To calculate the logarithm log P [ar in given multi-index æ = (221, R22, na)), we can use the basic formula 33 and then take the logarithm of the result. It can be shown that this operation has complexity O dd R because, roughly speaking,
we (d 1) times multiply a vector of length R by a matrix of size R x R to get the result. The corresponding method will be called tt log (x).
Sampling from the tensor in the TT-format. We use the approach proposed in the work [10] to generate a multi-index æ with a probability proportional to the corresponding value p = 9 of the TT-tensor P. The method is based the sequential calculation of univariate conditional densities with efficient integration the TT-format. The estimate for its complexity turns out to be
the following: o (K- d (N RR- R+ K d dd-ddoo aNN). where K is number of requested samples, and (x) is complexity of sampling from generalized Bernoulli distribution with outcomes. Note that the algorithm allows sampling in the case of the initially non non-normalized tensor, so we do not have calculate the normalization factor. We will refer to this method tt sample(PP, K).
4 Optimization method PROTES
With the formal scheme of the proposed approach given in Section 2 and the description of operations random, tt log, and tt sample given Section 3 we can formulate our method PROTES for gradient-free discrete optimization in the TT-formar, adam, a function that performs kyd steps of gradient ascent for the TT-tensor at multi-indices X by the well-known Adam method [18). In this case, the learning rate is X, the loss function is given in (2). and tt log with automatic differentiation support is used for the log-likelihood computation.
presented in Algorithm [ We denote as
Computational complexity of the method. Let estimate the complexity of the proposed algorithm, assuming that the number of requests to the target function (black-box) M is fixed. With the known estimate for the complexity of the tt sample function, K K-dd d ((N R) R R+ a(N)). Assuming that the complexity of one gradient step coincides with the complexity of calculating the differentiated function and estimate the total complexity of the tensor
obtain the complexity
of the sampling operations: O M
R
using the estimate for the tt sample function,
updates: 6 K k kind the method
R²") Combining the two above estimates we obtain the complexity of
(0M- dd ( kgd R², +²-00
R-+ (N R+ a(A)))
(4)

---

Algorithm 1 Method PROTES in the TT-format for multidimensional discrete black-box optimization
Data: the function f(2), that computes the value of target tensor 3 E R XM.........>. Na at multi-index at [Ti, -...-.. nal;; the maximum number of requests M; the number of generated samples per iteration K" the number of selected candidates per iteration &; the number of gradient ascent steps Kgdi
the gradient ascent learning rate A; the TT-rank of the probability tensor R.
Result: d-dimensional multi-ine index Zmin. which relates to minimum value of the tensor 'y.
Initialize target -ulti-index and tensor value: Unin None, Youi 00
Generate random non-negative rank- R T'-tensor: T ft_random (R, [N1, Ni,..,,
for iter to M/ K do
Generate K samples from P: 21, .2,,... 2x (ample K)
Compute related lensor values: y1 = f(z) 32 f(22)..... ........................ MK f(xx)
Find indices 8 = {81,.2,... $2,............................. 8x for top-k minimum items in the list [y1, 191,32.-.. va.....yK.. yx]
Collect k-top candidates: X { 4Pag..... ................................ ). 3 (Bla2 Mag..... May If 3 contains value less than Ymen then update Imii and yman
Update the TT-tensor: 3 adam(9, T,K, hig, A) // with the loss function and the method tt log
10 return Pmin-
However, it is important to note that this estimate does not consider the complexity of calculating M times the objective function f, which in practical applications can be significant and many times greater than the estimate (4)
The intuition behind the method. The proposed method PROTES, like most gradient-free op- timization approaches, is empirical; however, we can establish its connection with a well-known REINFORCE trick algorithm (43) Let make monotonic transformation (ar of the target func- tion f to be minimized that transtorms minimum to maximum. A reasonable choice for F(.) is the Fermi-Dirac function
FIf (x) exp((f(x))
Ymin E)/T) +1'
where 3min is an exact or approximate minimum of > 0 a parameter and E is some small threshold. With the function F(f) we can find maximum of the expectation maxØ Eg, Fif)(Eo), where a family of random variables &o has parameterised distribution function Pg (z). Using REINFORCE trick, we can estimate the gradient of the expectation by the following Konte-Carloy expression
M
VeE ~
Fff((.)Ve) log (x:))
(5)
M
where {x are independent realizations of the random variable Eg. If we find the optimal values of 0. then we expect the optimal distribution P, to peak at the maximum point for function F(f). Thus, we can obtain the argument of its maximum by sampling from this distribution. For very small values of T, only few terms contribute to the sum (5), namely those æ; for which f(ax,) Yman E is hold. For these values of F(f) is close to 1, while for the other samples its value is 0. Hence, we can discard all other samples and keep few with the best values. So, we come to the loss function (2). where instead of the parameter E we use fixed number l of the best samples, i. e., the samples for which the value of the target function f is the smallest.
Application of the method to constrained optimization. A very nice property of the proposed method is that it can be adapted to efficiently handle constraints such specified set of admissible multi-indices. One option is just to remove invalid samples from the top-k values, but in some cases, the probability of sampling multi-indices that are admissible very low, so this approach will not work. Instead, if the constraint permits, we use the algorithm from the work (32) for the constructive building of tensors the TT-format by a known analytic function, which defines the constraints. Once the indicator tensor (1 if the index is admissible and 0 if it is not) is built in the TT-format, we can just initialize the starting distribution by it, and it will be guaranteed that the samples almost always belong to the admissible set.
5

---

5 Related work
Below, we briefty analyze classical approaches for discrete optimization and then discuss the methods based on the low-rank tensor approximations, which have become popular in recent years.
Classical methods for gradient-free optimization. In many situations, the problem--peciffe target function is not differentiable, too complex, its gradients are not helpful due to the non-convex nature of the problem [T]. and standard well-known gradient-basee methods cannot be applied directly. selection, training neural networks with discrete weights,
The examples include and policy optimization in reinforcement learning. In all these contexts, efficient direct gradient- free optimization procedures are highly needed. As for high- nigh-dimensional black-box optimization, evolutionary strategies (ES) [9] are one of the most advanced methods. This approach aims to optimize the parameters of the search distribution, typically a multidimensional Gaussian, maximize the objective function. Finite difference schemes are commonly used to approximate gradients of the search distribution. Numerous works proposed techniques improve the convergence of ES [23] for example, second-order natural gradient (42) or the history of recent updates (Covariance Matrix Adaptation Evolution Strategy; CMA-ES) [15] may be used to generate updates. There is also large variety of other heuristic methods for finding the global extremum. In particular, we note such popular approaches as NoisyBandit [33) Particle Swarm Optimization (PSO) [7], Simultaneous Perturbation
Stochastic Approximation (SPSA) 22 Differential Evolution (DE)
38 and scrambled Haabbedddd-mmrrleey
(scr- Hoemersmersly)) (14).
Tensor-based methods for gradient-free optimization. Recently, the IT-decomposition has been actively used for multidimensional optimization. An iterative method TTOpt based on the maximum volume approach is proposed in the work (37). TTOpt utilizes the theorem of sufficient proximity of the maximum modulo element of the submatrix having the maximum modulus of the determinant to the maximum modulo element of the tensor. Based on this observation, tensor elements are sampled from specially selected successive unfoldings of the tensor. Dynamic mapping of the tensor elements is carried out to find the minimum element, which converts the minimum values into maximum ones. The authors applied this approach to the problem of optimizing the weights of neural networks in the framework of reinforcement learning problems in (37 and the QUBO problem in (24). A similar optimization approach was also considered in 33 and [35]. One more promising algorithm, Opt was proposed in recent work (5). This approach based on the probabilistic sampling from the TT-tensor and makes it possible to obtain a very accurate approximation for the optimum of the given TT-tensr.. However, this method intended for directly optimizing the TT-tensors, which means that its success strongly depends on the quality of the TT-approximation for the original multidimensional data array. Therefore, one of the related methods in the TT-format (TT-SVD, TT-ALS, TT-cross, etc.) should be additionally used for approximation.
6 Numerical experiments
To evaluate the effectiveness of the proposed method, we carried out a series of 20 numerical experiments for various formulations of model problems. The results are presented in Table [ where we report the approximation to the minimum value for each model problem (P-1, P-2,
P-20)
and all considered optimization methods (PROTES, BS- BS-2, ... BS-7). Taking into account the analysis of discrete optimization methods in the previous section as baselines we consider two tensor-based optimization methods: TTOp (BSI) and Optima-T TT (BS2), and five popular gradient- free optimization algorithms from the nevergrad framework (3)
T OnePlusOne (BS3), PSO (BS4), NoisyBandit (BS5), SPSA (BS6), and Portfolio approach (BS7), which based on the combination of CMA-ES, DE, and scr- sar-Hemerrsley methods. The model problems and obtained results will be discussed in detail below in this section.
For all the considered optimization problems, we used the default set of parameters for baselines, and for PROTES we fixed parameters K 100, k 10, lig 1, X 0.05. R (the description
'ee used implementation of the method from 'We used implementation from https: /tt::://ggttttttte com http://hhtt..ooooocccetttttttttr The TT-tensor for optimization was generated by the TT-cross method. "ee h thub. http:///...ccbbboeeeerrr//eergge

---

Table 1: Minimization results for all selected benchmarks 1:01-PM (P-1),-))), The values obtained by the proposed method PROTES and by all considered baselines (BS1 BS7) are reported.
PROTES
BS-1
BS-2
BS-3
B8-4
BS-5
BS-6
BS-7
P-01
1.3++01
1,3++01
1.3++01
1,3+001
1.3+001
2.18+01
1.8+00
1.38+01
P-02
6.5r-00
6.5
6.5r.400
6.9n-00
6.8n400
1.5:401
7.5:400
6.80400
P-03
-9.4:-01
-9.4..01
-9.4.-01
-9.x-01
-9.4x-01
-355-00
9.44-01
-9.4u-01
P-04
1.3+-00
1.35+00
1.3r-00
1,32+00
1.32+00
6.36+00
1,30+00
1.36+00
ANALYTIC
P-05
+3.7:400
43.7.400
-3.7r.00
-2600000
3.0m-00
-1.8::400
-1,2:400
3.710400
FUNCTIONS
4.28:01
P-07
6.21-006
2:2
3.12408
6,05+01
6.0E-01
6.0x+01
6.0E+01
1.25+02
1.0+002
6.06+01
P-09
2.7:400
2.7r400
2.71.400
3.0m+00
2.7x400
2,00000
5.40400
2.7x400
P-10
-8.7E-02
+8.7+002
-8.7+00
6. 6.16e02
6.9E+02
7.08+02
2,68+03
8.58+02
P-11
-3.6r-02
-3.54403
4m-02
2r-02
-3.40022
A.35402
3.35402
-3.6::02
P-12
-5.9
-5.9.403
+33
5.6eed8
5.9xe03
6.9
QUBO
P-13
-3.15+00
-3,6+00
2 268000
0.0E+00
1.5+01
2.85+02
-296+00
-3,0+00
P-14
-3.1r403
-2. 2.88
-3 3.0r.403
2.66403
-3 3.0n+03
-2.71:403
-300:403
-3.08403
P-15
6.7r-03
7.4p-03
2.3r-02
8.4n-03
8.9n-03
3.1n-02
8.7n-02
7,3r-03
CONTROL
P-16
1.4r-02
2. 66-002
3.56-02
1.7p-02
1,7p-02
5,36-02
5,26-02
1.4r-02
P-17
3.0r-02
5.76-01
1.5u-01
4.8u-02
3.6n-02
7.7t.02
5,36-02
3.76-02
P-18
1 41-002
1.1r.02
1.40-02
3.4n-02
6.2n-02
2.8t:01
6.41-02
2.11-02
CONTROL
P-19
6.4E-02
5.7E:01
6.7p: 02
FAIL
FAIL
FAIL
FAIL
FAIL
+CONSTR.
P-20
1,5r-01
FAm.
2.0 0c-01
FAIL.
FAIL
FAIL
FARL.
FAR
of these parameters was presented in Algorithm Dp. requests to the objective function was fixed at the value M = 10". As seen from Table I PROTES, contrast to alternative approaches, gives a consistently top result for almost all model problems (the best result for 19 of the 20 problems considered).
For all methods, the limit on the number of
6.1 Multivariable analytic functions
First, we consider the optimization task for various tensors arising from the discretization of mul- tivariable analytic functions. We select 10 popular benchmarks: Ackley (P-01), Alpine (P-02), Exponential (P-03), Griewank (P-04), Michalewiez (P-05), (P-66)) Qing (P-07), Rastrigin (P-08), Schaffer (P-09) and Schwefel (P-10). These functions have a complex landscape and are often used in problems to evaluate the effectiveness of optimization algorithms [8 [6]. including tensor-based optimizers 39). We consider the 7-dimensional case (since this is the dimension of the Piston function) and discretization on a uniform grid with 16 nodes.
As follows from Table [ (benchmarks P-1, P-2, P-10), our method, like the other two tensor approaches (BS-1 and BS-2), gave the most accurate solution for all model problems. The most sophisticated approach from the nevergrad package (BS-7) was the next in accuracy (the method did not converge only in two cases out of ten).
6.2 Quadratic unconstrained binary optimization
QUBO is a widely known NP-hard problem 12] which unifies a wide variety of combinatorial optimization problems from finance and economics applications to machine learning and quantum computing. QUBO formulation in a very natural manner utilizes penalty functions, yielding exact model representations in contrast to the approximate representations produced by customary uses of penalty functions. The standard QUBO problem can be formulated as follows
(z)) T Tx min, m
s.l. x E (0,
where x is a vector of binary decision variables of the length and Q E RRdx is a square matrix of constants. In all our experiments, we fixed the number of dimensions
50.
'This function corresponds to the problem of modeling the time that takes a piston to complete one cycle within a cylinder; the description of its parameters can be found in 44

---

Table 2: Minimization results for benchmark P-12 for the 1000 dimensional case with different random initializations. Results for each value of p are averaged over 5 runs of optimizers.
PROTES
BS-1
BS-2
BS-3
BS-4
BS-5
BS-6
BS-7
P 0.1
-4,868+05
+.,00+05
FAIL
-4.10
44 4.356=05
-3,94E+05
-3.44++05
+..7
p 0.2
-9.70x+05
-8.88:405
FAII.
-8 24m-05
-8. 71n-0
-7,87m405
-7.644405
-9.38:405
P-12
p 0.3
-1.462+06
-1.322406
PAIL
-1.23e-06
-1.30c006
-1.182006
-1.44206
-1.412105
p = 0.4
-1.94x+06
-1,758+06
FAIL
-1,44-066
-1.726-06
-1.577+06
-1.528406
-1.888+06
P =
-2.43x446
-2.1914466
FAIL.
-2.044-06
-2. -2.166406
-1.961466
-1.900406
-2.35::06
We consider the following QUBO problems from the qubogen package Max-Cut Problem (P-11: which refers to finding partition of an undirected graph into sets such that the number of edges between the two sets is as large as possible), Minimum Vertex Cover Problem (P-12: which refers to finding a cover with minimum number of vertices in the subset of the graph vertices such that each edge in the graph is incident) and Quadratic Knapsack Problem (P-13; which refers to finding a subset of maximum profit that satisfies the budget limitations from set of potential projects with specified interactions between pairs of projects).
We also consider one more benchmark (P-14) from the work (IT) (problem K3; d 50), where angle-modulated bat algorithm (AMBA) was proposed for high-dimensional QUBO problems with engineering application to antenna topology optimization. This is the ordinary binary knapsack problem with fixed weights vc; E [5, 20). profits E [50, 100) (i = 1.2...., d), and the maximum capacity C' 1000. In experiments, we used the same values of the weights and profits as in (TT..
For all four considered problems (P-11, P-12, P-13, P-14), the proposed method PROTES gives the best result, as seen from Table and the baseline BS-7 again turned out to be the next in accuracy. We also note that several optimization methods were compared in [11] for the P-14 problem: BPSO (with the result- - 2854),, BBA (with the result- 2976), AMBA (with the result -2956), A-AMBA (with the result- - 2961), P-AMBA (with the result 2989),, and the solution obtained using the PROTES method (the result- 3079) turns out to be more accurate.
To study the stability of the proposed optimization method in the essentially multidimensional case, also consider d' = 1000 for the benchmark 12. Additionally, vary the value of the vertex connection probability in the generated random graph (i.c., the parameter in the corresponding function from the well-known networkx package), getting 5 different optimization problems with p = 0.1, p = 0.2, p 0.3, p 0.4, and p = 0.5 (note that in Table T we used the value = 0.5). We performed the calculations with PROTES and baselines BS1 BS7 with the same settings as above. In Table 2 we report the results averaged over 5 runs of each optimizer for each value of p. In all cases, PROTES& demonstrates the best result. We note that BS2 failed for this problem because the selected dimension value is too large to construct the TT-approximation of the tensor.
In Figure 3 we plotted the convergence curves for p 0.1, p 0.3, p 0.5, and all optimizers (except the BS2, since it failed), showing the average result and the spread of values over restarts for 1- (d = 1000). Our optimizer PROTES has slight variance of the solution and consistently shows the best result when the number of iterations exceeds 2000.
6.3 Optimal control
Suppose have state variable R controlled by binary variable z called control (i.e., it is
just a switch with modes "off" = O and "on" = 1) over some discrete interval of time [0, 7). The state z(1 1) at time / depends the control stt) time t and is obtained from the solution of the following differential equation å(r) g(=(T), (t)), S < f + 1, where the function is called an equation function. The optimal control problem is to find such sequence of controls x² [x* (0), x=' (1), ..... x² (T') (optimal solution) over the given time interval [0, 2'] that minimizes the given objective function G.
See https:://iith// httpg/ttwwb.comtmmuueee/uuuuoeee
8

---

BS-1
PROTES
BS-1
-17
D8-1
PROTES
BS-3
-1.1
-18
BS
BS
BS
-1.2
RS-5
RS-5
B8-6
B8/
2.1
-425
-22
-4.00
-14
4.75
-44
-500 10'
-1.5
-25
16
Number requests
Number of requests
Number requeets
(a) =
0.1 (scaled x 10
(b) = 0.3 (scaled x 10
(c) p 0.5 (scaled x 10
Figure 3: Minimization results benchmark P-12 the 1000-dimensional case. The results for
various values of the vertex connection probability p of the generated random graph are presented. For each of the optimizers PROTES, BSI, and BS3 BS7 (BS2 failed for this problem), we plot the value of the solution averaged over 5 runs with solid line and fill in the area between the worst and best result with the same color.
Formulating the problem mathematically, we need to find such a solution
2(0) =
G(2, x)
min, 2,
s.t.
((T) g(=(r) sxt), 15t<t+1, xxt) E (0, 1]. t= 0..,,
where = |2(0), 21).....,, xT)) is a state variable path. numerical experiments, we consider the nonlinear equation function g(2,x)) 23 x, and since it is nonlinear, finding an optimal solution
T
raises a lot of difficulties. We take the objective function G in the form G(z, a)
Σ (=(1)= -(e)))=rer)) = 1=0
The initial and the reference state are fixed at values 2g = 0.8, 2 = 0.7. For a fixed initial value 20 and fixed equation function g, the objective function G can be represented as binary multidimensional tensor, whose elements are calculated using the following function: f(x:) = G(z(xz), 2;, hence the optimal solution x=
apply discrete optimization methods to find Emin, which approximates
We considered several values for variable T, such 25, 50. and 100 (benchmarks P-15, P-16, and P-17, respectively). As follows from the results presented Table PROTES gives the most accurate solution in all three cases. Note that result comparable in accuracy to our method is obtained only in one case when using basclines (i.c., P-16, BS-7).
6.4 Optimal control with constraints
In practical applications, some conditions constraints may be imposed on the solution of the optimal control problem. In this work, we consider the following control constraint P: the control variable {0, 1 can take the value "I" less than 3 times in row during the whole time interval. Formally, this can be written as follows:
x f 2 x
z|t =
P
V1 2, we let lct 0 for t and N/
ze at - -/]--]]
To account for this condition in the PROTES algorithm, we constructively build the initial distribution in the form of an indicator tensor as described in Section n the constrained optimization subsection.
The details of this construction are presented in the Appendix. The numerical results for = 25 (P-18), T = 50 (P-19), and T 100 (P-20) are reported in Table In two cases out of three (P-19, P-20), our method showed the best result, and in one case (P-18) slightly yielding to the TTOpt method (BS-1) 1), which, however, in two other cases, gave significantly worse result.

---

TT-rank (R) 3
TT-rank (R)
TT-rank (R) =
10
10
Beich size (K) - 50
Batch size (K) - 100
Batch size = 200
Betch size (K) = 250
10
10"
10
10"
10"9
10"
20
25
10
10 15 Top candidates (k)
20
Top candidates (k)
Top candidates (K)
Relative error of the optimization result with PROTES method for the P-14 model problem with known exact minimum for different values of the hyperparameters K, k and R.
Figure 4:
6.5 Robustness and performance of the PROTES
The results in Table D relate to the "intuitive" selection of the hyperparameters for the PROTES method (as mentioned above, have used the values: K 100, k = 10, kig 1, X 0.05, and R = 5). In Figure 4 we present an analysis of the dependence of the optimization result for the benchmark P-14 the choice of hyperparameters K, k, and R. with fixed kyd 1 and 0.05. We report the relative error of the result for all combinations K = 50, 100, 150, 200, 250; 5, 10, 15, 20, 25; R 3,5,7 As main calculations are not optimal for this particular problem, that is, additional fine-tuning of the method for specific problems classes of problems is possible. At the same time, according to the results in Figure 4 the method remains stable over wide range of byperarrameter values. We also note that all computations were carried out a regular laptop, while the operating time of the considered optimizers was commensurate; for example, for the benchmark P-17, the measured operating time (in seconds) is: PROTES (641), BS-1 (607), BS-2 (4245), BS-7 (780).
see from the plots, the hyperparameters used in the
more detailed analysis of the PROTES performance and the dependence of optimization results on presented in the Appendix. Also, in the Appendix, consider
the values of hyperparameters another promising application of our optimizer for several popular reinforcement learning problems from Mujoco/ / OpenAl-GYM collection.
7 Conclusions
In this work, we presented an optimization algorithm PROTES based on sampling from the probability density defined the tensor train format. We used the same set hyperparameters for all considered numerical experiments, so our algorithm is rather universal To take into account the constraints, as in the problem of optimal control with constraints, we only considered them in the form of specially selected initial approximation (a special form of an indicator tensor in the tensor train format); further on, the algorithm did not consider the constraints explicitly. This approach allows us to extend the algorithm's 's capabilities by using the properties of the tensor train representation. Numerical experiments show that we outperform many popular optimization methods.
The main direction in our future work is scaling the method to large dimensions. For d 2 1000, we have encountered numerous technical difficulties, which can be alleviated by other tensor formats (such as hierarchical Tucker, which can be parallelized over d) and more efficient implementations of the optimization method (now we use standard automatic differentiation without special tensor optimization methods such as Riemannian optimization).
Acknowledgements
The work is supported by the Ministry of Science and Higher Education the Russian Federation (Grant 075-10-2021-066
10

---

References
[1] Stephane Alarie et al. "Two decades of blackbox optimization applications". In:
EURO Journal Computational Optimization 9 (2021), p. 100011.
[2] Boian Alexandrov al. "Challenging the of dimensionality in multidimensional merical integration by using a low-rank tensor-tramt train format". In: Mathematics 11.3 (2023),
p. 534.
[3] Pauline Bennet et al. "Nevergrad: black-box optimization platform". In: SIGEVOlution 14.1
(2021), pp. 8-15.
[4] Andrei Chertkov, Gleb Ryzhakov, and Ivan Oseledets. "Black box approxi- mation in the tensor train format initialized by ANOVA decomposition". In: SIAM Journal Scientific Computing 45.4 (2023), A2101-A2118.
[5] Andrei Chertkov et al. "Optimization of functions given in the tensor train format". In: arXiv preprint arXiv:2209. 14808 (2022).
[6] Andrzej Cichocki et al. "Tensor networks for dimensionality reduction low-rank tensor decompositions.. In: Foundations and Trends in Machine Learning 9.4-5 (2016), pp. 249-429.
and large-scale optimization: Part
[7] Andrzej Cichocki et al. "Tensor networks for dimensionality reduction and large-scale optimization: Part 2 applications and future perspectives". In: Foundations and Trends in Machine Learning 9.6 (2017), pp. 431-673.
[8] Johannes Dieterich and Bernd Hartke. "Empirical review standard benchmark functions using evolutionary global optimization", In: Applied Mathematics 3.10 (2012), pp. 1552- 1564.
[9] Benjamin Doerr and Frank Neumann. "A survey on recent progress in the theory of evolutionary algorithms for discrete optimization". In: ACM Transactions on Evolutionary Learning and Optimization 1.4 (2021), pp. 1-43.
[10] Sergey Dolgov et al. "Approximation and sampling of multivariate probability distributions in the tensor train decomposition". In: Statistics and Computing 30 (2020), pp. 603-625.
[11] Jian Dong. Zhiyu Wang, and Jinjun Mo. "A phase angle-mon modulaedd bat algorithm with appli- cation to antenna topology optimization". In: Applied Sciences 11.5 (2021), p. 2243.
[12] Fred Glover et "Quantum bridge analytics a tutorial on formulating and using QUBO models". In: Annals of Operations Research (2022), pp. 1-43.
[13] Lars Grasedyck, Daniel Kressner, and Christine Tobler. "A literature survey of low-rank tensor approximation techniques". In: GAMM-MMiteillngenn 36.1 (2013), pp. 53-78.
[14] John Hammersley. "Monte Carlo methods for solving multivariable problems". In: Annals of the New York Academy of Sciences 86.3 (1960), pp. 844-874.
[15] Nikolaus Hansen. "The CMA Evolution Strategy: Towards new evolutionary computation: advances in the estimation of distribution algorithms. Springer Berlin Heidelberg, 2006, pp. 75-102.
Comparing Review". In:
[16] Momin Jamil and Xin-She Yang. "A literature survey of benchmark functions for global optimization problems". In: Journal of Mathematical Modelling and Numerical Optimisation 4.2 (2013), pp. 150-194.
[17] James Kennedy and Russell Eberhart. "Particle swarm optimization". In: Proceedings of ICNNICCN..... 95-interrnation conference on neural networks. Vol. 4. IEEE. 1995, pp. 1942-1948.
[18] Diederik P Kingma and Jimmy Ba. "Adam: a method for stochastic optimization". In: arxiv preprint arviv: urXiv:1412.668 141..6922988 (2014).
[19] Kirandeep Kour et al. "Efficient structure-preservingg support tensor train machine". In: Journal of Machine Learning Research 24.4 (2023), pp. 1-22.
[20] Caio Kalil Lauand and Sean Meyn. "Approaching quartic convergence rates for quasi-stochastic approximation with application to gradient-free optimization". In: Advances Neural Information Processing Systems. 2022.
[21] Tianyi Lin, Zeyu Zheng. and Michael Jordan. "Gradient-free methods for deterministic and stochastic nonsmooth nonconvex optimization". In: Advances in Neural Information Processing Systems 35 (2022), pp. 26160-26175.

---

[22] John Maryak and tion by simultaneous Proceedings of the 2001 American Control Conference.(Cal No. 01CH37148).
Daniel perturbation stochastic
Chin.
"Global
random approximation".
optimiza-
In:
Vol. 2.
IEEE. 2001, pp. 756-762.
[23]
Yurii Nesterov and Vladimir Spokoiny. "Random gradient- minimization of convex func- tions". In: Foundations of Computational Mathematics 17.2 (2017), pp. 527-566.
Artyom Nikitin et al. " "Are quantum computers practical yet? case for feature selection in recommender systems using tensor networks". In: arXiv preprint arXiv:2205.08 (i4::22000044900 044900 (2022). [25] Georgii Novikov, Maxim Panov, and Ivan Oseledets. "Tensor-train density estimation", In: Uncertainty in artificial intelligence. PMLR. 2021, pp. 1321-1331.
[24]
[26] Ivan Oseledets. (Tensor-train decomposition". In: SIAM Journal on Scientific Computing 33.5
(2011), pp. 2295-2317.
[27] Ivan Oseledets and Eugene Tyrtyshnikov. "Breaking the curse of dimensionality, or how to use SVD in many dimensions". In: SIAM Journal on Scientific Computing 31.5 (2009), pp. 3744- 3759.
[28] Ivan Oseledets and Eugene Tyrtyshnikov. "TT-sross approximation for multidimensional arryys". In: Linear Algebra and its Applications 432.1 (2010), pp. 70-88.
[29] Olivier Parcollet et al. "Learning Feynman diagrams with tensor trains". In: Bulletin of the American Physical Society (2023).
[30]
Gary Parker and Ronald Rardin. Discrete optimization. Elsevier, 2014.
[31]
Jun Qi et al. "Exploiting low-rank tensor-train deep neural networks based on Riemannian gradient descent with illustrations speech processing". In:
IEEE/ACM Transactions Audio, Speech, and Language Processing 31 (2023), pp. 633-
642.
[32] Gleb Ryzhakov and Ivan Oseledets. "Constructive TT-representaation of the tensors given as index interaction functions with applications". In:
11th International Conference on Learning Representations, ICLR. 2023.
[33]
Jonathan Scarlett, Ilija Bogunovic, and Volkan Cevher. "Lower bounds regret for noisy gaussian process bandit optimization". In: Conference on Learning Theory. PMLR. 2017, pp. 1723-1745
[34] Cheryl Selvanayagam et "Global optimization of surface warpage for inverse design of ultra-thin electronic packages using tensor train decomposition". In: IEEE Access 10 (2022),
pp. 48589-48602.
[35] Suhan Shetty et al. "Tensor train for global optimization problems robotics". In: arXiv preprint S.iv:2206.050777 (2022).
[36] Micheline Soley al. "Global optimization with the iterative power algorithm via quantum computing and quantics tensor trains". In: Bulletin of the American Physical Society (2023).
[37] Konstantin Sozykin et al. "TTOpt: train-based optimization and its application to reinforcement learning". In: Advances in Neural Information Processing Systems. 2022.
maximum volume quantized tensor
[38] Rainer Storn and Kenneth Price. "Differential evolution a simple and efficient heuristic for global optimization over continuous spaces". . In: Journal Global Optimization 11.4 (1997), pp. 341-359.
[39] Christoph Strössner, Bonan Sun, and Daniel Kressner. "Approximation in the extended func- tional tensor train format". In: arXiv preprint XrXiv:22 1. 1338 (2022).
[40] Minghua Wang et al. "Tensor decompositions for hyperspectral data processing in remote sens- ing: a comprehensive review". In: IEEE Geoscience and Remote Sensing Magazine (2023).
[41] Xiaoxing Wang et al. "ZARTS: on zero-order optimization for neural architecture search". In: Advances Neural Information Processing Systems 35 (2022), pp.
[42] Daan Wierstra et al. "Natural evolution strategies". In: Journal of Machine Learning Research
15.27 (2014), pp. 949-980.
[43] Ronald Williams. "Simple statistical folllwwig algorithms for connectionist reinforce- ment learning". In: Machine learning 8.3-4 (1992). pp. 229-256.
[44] Vitaly Zankin, Gleb Ryzhakov, and Ivan Oseledets. "Gradient descent- D-optimal design for the least-squares polynomial approximation". In: arXiv preprint
(2018).
12

---

Supplementary Material
A1 Stability of the PROTES
To check the stability of the optimization result, series of 10 calculations were performed for each method (PROTES, BS1 BS7) with random initializations. We consider the binary knapsack problem from the work [1] (benchmark P-14, described in detail in the main text of the work) with the known exact minimum 3103. For all methods, the limit on the number of requests was fixed at the value M = 10 All other parameters were the same as in the computations from the main text, i.e., the PROTES parameters are K = 100, & = 10, kgd = 1, À = 0.05, R 5. In Table we present the average and best results over 10 runs for each optimization method. As follows from the reported results, only PROTES and Portfolio method (BS-7) managed to successfully find the exact optimum, while the average result for PROTES is significantly better than that of all the baselines.
Table 3: Average and best result for 10 independent runs for the P-14 benchmark.
PROTES
BS-2
BS-3
BS-5
BS-6
BS-7
MEAN
-3095
2992
-3048
-2650
-2937
-2701
-3064
-3075
P-14
BEST
-3103
33044
-3084
-2825
-2996
-2752
-394
-3103
A2 Reinforcement learning with PROTES
As in the work (37], we consider maximization of the cumulative reward of the episode for several popular reinforcement learning problems from Mujoco/ OpenAl-GYM collection: InvertedPendulum (P-21), Swimmer (P-22), Lunar Lander (P-23), and Half Cheetah (P-24). The policy represented by neural network with three layers and tanh activations. Parameters of the neural network discretized on a grid with limits from -1 to +1 and 256 nodes (i.c., mode size N). The total number of neural network parameters for the problems under consideration turns out to be from 50 to 68 (i.e., dimension d). Thus, as in the work (37). we obtain a discrete on-policy learning problem (search for the values of the neural network parameters that lead the maximum reward). To speed up the calculations, we took only two baselines (BS1, and BS7) for comparison, which showed the best results in the main experiments presented in the paper. We used the same optimizer settings in the main calculations described in the text of the paper; however, in this case, we chose the budget value equal to M 10 in (37). For PROTES, BSI, and BS7, we repeated the calculations three times and reported the average result for found maximum and computation time in Table 4 PROTES shows the best result in terms of accuracy in 3 cases out of 4, and once it in second place, slightly inferior to the BS7.
Table 4: Maximization results for InvertedPendulum (P-21), Swimmer (P-22), Lunar Lander (P-23) and Half Cheetah (P-24) benchmarks from Mujoco collection. We report the reward and calculation time (seconds), averaged over 3 runs, for the proposed method PROTES and baselines BS1 and BS7.
REWARD
TIME
PROTES
BS-1
BS-7
PROTES
BS-1
BS-7
P-21
1.00k+03
7.28u-02
3.492402
8.1be02
6.3u102
1.3
P-22
3.618+02
3.436-02
3,55E+02
9.76+03
7.7E+03
2.36+04
RL AGENTS
P-23
2.800x422
-1 -1.874401
2.64::402
3,50403
8.1m403
1.51:44
P-24
2.26ex03
7.78e-02
2.662+03
8.5ee03
5.5er03
2.28404
A3 Choice of hyperparameters for the PROTES
We conduct additional experiments with varying the values of the hyperparameters of the PROTES method for benchmarks P-01 (Figure and Figure we report the optimization result) and P-14 (Figure 4 in the main text and Figure 7 below; we report the relative error of the optimization result). In the first series of experiments (Figure (Figure3) below and Figure 4 in the main text), we fixed the values kgd- = 1 and X 0.05 and tried all combinations of the remaining hyperparameters:
13

---

TT-rank
TT-rank
TT-rank (R)
10
8406
X1mml
8200
<00
8.300
7.8
7.8
7.0
7.6-
7.4
(K)
7.4
7.4
Batch 100
Batch 150
72
Batch vice (K) 200
7.2
722
Batch size (K) 250
7.0
7.0
15
20
10
15
Top (k)
Top candidates (k)
Top candicites
Figure 5: Optimization result (approximation of the minimum value) with PROTES method for the P-01 for different values of the hyperparameters: K, and
TT-aakk (R)
TT-rank (R) =
TT-rank (R) =
15
15-
03
LR 10c 02
5.e-02
15
1.0-21
13-
13
-..
LR 5.0e-
11
11-
04
0-
1 10
25
30
200
1 30
25
50
100
10
25
50
100
Gradient (k-ge)
Gradient stops (k-gd)
Gracient steps (k.gd)
Figure 6: Optimization result (approximation of the minimum value) with PROTES method for the
P-01 model problem for different values of the hyperparameters: learning rate (LR), kyd. R.
(R)
(R)
10
10"
10"
LR
LR 1.0
UR = 5 5.0e-07
RR - 1.0e-01
5.0e-01
10
10
10 25
5D
100
1 10
25
50
100
1
10
25
50
100
Gradient steps (K_gd)
Gradient steps (k_gd)
Gradient steps (K_gd)
Figure 7: Relative error of the optimization result with PROTES method for the P-14 model problem with known exact minimum for different values of the hyperparameters: learning rate (LR), kyd, R.
K = 50, 100, 150, 200, 250; = 3, 10, 15, R = 3,5,77 the second series of experiments (Figure 6 and 7) we fixed the values K 100 and k = 10, and tried all combinations: kgd 1, 10, 25, 50, 100, = 0.005, 0.01, 0.05, 0.5, R 3,
As follows from the presented results, for problem P-01, the dependence the choice of hyper- parameters turns out to be extremely weak, except for outliers for the small values of the learning rate À 0.005 and À 0.01. For more complex problem P-14, the dependence of the result on the
14

---

Table 5: Computation time in seconds for all selected benchmarks (P-01 - P-20) and for all used optimization methods (PROTES, BS1 BS7).
PROTES
BS-1
B5-2
BS-3
BS-4
BS-5
BS-6
BS-7
P-01
3.28
0.06
0.11
23.45
25.31
22.78
22.03
67.23
P-02
2.25
0.05
0.06
22.42
25.68 68
20.6
19.78
77.01
P-03
2,36
0.05
0.03
26.11
23.62
20.55
19.71
65.62
2.32
0.05
0.07
22.04
24.04
20.67
19.9
65.31
ANALYTIC FUNCTIONS
0.05
26.47
20.78
21.58
26.12
21.66
21.26
2.31
0.06
P-10
2.35
0.05
0.03
23.56
33.66
21.25
22.61
0.06
17
17.98
0.05
16.77
QUBO
P-13
2.29
0.05
0.34
18.78
21.09
20.51
17.97
77.22
P-14
2.38
00.13
0.4
16.41
27.42
17.61
17.59
74.77
P-15
513.6
1256 0
1839 0
550.0
545.0
620.3
$30.4
707.5
CONTROL
P-16
542.4
969.4
3007.0 0
578.0
595.0
570.7
573.7
697.5
P-17
640.7
607..
4245.0 0
673.6
687.66
661.0
687.1
779.5
P-18
528.1
92.9
588.1
319.1
516.2
202.9
533.9
474.5
CONTROL
P-19
8.74
69.49
912.0
17.13
17.45
16.16
16.84
47.43
+CONSTR.
P-20
9.23
0.53
931.8
21.16
21.42
20.28
20.84
62.26
choice of hyperparameters is more complex, but it can be seen that the method remains stable wide range of hyperparameter values.
A4 Performance comparison of the optimization methods
All the results described in the main text and presented there in Table 1 were obtained on regular laptop. Table report the related computation time for each method (PROTES, B$1 BS7) and each model problem (PO1 P20). As follows from the results, the PROTES works much faster than classical optimization methods (BS3 BS7), as well as faster than tensor-based methods (BS1 and BS2) for several optimal control problems (P-15 P-20). However, the TTOpt (BS-1) and Optima-TT (BS-2) methods are faster for simpler analytic (P-01 P-10) and QUBO (P-11 - P- 14) problems. This is due to the fact that within the framework of the dynamic TT-rank refinement procedure in these methods, the TT-rank, and hence the computation time, turn out to be significantly higher for the optimal control problems. We note that the very short running time of the TTOpt method for P-20 because most the requests the method did not satisfy the imposed constraints, and in this case, the differential equation not solved.
A5 Derivative functions in the optimal control problem
We use the following derivative functions for the constructive building of the tensor described in
0, # or # otherwise,
f²(2) = None,
(+()) min (1, + 1),
for all TT-cores except the last one 1, 2, ... d 1), and for the last TT-core
(d(x)
1,
x =
0,
otherwise,
1,
2 1
fl(x)
otherwise.
tensor in the TT-format built on such derivative functions equal if there are less than ones among its vector argument in row, and is equal to 1 in all other cases.
15

---

Let us briefly explain why such derivative functions give tensor of the restriction condition (i.e., the constraint in the considered optimal control problem). Recall that the upper index in the derivative function notation corresponds to the index number of the tensor argument, and the lower index corresponds to the value of this argument index. our scheme, the argument of the derivative functions has the meaning of the number of ones that already stand to the left of the given index.
First, let focus on the functions for all indices except the rightmost one =, ......, d If the current index is one, then we simply increase the value of the argument by 1 and pass it on, to the input of the next derivative function. This is what the function f does (we take the maximum, so if the number of ones has already reached 1, we do not care if it is greater than or equal to I; note that the maximum operation reduces the TT-rank but does not affect the result). For the zero value of the current index, if the argument is zero (which means that the previous index is also zero), the function fö also returns zero as nothing have changed. the argument of this function is I, means that there are or more ones in a row to the left of the considered index, so the given condition is not violated, and the function simply returns zero, which means that there are no ones. If the argument greater than zero but less than I, means that the condition is violated because the previous sequence of ones of length less than cut off at the current index. this case, the function fö returns None, which means that the value of the tensor will be 0 regardless of the subsequent indices.
The functions corresponding to the last (k d) index behave similarly. Namely, if the last index is zero, the function fo returns 1 if its argument is O or since it does not violate the condition as just explained. Otherwise, it returns O. If the last index is 1, the function f returns only if its argument equal (I 1) or since this means that the current one an element of sequence of ones of length least I. Otherwise, the function fl returns 0.
16


Reviews:
[{"review_id": "WTfOuRWxG6", "reviewer": "Reviewer_7NgC", "summary": "The paper targets the challenging gradient-free discrete multidimensional optimization. It is formalized as finding a vector of d indices whose associated entry in a d-dimensional tensor is minimized, where the tensor contains all the possible values for the black-box objective function. There are an exponential number of entries in the tensor for brute-force search. Instead of enumerating all of them, the paper proposes to optimize a low-rank tensor whose each entry represents the probability of selecting the corresponding entry as the minimum in the objective value tensor. The algorithm starts from a random probability tensor and updates it using a similar trick as REINFORCE algorithm. In experiments on benchmark tasks, it performs comparably or better than the baselines evaluated.  ", "questions": "1. Can you elaborate on when and why the combinatorial optimization problem's structure is consistent with the low-rank probability assumption? Some concrete examples can be helpful.\n2. How to balance the exploration and exploitation in the black-box optimization problem? It seems like the proposed algorithm mainly focuses on exploitation (always selecting the top-k candidates). \n", "limitations": "The author did not discuss the potential limitations. An efficient solver for combinatorial black-box optimization can be used for attacks such as threats to cybersecurity. It would be nice to have some discussions on that. ", "rating": 4, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "1. It targets an important and challenging problem interesting to many practical problems.\n2. The idea of using a low-rank tensor can reduce the cost of brute-force search. \n3. Evaluation of different tasks and robustness to the hyperparameters.", "weaknesses": "1. It is a combinatorial search problem from an exponential number of candidates so some structural information of the problem has to be used. However, it is not clear why and when the low-rank probability tensor assumption holds, and if so, for which classes of problems.\n2. The technical contribution is ad-hoc to some extent by combining tensor-train decomposition with REINFORCE trick.\n3. I like the explanation of the intuition starting at line 122 but the approximation in Eq. (5) does not hold rigorously without using the true Monte-Carlo samples as in REINFORCE (the proposed alg instead uses the best k samples)\n4. The complexity analysis assumes a fixed number of queries to the black-box function and does not study the number of iterations required to reach a sufficiently good local minimum. \n5. The proposed method does not show advantages to BS-1 and BS-2 in Table 1 on analytic functions. Figure 3 does not show a clear advantage in the sample efficiency. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you very much for your detailed comments and questions, which allowed us to refine the presentation of the proposed approach. Below we present our respective responses.\n\n1. ... when and why the combinatorial optimization problem's structure is consistent with the low-rank probability assumption?\n\n    - We consider the case of black-box optimization, that is, knowledge of the internal structure of the target function is not assumed; and our optimizer shows good results on several classes of problems at once, unlike baselines. But we agree that the corresponding intuition, based on the knowledge of the internal structure, can be useful. As, for example, the authors of TTOpt and OptimaTT (BS1 and BS2 in our work), are expected their methods to work especially well if the target function has a low-rank tensor structure. However, unlike these methods, we work with a low-rank approximation not of the objective function, but of the probability tensor (P). We expect that the tensor P will be close to the \"delta function\" as a result of optimization, and this function has the exact TT-rank equals 1. At the beginning of the optimization process, the value of the selected rank determines the level of flexibility of the method, however, as shown in the Appendix, varying its value does not lead to a significant deterioration or improvement.\n    \n2. ... ad-hoc ... combining tensor-train decomposition with REINFORCE trick.\n\n    - Thank you for this comment, but we want to emphasize that our algorithm is not based on the REINFORCE trick. After the development and successful testing of the PROTES method, we saw some of its external similarities with the REINFORCE trick and discussed it in the section \"The intuition behind the method\". In particular, the loss function we use (please, see Formula 2) differs from formula 5 used in the REINFORCE trick.\n\n3. ... does not hold rigorously without using the true Monte-Carlo samples as in REINFORCE ...\n\n    - In Formula 5, there is a functional F(f) containing parameters E and T. As we tried to explain in the text after this formula, with an appropriate choice of these parameters in the sum, just the top-k samples will make a noticeable contribution. On the other hand, indeed, with fixed values of the parameters E and T, the optimal number of selected top-k samples can change from iteration to iteration, but we did not see an increase in optimizer performance when using such a modification. We emphasize that our method is not directly based on the REINFORCE trick (the PROTES algorithm is significantly different), at the same time, of course, the study of the connection between the two approaches can give additional intuition.\n    \n4. ... not study the number of iterations required to reach a sufficiently good local minimum.\n\n    - We considered the case of a fixed budget as the most common one in our work (the value of \"sufficient\" accuracy is difficult to formalize and depends on the problem). We note that Table 1 in the Appendix, shows the results of multiple (10) runs of all optimizers for a specific task with a known exact optimum. As follows from the reported results, only our method and BS7 successfully find the exact optimum with 10 restarts, while the average result of our method is significantly better.\n    \n5. ... not show advantages to BS-1 and BS-2 in Table 1 on analytic functions...\n\n    - Yes, indeed, for this class of problems, all three tensor methods (PROTES, TTOpt, and Optima-TT) showed the best results, but for the other three classes of considered problems, the advantages of our approach, as follows from Table 1, are significant. The good result of TTOpt and Optima-TT for analytic functions may be because they have a small tensor rank, in contrast to more complex problems from other classes.\n    \n6. Figure 3 does not show a clear advantage in the sample efficiency.\n\n    - The demonstration on these graphs probably turned out to be not clear enough, since we chose not the best set of benchmarks for demonstration. We conducted a series of additional numerical experiments and described them in detail, as well as we attached a PDF file with plots and tables in our \"global response to all reviewers\". In the PDF file, we present new graphs of optimizer convergence (averaged over 5 restarts) for the considered more complex 1000-dimensional QUBO problem (like P-12 in Table 1). We hope that on these graphs, the advantage of our approach is more noticeable.\n    \n7. How to balance the exploration and exploitation?\n\n    - Thank you for this question! Within the logic of our algorithm, the balance is achieved through the optimal choice of the top-k parameter (that is, the number of good samples selected from the batch). If the value of the parameter is small, then mainly the points close to the current found local optimum are collected (exploitation). On the contrary, if the parameter value is large, then samples from a larger range can be selected, that is, the search area is expanded (exploration). We note that the Appendix contains an analysis of the dependence of the result of our optimizer on the choice of this parameter (as well as other parameters of the method).\n\n8. ... potential limitations ... can be used for attacks such as threats to cybersecurity.\n\n    - We are categorically against the use of our developments for malicious purposes! Unfortunately, like the creators of other powerful tools in the field of optimization, machine learning, quantum physics, etc., we cannot significantly influence this process, since the principle of open source software is also important. To partially remove the risk you indicated, we have already limited in the license the possibility of commercial use of our software product without the explicit approval of the developers, with no restrictions on use for scientific purposes."}, {"review_id": "AqfSQxIDRn", "reviewer": "Reviewer_e361", "summary": "This paper proposes a novel algorithm to estimate global extrema of discrete black-box target functions. The approach proposes learning a probability distribution, represented as a multidimensional tensor and efficiently stored using a TT-decomposition, that is trained to favor multi-indices that achieve high/low values in the target function. The approach, PROTES, draws inspiration from the REINFORCE algorithm, and can be trained using gradient-based optimization techniques. Extensive experiments are conducted to demonstrate the competitive performance and robustness of PROTES.", "questions": "- What are the computational complexities of the baselines?\n- Does this method work especially well when there is a particular structure in the target function? Clearly this method, and any of the baselines would fail if $\\mathcal{Y}$ is random.\n- How is the probability tensor parametrized? Are the values of the tensor learnt directly? Could there be a benefit to alternative parametrizations (i.e. having a multilayer structure)?", "limitations": "Discussed in detail", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The paper is quite self-contained and well written, concepts are introduced gradually and are easy to understand\n- The connection to REINFORCE is informative and helpful\n- Experiments are compelling, detailed ablation analyses are included demonstrating the robustness of PROTE\n", "weaknesses": "Clarity/Notation:\n- \"Gradient-free\" is defined very late in the paper but is used frequently throughout the earlier sections. This is confusing since the optimization of the probability tensor is done with gradient based optimization. I would suggest providing a clear description of this in the introduction.\n- It is stated in the footnote on page 1 that upper case letters are used to denote matrices, but they are used several times in the paper to refer to scalars i.e. $N_i$, $R_i$, etc.\n- In figure 2, the tensor network diagram needs further explanation. It is not clear what is being conveyed from the figure and the current caption.\n\nComparison to Existing Approaches:\n- While the experiments empirically demonstrate the efficacy of the proposed approach, the computation complexities of prior works is not stated anywhere in the paper. Is it the case that the performance benefits are attained by adding significant computational costs? This should be discussed in the comparison.\n- Since the proposed approach is probabilistic, it would be useful to include error bars in the experimental section.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you very much for the high appreciation of our work, very detailed comments and useful suggestions to improve its quality! Below we present our responses to the main points of your review.\n\n1. \"Gradient-free\" is defined very late in the paper...\n\n    - You are absolutely right that this point should be clarified in the text. We added a Footnote in the introduction on Line 34 (\"...and updating its parameters to approximate the optimum...\") explaining that our optimization method is gradient-free, i.e., it does not use the gradient of the objective function (black box) and also does not build an approximation of its gradient. At the same time, to update the internal parameters (elements of the TT-cores of the expectation tensor presented in the TT-format), we use the gradient method. However, in our opinion, this fully preserves the correctness of classifying our method as gradient-free. We also especially note that in some cases, non-gradient approaches are possible for updating the expectation tensor (for example, we can use explicit formulas for the updates of the TT-cores), but we left this minor point outside the scope of this work.\n\n2. ... upper case letters are used to denote matrices... used several times in the paper to refer to scalars...\n\n    - Thank you for noting this! We have added a comment in the Footnote 1 that lowercase and uppercase (where this does not lead to confusion) notation is used for scalars. Renaming uppercase letters like $N_i$ and $R_i$ is difficult since this is a common practice in tensor applications (e.g., $N_i$ is the maximum possible value, and $1 \\leq n_i \\leq N_i$ is a specific value). If you think that the current notation can make the text difficult to read, then we can, for example, use bold capital letters for matrices.\n\n3. ... the tensor network diagram needs further explanation.\n\n    - Unfortunately, due to the limitation on the length of the text, we did not have the opportunity to make the necessary detailed explanations. We will add relevant explanations to the Appendix.\n\n4. ... demonstrate the efficacy of the proposed approach ... What are the computational complexities of the baselines?\n\n    - Thank you for bringing attention to this important point. An estimate for the computational complexity of the proposed method was given in Formula 4 (Section \"Computational complexity of the method\"), and, under certain assumptions, it is close to the complexity estimates for the other two tensor methods chosen as baselines. Also, the Appendix contains the section \"3 Performance comparison of optimization methods\" with Table 2, which presents the measured running time of all methods for the considered benchmarks, and the corresponding discussion. It is important to note that the complexity of the optimization algorithm is largely determined by the complexity of the black box calculation, and if it is significant (and in most practical cases it is), then the determining factor is the number of performed requests to the black box, and not directly the complexity of the optimization algorithm inner operations. Since you have brought this to your attention, we will add a more detailed discussion in the appendix.\n\n5. ...include error bars in the experimental section.\n\n    - Thank you for this comment! We conducted a series of additional numerical experiments and described them in detail, as well as attached a PDF file with plots and tables in our \"global response to all reviewers\". In particular, we presented plots for the mean and uncertainty in a PDF file for a more complex 1000-dimensional QUBO optimization problem that we additionally considered.\n\n6. Does this method work especially well when there is a particular structure in the target function?\n\n    - This is a very important and interesting question, but, unfortunately, it is difficult to give an exact answer to it. As, for example, the authors of tensor optimizers TTOpt and Optima-TT, chosen by us as baselines, are expected their methods to work especially well if the function being optimized has a low-rank tensor structure. However, unlike these methods, we work with a low-rank approximation not of the objective function, but of the probability tensor, and it seems the constraint on the small rank is easier to fulfill in this case (i.e., we expect that the probability tensor will be close to the \"delta function\" as a result of optimization, and this function has the exact TT-rank equals 1).\n\n7. How is the probability tensor parametrized?\n\n    - The probability tensor is represented in the form of a tensor train (TT) decomposition and its parameters are elements of three-dimensional TT-cores (as shown in Figure 2). The use of this parametric representation is the cornerstone of our method because, as discussed in the main text: 1) the number of parameters depends linearly on the dimension; 2) for the TT-decomposition, there is an efficient sampling method with complexity linear in dimension; 3) at the same time, tensor values can also be computed with complexity that is linear in dimension.\n\n8. Could there be a benefit to alternative parametrizations?\n\n    - In view of our comment on the previous point, we consider the TT-decomposition to be an extremely good choice for the parametric representation of the probability tensor. However, you are right, there may be other options as well. In our opinion, it would be interesting to use orthogonalized TT-cores with a special update procedure (this will reduce the number of parameters) or the hierarchical Tucker decomposition, which, like the TT-decomposition, is a low-rank tensor approximation, but has a more complex (tree-like) structure, which can lead to greater expressiveness of the method."}, {"review_id": "OlgXVfqJo9", "reviewer": "Reviewer_U74k", "summary": "The paper proposes PROTES, a new method for multi-dimensinoal discrete black-box optimization based on adaptive probablistic sampling in the tensor train representation. ", "questions": "N/A", "limitations": "The authors discussed the limitations of the paper which include the inability to handle large dimensions (>1000) and explicit constraint consideration.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The technical method is simple and intuitive, with a clear presentation.\n- Diverse benchmark problems and comprehensive comparisons to relevant SOTA baselines.\n- Good empirical performance under the tested experimental setup.", "weaknesses": "- Besides presenting and describing the results, lack of intuition/discussion on why the proposed approaches outperform baselines (both classical and tensor-based). \n- It’s encouraged to have a wider spectrum of dimensionalities of the test problems besides 7 and 50 used.\n- It's not the most efficient method for any given number of requests (evaluation budgets), as seen in Fig 3. The method gains more advantage as requests become closer to 10^4.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you very much for your positive evaluation of our work and your comments! Below we present our responses; we hope that they will eliminate the shortcomings you have noted.\n\n1. ... lack of intuition/discussion on why the proposed approaches outperform baselines ...\n\n    - A potential limitation of the considered tensor-based baselines (TTOpt and Optima-TT) is the requirement for the presence of a low-rank structure for the function being optimized. Since, within the framework of our method, a low-rank tensor approximation (tensor train decomposition) is constructed not for the objective function, but for the probabilistic expectation tensor of the optimum, our method potentially does not have this limitation (at the end of the optimization process, the probabilistic tensor is expected to be close to the \"delta\" function, which has exact TT-rank 1). At the same time, the classical baselines we have considered (PSO,  SPSA, etc.) are predominantly empirical methods, and our advantage can be explained by the successful sampling model used, which effectively combines exploration and exploitation. We also emphasize that the use of the low-rank representation (tensor train) is the cornerstone of our method since the number of parameters depends linearly on the dimension; for the TT-decomposition, there is an efficient sampling method with complexity linear in dimension; at the same time, tensor values can also be computed with complexity that is linear in dimension.\n\n2. It’s encouraged to have a wider spectrum of dimensionalities of the test problems...\n\n    - Thank you for this comment! Indeed, in the original version of the manuscript, the dimensions of the problems did not exceed 100 (we considered 7-dimensional analytic functions; 50-dimensional QUBO problems; 25, 50, and 100-dimensional optimal control problems). We conducted a series of additional numerical experiments and described them in detail, as well as attached a PDF file with plots and tables in our \"global response to all reviewers\". This includes the more complex 1000-dimensional QUBO problem (like P-12 in Table 1), for which we performed computations for various values of the vertex connection probability in the generated random graph. In all calculations, our method outperforms the baselines.\n\n3. The method gains more advantage as requests become closer to $10^4$.\n\n    - In the attached PDF file, we present convergence graphs for the new problem considered (1000-dimensional QUBO problem), from which the advantage of our approach is perhaps more noticeable. We also note that Table 1 in the Appendix, shows the results of multiple (10) runs of all optimizers for a specific problem with a known exact optimum value. As follows from the reported results, only our method and BS7 successfully find the exact optimum with 10 restarts, while the average result of our method is significantly better.\n\n4. ... the inability to handle large dimensions (>1000) and explicit constraint consideration.\n\n    - As follows from the additional numerical experiments, our method remains efficient even when the dimension of the problem is 1000. However, in fact, for the case of higher dimensions, further development of our approach may be needed; in particular, the used TT-decomposition can be replaced by potentially more expressive and compact types of tensor networks, such as hierarchical Tucker decomposition, etc."}, {"review_id": "brq751Meot", "reviewer": "Reviewer_SmA2", "summary": "This paper presents a study on black-box optimization. The authors propose a new method based on sampling from the probability density defined in the tensor train (TT) format (extension effort). The experiments have shown that PROTES outperforms existing efforts.", "questions": "1. Does dimension d have an impact on the experiment? Why was it fixed at 50?", "limitations": "No negative societal impact.", "rating": 6, "confidence": 2, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "(Overall) This paper is well organized. The technical soundness and practical performance are strong. \n\n1. This paper has given a detailed motivation and an analysis of the computational complexity of the method. This makes the work appear comprehensive.\n\n2. The research on gradient-free optimization plays a crucial role in various fields. Despite its significance, it often receives relatively low attention compared to gradient-based optimization methods. This effort tries to add distribution sampling for new tech, \"TT-tensor.\"", "weaknesses": "1. The method section \"Optimization method PROTES\" is hard to follow and what's the relationship between this section and others?\n\n2. The experiment (P1-P20) seems too simple to fully demonstrate the method's effectiveness.\n\n3. This paper also uses the gradient ascent for the TT-tensor, so there lacks more discussion about the model design and gradient-free optimization.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your comments, which allowed us to improve the presentation of the proposed approach! Below we present our responses to your questions on the theoretical part of the work and comments on additional numerical experiments performed.\n\n1. The method section \"Optimization method PROTES\" is hard to follow...\n\n    - We have done our best to make the presentation of the proposed method consistent: in section 2 we have given a simplified description of our approach with a graphical illustration; in section 3 we have briefly described the properties of the tensor train (TT) decomposition that underlies our method; and then in section 4 (\"Optimization method PROTES\"), we have given a formal description of the algorithm, an estimate of its computational complexity, and a possible connection to the REINFORCE trick. Unfortunately, due to page size limitations, we may not have been able to provide a very detailed description of all aspects. If the paper is accepted for the conference, we will of course try to improve the presentation of the method and will certainly take into account your recommendations.\n\n2. The experiment (P1-P20) seems too simple...\n\n    - Thank you for this note! We conducted a series of additional numerical experiments and described them in detail, as well as attached a PDF file with plots and tables in our \"global response to all reviewers\". We consider four challenging problems of maximization of the cumulative reward of the episode for RL agents like InvertedPendulum, Swimmer, Lunar Lander, and Half Cheetah from Mujoco / OpenAI-GYM collection, and we consider the 1000-dimensional case for our benchmark P-12 (Minimum Vertex Cover problem). \n\n3. This paper also uses the gradient ascent for the TT-tensor ... gradient-free optimization.\n\n    - Thank you to note this. We added a Footnote in the introduction on Line 34 (\"...and updating its parameters to approximate the optimum...\") explaining that our optimization method is gradient-free, i.e., it does not use the gradient of the objective function (black box) and also do not build an approximation of its gradient. At the same time, to update the internal parameters (elements of the TT-cores of the expectation tensor to the optimum presented in the TT-format), we use the gradient method. However, in our opinion, this fully preserves the correctness of classifying our method as gradient-free. We also especially note that in some cases, non-gradient approaches are possible for updating the expectation tensor (for example, we can use explicit formulas for the updates of the TT-cores).\n\n4. Does dimension d have an impact on the experiment?\n\n    - We considered three classes of problems (the corresponding results were presented in Table 1): analytic functions (P-01 - P-10), QUBO problems (P-11 - P-14), and optimal control problems (P-15 - P-20; with and without constraints). For analytical functions, we have chosen dimension 7, since among them there is a Piston function (P-06) that is interesting from a practical point of view and has exactly this dimension. For optimal control problems, we considered the dimensions 25 (P-15 and P-18), 50 (P-16 and P-19), and 100 (P-17 and P-20); since our method significantly outperforms the baselines already at a dimension value of 100, we did not consider large values of the dimension. However, for QUBO problems, the chosen dimension value of 50 may indeed seem too low. We conducted a series of additional experiments for a dimension value of 1000 for QUBO problems. The results are presented in the attached PDF file. Due to the page limit, we report the result for only one of the problems, i.e., for P-12, but we considered various formulations of this problem corresponding to various options for initializing a random graph."}, {"review_id": "XerQoWd06j", "reviewer": "Reviewer_7wCj", "summary": "The authors present the algorithm PROTES for solving discrete optimization problems. The method uses probabilistic sampling to obtain a probability function of minimal (or maximal) values of a target function. This probability function takes the same shape as the target tensor. Alternatively, this could be described as a likelihood score for the optimal values of the discrete variables. To make the algorithm more efficient, the authors use the tensor train (TT)-representation for tensors. This is a low-rank approximation of a tensor, which makes the iterative computations faster. This TT-representation is used for the probability function.\n\nThe authors test their method on 20 different problem instances and compare their method to seven other methods. The results show that PROTES works well, and is stable for the instances it was tested on.", "questions": "- What is the reason for using gradient ascent for updating the the probability distribution? Cannot alternatively a cheaper updating rule be used? Think for instance of increasing the likelihood of the selected candidates by some value in each iteration (similar to weighing).\n\n- In line 62 of the paper, it is mentioned that the probability distribution is expected to present a Kronecker delta-function. It sounds more like an indicator function with a peak in the value of the minimum target function. Moreover, why is this function expected?\n\n- Please address why only one instance was tested per problem (size), and why this is enough evidence of the efficacy of the method for the type of problem.\n\n- In line 182, it is mentioned that related methods in the TT-format should be used for approximation. Is the approximation of the optimum meant here? It is a bit confusing since in the sentence before the (TT-)approximation of the original tensor is mentioned.\n\n- What are the specifications (CPU etc) of the ‘regular laptop’?", "limitations": "Not applicabel.", "rating": 5, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "The paper has a clear structure and is well written.", "weaknesses": "The method was not well motivated and the update rule was unclear to the reader. This makes is it very hard to assess whether there is sufficient contribution to the research area.\n\nThe comment on the application of the method on constrained optimization is too short-sighted. In other words, such an important topic deserves more explanation. The current version does not clearly explain how general constrained optimization problems can be handled, particularly the feasibility concerns.\n\nThe instance sizes of the QUBO problems are very small compared to what is now possible with the constrained version of these problems.\n\nAlthough the method was tested on 20 different problems, compared to seven benchmarks, only one instance per problem (size) was tested. For instance, for P-11 to P-13, it looks like only one instance was selected, which was not specified whether random problem parameters or something else was used as input. Therefore, the reader cannot tell if the proposed method works well in general for the tested problems or only on the tested instance.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for the detailed analysis of our work and extensive comments and remarks! We tried to take them into account and make appropriate improvements. Below we present (within the specified limit on the length of the sent message) our responses to the main points of your review.\n\n1. The method was not well motivated and the update rule was unclear to the reader.\n\n    - Indeed, as in the case of many other popular optimization methods, our approach is predominantly heuristic, so we confirm its effectiveness by comparing it with alternative known methods on a wide class of problems. Due to text size restrictions, we, unfortunately, could not describe in great detail all aspects. The scheme for updates we use (line 9 of Algorithm 1) is based on the automatic differentiation of our loss function from Formula 2, the rationale for which and connection with the REINFORCE trick is given in the section \"The intuition behind the method\". We note that due to the use of the tensor train (TT) decomposition, we can simultaneously efficiently (with the complexity linear in dimension) compute the likelihood for the automatic differentiation and sample from the probability tensor.\n\n2. The comment on the application of the method on constrained optimization is too short-sighted.\n\n    - Due to the text size limitations, we could not describe this point in more detail, however, we will add additional comments to the Appendix. We note that even without explicitly taking into account the constraints (according to the scheme we described in the paper), our method is effective in comparison with baselines.\n\n3. The instance sizes of the QUBO problems are very small... \n\n    - Thank you for this comment! We conducted a series of additional numerical experiments and described them in detail, as well as attached a PDF file with plots and tables in our \"global response to all reviewers\". In particular, we consider a 1000-dimensional QUBO problem (like P-12 in Table 1) in five different settings. As follows from the results in the PDF file, our optimizer PROTES demonstrates the best result. At the same time, we understand that there are specialized solvers for this class of problems that may be potentially more effective (at least in terms of speed), but we compared our general-purpose optimization method with similar non-specialized approaches.\n\n4. ... only one instance per problem (size) was tested ...\n\n    - Thank you! We tried to take this remark into account by conducting a series of additional experiments (the results and a detailed discussion are presented in our \"global response to all reviewers\"). In particular, for the 1000-dimensional QUBO problem (like P-12 in Table 1) we consider 5 different options for initializing a random graph and also averaged the results for 5 independent runs of all optimizers.\n\n5. What is the reason for using gradient ascent for updating the probability distribution?\n\n    - Since we are using a parametric representation (TT-decomposition) of the probability tensor, we can change the parameters, but not the value at one point or set of points (in this case, the change in the tensor turns out to be global, that is, a certain \"learning\" process takes place). We considered several variants close to that proposed by you (e.g., addition in the TT-format with subsequent rounding of the result in the TT-format, since the rank increases with addition), but we did not find significant improvement in the optimization results. Thus, we use automatic differentiation to update the TT-tensor parameters, since this is the most simple and common option, however, you are absolutely right that finding more efficient alternatives is a good direction for further research.\n\n6. ... the probability distribution is expected to present a Kronecker delta-function ... why is this function expected?\n\n    - Thank you for this clarification about the indicator function. We expect this behavior, since when sampling, values closer to the optimum will more often be selected (after the first successful choice), under certain conditions for smoothness. We have tested this intuition with a number of examples and it was carried out, but it is obvious that certain counterexamples can be also found (strongly isolated optimum, multiple global optima, etc.). Also, due to budget constraints, this behavior is not fully realized. We have added relevant clarifications to the text of the paper.\n\n7. In line 182, it is mentioned that related methods in the TT-format should be used for approximation...\n\n    - Thank you for this note! We have corrected the wording in the text of the paper. The Optima-TT essentially is a method for optimizing a surrogate model presented in the TT-format, that is, first, the TT-approximation (surrogate model) of the objective function is built using one of the existing approximation methods (TT-SVD, TT-ALS, TT-cross), and then the optimum of this model is searched with the Optima-TT. As we indicated in Footnote 5, when using this method as a baseline, we additionally applied the TT-cross method to build an approximation (surrogate model).\n\n8. What are the specifications (CPU etc)?\n\n    - Numerical experiments presented in the paper were performed with the MacBook Pro (Intel Core i7, RAM 16 Gb). The wording \"regular laptop\" in the text meant that our method is not resource-intensive and does not require special computing equipment (however, it may certainly be needed in case of too complex objective functions for their efficient calculation).\n"}]

Author rebuttal:
Thank you very much for your comments, which allowed us to correct the shortcomings and refine the presentation of the proposed approach. We have improved the text of the manuscript, however, according to the rules of the conference, we are now unable to attach the updated version.

1. We practically did not meet major remarks from the Reviewers on the algorithmic part of the proposed new method PROTES for the multidimensional gradient-free optimization, but taking into account the comments, we will make three general replies on the theoretical part:

    - 1.1. In our approach, the probability tensor for the optimum expectation (P) is represented in the form of a low-rank (low-parametric) tensor train (TT) decomposition, and we iteratively optimize its parameters (i.e., the elements of three-dimensional TT-cores) to sample candidates for the optimum from this tensor. The use of this parametric representation is the cornerstone of our method since the number of parameters depends linearly on the dimension; there is an efficient sampling method in the TT-format with linear complexity; values of the TT-tensor can also be computed with linear complexity.

    - 1.2. Each step of updating the probability tensor turns out to be non-local: the probability increases not only for the selected best candidates ("top-k") but also the probabilities for all other multi-indices are redistributed within the framework of the parametric structure used. In this case, with a wide range of values of the "top-k" parameter, the balance between exploration and exploitation is maintained.

    - 1.3. We demonstrate the formal connection of our approach to the REINFORCE trick in the text, however, PROTES is not based on it. The algorithm and the loss function used are significantly different.

2. We tried to take into account the comments of the Reviewers on the experimental part of our work (i.e., insufficiently complex benchmarks; low values of dimensions for QUBO problems) by conducting a series of additional experiments. We report the results of our experiments in the attached PDF file (Table 1, Table 2, and Figure 1) and discuss them below:

    - 2.1. As in the work "TTOpt: A maximum volume..." (NeurIPS-2022), we consider four challenging problems of maximization of the cumulative reward of the episode for RL agents: InvertedPendulum (P-21), Swimmer (P-22), Lunar Lander (P-23), and Half Cheetah (P-24) from Mujoco / OpenAI-GYM collection. The policy is represented by a neural network with three layers and tanh activations. Parameters of the neural network are discretized on a grid with limits from -1 to +1 and 256 nodes (i.e., mode size N). The total number of neural network parameters for the problems under consideration turns out to be from 50 to 68 (i.e., dimension d). Thus, as in the work "TTOpt", we obtain a discrete on-policy learning problem (search for the values of the neural network parameters that lead to the maximum reward). To speed up the calculations, we took only two baselines (BS1, and BS7) for comparison, which showed the best results in the main experiments presented in the paper. We used the same optimizer settings as in the main calculations described in the text of the paper, however, in this case, we chose the budget value equal to M = 100 K, as in the work "TTOpt". For PROTES, BS1 and BS7, we repeated the calculations three times and reported the average result for found maximum and computation time in Table 1. As we can see, PROTES shows the best result in terms of accuracy in 3 cases out of 4, and once it is in second place, slightly inferior to the BS7. We propose to add Table 1 with a detailed text description to the new section-1 of the Appendix.

    - 2.2. We consider the 1000-dimensional case for our benchmark P-12 (Minimum Vertex Cover problem). Additionally, we vary the value of the vertex connection probability $p$ in the generated random graph (i.e., the parameter in the corresponding function from the well-known networkx package), getting 5 different optimization problems with p=0.1, p=0.2, p=0.3, p=0.4 and p=0.5 (please note that in the original text of the paper, we used the value p=0.5). We performed the calculations with PROTES and baselines BS1 - BS7 with the same settings as in the main computations described in the text of the paper. In Table 2 we report the results averaged over 5 runs of each of the optimizers for each value of p (we note that BS2 failed for this problem because the selected dimension value is too large to construct the approximation of the tensor). As we can see, in all cases our optimizer PROTES demonstrates the best result. We propose to add Table 2 with a detailed text description to the new section-2 of the Appendix.

    - 2.3. Using the numerical results from the previous paragraph (2.2), we plotted the convergence curves for p=0.1, p=0.3, p=0.5, and all optimizers (except the BS2, since it failed), showing the average result and the spread of values over 5 restarts. As we can see, our optimizer PROTES has a small result variance and consistently shows the best result when the number of iterations is more than 1000. We propose to add new Figure 1 to the new section-2 of the Appendix or replace Figure 3 in the main text with it.

We once again thank all the Reviewers for the positive assessment of our work, as well as for useful comments and advice that allowed us to improve its presentation! We will be happy to answer any additional questions you may have.

Meta-review: This paper proposes a new method for black-box optimization over a sizable number of combinatorial (e.g., binary) parameters. Different to earlier work, the authors use an efficient tensor train (TT) representation, not for the unknown function, but to parameterize a distribution over it (i.e., a surrogate model).

This paper tackles an important and understudied problem in ML, where many BO methods use Gaussian process surrogate models, which do not work well for combinatorial spaces. The reviewers are also convinced about the empirical evaluations. Reviewers criticize that it is not clear which type of combinatorial problems are best approximated by the surrogate model here, but the authors convincingly argue that this is not an easy question to answer, but that the assumptions are certainly weaker than for prior work, where the function itself has to follow a TT structure (here, it is only the distribution over functions). Another criticism is that the problems in the experiments are fairly small in terms of number of parameters, but the authors post favourable results for another problem with 1000 variables. They also address criticism that the original submitted results are for single problem instances. I'd like to highlight the high quality and measuredness of the feedback, which seems to have satisfied the reviewers, leading to score increases (and myself as AC). If this paper gets accepted, the authors should carefully take the recommendations of the reviewers for more clarity into account. For example, there was some confusion about REINFORCE, whereas the authors seem to use automatic differentiation, so the relation to REINFORCE does not seem close.
Decision: Accept (poster)