Black-Box Approximation and Optimization with Hierarchical Tucker Decomposition

## Abstract

Abstract We develop a new method HTBB for the multidimensional black-box approximation and gradient-free optimization, which is based on the low-rank hierarchical Tucker decomposition with the use of the MaxVol indices selection procedure. Numerical experiments for 14 complex model problems demonstrate the robustness of the proposed method for dimensions up to 1000, while it shows significantly more accurate results than classical gradient-free optimization methods, as well as approximation and optimization methods based on the popular tensor train decomposition, which represents a simpler case of a tensor network.

## 1 Introduction

Many physical and engineering models can be represented as a real function (output), which depends on a multidimensional argument (input) and looks like

$$ $y=\textsf{f}(\bm{x})\in\mathbb{R},\quad\bm{x}=[x_{1},\,x_{2},\,\ldots,\,x_{d}]^{T}\in\Omega\subset\mathbb{R}^{d}.$ (1) $$

Such functions often have the form of a black-box (BB), i. e., the internal structure and smoothness properties of f remain unknown.
Its discretization on some multi-dimensional grid results in a multidimensional array (tensor(^3^33 By tensors we mean multidimensional arrays with a number of dimensions $d$ ($d\geq 1$). A two-way tensor ($d=2$) is a matrix, and when $d=1$ it is a vector. For scalars we use normal font, we denote vectors with bold letters and we use upper case calligraphic letters ($\mathcal{A},\mathcal{B},\mathcal{C},\ldots$) for tensors with $d>2$. The $(n_{1},n_{2},\ldots,n_{d})$th entry of a $d$-way tensor $\mathcal{Y}\in\mathbb{R}^{N_{1}\times N_{2}\times\ldots\times N_{d}}$ is denoted by $\mathcal{Y}[n_{1},n_{2},\ldots,n_{d}]$, where $n_{k}=1,2,\ldots,N_{k}$ ($k=1,2,\ldots,d$) and $N_{k}$ is a size of the $k$-th mode. Mode-$k$ slice of such tensor is denoted by $\mathcal{Y}[n_{1},\ldots,n_{k-1},:,n_{k+1},\ldots,n_{d}]\in\mathbb{R}^{N_{k}}$.)) $\mathcal{Y}\in\mathbb{R}^{N_{1}\times N_{2}\times\ldots\times N_{d}}$ that collects all possible discrete values of the function ([1](#S1.E1)) inside the domain $\Omega$, i.e.,

$$ $\mathcal{Y}[n_{1},\,n_{2},\,\ldots,\,n_{d}]=\textsf{f}\bigl{(}x_{1}^{(n_{1})}\!,\,x_{2}^{(n_{2})}\!,\,\ldots,\,x_{d}^{(n_{d})}\bigr{)}.$ (2) $$

Storing such a tensor often requires too much computational effort, and for large values of the dimension $d$, this is completely impossible due to the so-called curse of dimensionality (the memory for storing data and the complexity of working with it grows exponentially in $d$).
To overcome it, various compression formats for multidimensional tensors are proposed: Canonical Polyadic decomposition aka CANDECOMP/PARAFAC (CPD) (Harshman et al., ), Tucker decomposition (Tucker, ), Tensor Train (TT) decomposition (Oseledets, ), Hierarchical Tucker (HT) decomposition (Hackbusch & Kühn, ; Ballani et al., ), and their various modifications.
These approaches make it possible to approximately represent the tensor in a compact low-rank (i.e., low-parameter) format and then operate with the compressed tensor.

The TT-decomposition is one of the most common compression formats (Cichocki et al., , ).
There is an algebra for tensors in
the TT-format (i.e., TT-tensors):
we can directly add and multiply TT tensors,
truncate TT tensors (reduce the so-called TT-rank, i. e, the number of storage parameters), integrate and contract TT tensors.
It is important that effective algorithms have been developed (Kapushev et al., ; Ahmadi-Asl et al., ; Chertkov et al., ) for approximating BB like ([1](#S1.E1)) and ([2](#S1.E2)) in the TT-format, that is, for constructing an approximation (surrogate model) using only a small number of explicitly computed BB values.
In recent years, new efficient algorithms have also been proposed (Sozykin et al., ; Nikitin et al., ; Chertkov et al., ) for the second important problem associated with gradient-free optimization of such BB, that is, finding an approximate minimum or maximum value based only on queries to the BB.
Such TT-based methods of surrogate modeling (in particular, the TT-cross algorithm (Oseledets & Tyrtyshnikov, )) and gradient-free optimization (in particular, the TTOpt algorithm (Sozykin et al., )) have shown their effectiveness for various multidimensional problems, including compression and acceleration of neural networks, data processing, modeling of physical systems, etc.

However, the TT-decomposition is one of the simplest special cases of a tensor network: it is a linear network or a degenerate tree, and it has a number of limitations related to weak expressiveness and instability for the case of significantly large dimensions.
The HT-format is potentially more expressive and robust (Buczyńska et al., ); thus, it makes it possible to approximate more complex functions with fewer parameters.
Taking into account TT-Cross and TTOpt algorithms which use the well-known MaxVol approach (Goreinov et al., ; Mikhalev & Oseledets, ), in this work we develop new methods of surrogate modeling and gradient-free optimization based on the HT-format,
and our main contributions are the following:

- •
we develop a new black-box approximation method HT-cross based on the HT-decomposition and the rectangular MaxVol index selection procedure;
- •
we develop a new gradient-free optimization method HTOpt based on the HT-decomposition and the rectangular MaxVol index selection procedure;
- •
we implement the proposed HT-cross and HTOpt algorithm as a unified method HTBB for surrogate modeling and optimization of multidimensional functions given in the form of a black-box and share it as a publicly available python package;(^4^44 The program code with the proposed approach and numerical examples, given in this work, is publicly available in the repository [https://github.com/G-Ryzhakov/htbb](https://github.com/G-Ryzhakov/htbb).)
- •
we apply our approach HTBB to $14$ different complex model functions with input dimensions up to $1000$ and demonstrate its significant advantage in the accuracy and robustness for the same budget in comparison with the TT-cross method for approximation problems and with the TTOpt, and well-known classical gradient-free SPSA and PSO methods for optimization problems.

## 2 Hierarchical Tucker Decomposition

By Hierarchical Tucker (HT), we mean a tensor tree that is not necessarily balanced (Ballani et al., ).
Let us describe this concept in detail in the context of our work.
HT is such a low-parameter decomposition of a $d$-way tensor, which is a hierarchical contraction of $3$-way tensors and $2$-way tensors, ordered in the form of a binary tree.

Consider a binary tree — a graph without cycles,
where every *node* (except the root one) has a *parent* and at most two *children*.
In what follows, we consider trees where each node has either 2 or 0 children.
We call a node without children a *leaf*.
We denote the depth of the tree by $L$, and the number of nodes at level $l$ (starting from the root node) by $\lambda_{l}$; note that for a balanced tree, $\lambda_{l}=2^{l-1}$ is satisfied.
With each node, we associate a *core* tensor, i.e., a 2-way tensor with the leaves,
and 3-way tensors with all others
(for the root core we add a dummy dimension of the length $1$).

The number of leaves $d$ determines the dimensionality of the considered tensor $\mathcal{Y}$,
which is represented in the described tree structure, i.e., $\mathcal{Y}\in\mathbb{R}^{N_{1}\times N_{2}\times\cdots\times N_{d}}$, where $N_{j}$ is the size of $j$th mode.
The dimensions of the cores are as follows.
Leaves dimensions correspond to the dimensionality of the tensor $\mathcal{Y}$:
each core ${\mathcal{G}}^{(L)}_{j}$ that is associated
with a leaf node with number $j$ satisfies ${\mathcal{G}}^{(L)}_{j}\in\mathbb{R}^{r^{(L)}_{j}\times N_{j}}$.
The dimensions of the non-leaves cores match the dimensions of their children:
if core ${\mathcal{G}}^{(l)}_{j}\in\mathbb{R}^{r^{(l+1)}_{j_{1}}\times r^{(l)}_{j}\times r^{(l+1)}_{j_{2}}}$
for $1\leq l<L$, then its child ${\mathcal{G}}^{(l+1)}_{j_{1}}$ and ${\mathcal{G}}^{(l+1)}_{j_{2}}$ have such dimensions that
${\mathcal{G}}^{(l+1)}_{j_{1}}\in\mathbb{R}^{r_{1}\times r^{(l+1)}_{j_{1}}\times r_{2}}$
and
${\mathcal{G}}^{(l+1)}_{j_{2}}\in\mathbb{R}^{r_{3}\times r^{(l+1)}_{j_{2}}\times r_{4}}$.
The numbers $r^{(i)}_{j}$
are called *ranks* of the HT decomposition.
For the core of the root node it is hold ${\mathcal{G}}^{(1)}_{1}\in\mathbb{R}^{r^{(2)}_{j_{1}}\times 1\times r^{(2)}_{j_{2}}}$.
In the case of notations related to tree nodes,
the index at the top in parentheses denotes the level of the tree $l$, it varies from $1$ to $L$ ($L=\ln d$ for the balanced tree),
and with the index at the bottom we denote the numbering within this level of the tree, this numbering is not fixed and is arbitrary.

To calculate the value of the tensor $\mathcal{Y}$ represented in the HT-format at a given index $I$, we perform the following iterative procedure.
We associate a vector ${b}^{(l)}_{j}$ with each node, which is recursively defined as

$$ ${b}^{(l)}_{j}=\sum_{i=1}^{r_{1}}\sum_{k=1}^{r_{2}}{\mathcal{G}}^{(l)}_{j}[i,\,:,\,k]\cdot{b}^{(l+1)}_{j_{1}}[i]\cdot{b}^{(l+1)}_{j_{2}}[k],$ $$

where the vectors
${b}^{(l+1)}_{j_{1}}$ and ${b}^{(l+1)}_{j_{2}}$ are vectors associated with children of the current node;
$r_{1}={r}^{(l+1)}_{j_{1}}$ and
$r_{2}={r}^{(l+1)}_{j_{2}}$ are the corresponding ranks.
For a leaf node,
its corresponding vector ${b}^{(L)}_{j}$ depends on the given index $I$ as ${b}^{(L)}_{j}={\mathcal{G}}^{(L)}_{j}\bigl{[}:,\,I[j]\bigr{]}.$ Finally,
the resulting tensor value at index $I$ is equal
to the value of the single element of the vector ${b}^{(1)}_{1}$
associated with the root node, i.e.,
$\mathcal{Y}[I]={b}^{(1)}_{1}[1].$ Note that this procedure is easily parallelized naturally since vectors $b$ of the same level in different parts of the tree are calculated independently.
Moreover, the HT-format is more expressive and robust (Buczyńska et al., ) than simpler forms of tensor networks (for example, the well-known TT-decomposition), which makes it potentially possible to build the low-rank approximation for complex functions.

## 3 Proposed Approach

Due to the potential strengths of the low-rank HT-decomposition for high-dimensional applications described in the previous section, it seems important to develop new approximation and optimization methods based on it.
We are inspired by a simpler, but carefully designed TT-format and implement analogues of the known methods TT-cross and TTOpt on its basis for the HT-decomposition.
The TT-cross algorithm (Oseledets & Tyrtyshnikov, ) adaptively calls the BB and iteratively builds the TT-surrogate until a given accuracy is reached or the BB access budget is exhausted.
During this construction, the so-called Maximum Volume submatrix search (MaxVol) procedure (Goreinov et al., ) is used to find a close to the dominant matrix of the tensor unfolding.
Thus, this matrix contains values close to the maximum modulus values of the tensor.
This effect can be used to find the quasi-maximal element in the tensor, and the corresponding algorithm is called TTOpt (Sozykin et al., ).
Further in this section, we successively describe our algorithms for approximation (HT-cross) and optimization (HTOpt) in the HT-format.

### 3.1 Upper and Down Indices

The key concept that is used for both the approximation and optimization algorithm is to associate index sets with each link between nodes.
Each link between node ${D}^{(l-1)}_{m}$ and its child ${D}^{(l)}_{j}$ have *down* $i^{\text{\scriptsize\tt down}}_{l,j}$ and *upper* $i^{\text{\scriptsize\tt up}}_{l,j}$ indices and corresponding values $v^{\text{\scriptsize\tt down}}_{l,j}$ and $v^{\text{\scriptsize\tt up}}_{l,j}$ of this indices.
Since each link is unambiguously defined by the child node ${D}^{(l)}_{j}$ it is part of, the index notations are similar to this children node notation and sometimes we refer to these indices as being associated with the child node rather than a relation.

Down $i^{\text{\scriptsize\tt down}}_{l,j}$ and upper $i^{\text{\scriptsize\tt up}}_{l,j}$ indices depend only on their position and are fixed during initial tree construction according to the following recursive rule.
Each leaf node ${D}^{(L)}_{j}$ has upper index $i^{\text{\scriptsize\tt up}}_{L,j}=\{j\}$ containing one element equal to the element number of the tensor index element that is associated with this leaf node.
Each non-leaf node ${D}^{(l)}_{j}$ except the root one has an upper index consisting of the union of the elements of the upper indices of all its children: $i^{\text{\scriptsize\tt up}}_{l,j}=i^{\text{\scriptsize\tt up}}_{l+1,j_{1}}\cup i^{\text{\scriptsize\tt up}}_{l+1,j_{2}}$.
For all the cases, down indices are equal to the set difference between all tensor indices and upper indices: $i^{\text{\scriptsize\tt down}}_{l,j}=\{1,\,2,\,\ldots,\,d\}\setminus i^{\text{\scriptsize\tt up}}_{l,j}$.
Since the root node is not a child, we do not associate indexes with it.
The down indices of the left child of the root node and their values are equal to the upper indices of the right child and their values, respectively, and vice versa.
Please see [Figure 1](#S3.F1) for relevant illustration.

Figure: Figure 1: Examples of upper and down indices and their values for $\mathcal{Y}\in\mathbb{R}^{N_{1}\times N_{2}\times N_{3}\times N_{4}\times N_{5}\times N_{6}\times N_{7}\times N_{8}}$ with $N_{1}=N_{2}=N_{5}=N_{6}=N_{7}=2$, $N_{3}=N_{4}=3$, and $N_{8}=10$.
Refer to caption: /html/2402.02890/assets/x1.png

Figure: (a) Update upper indices values.
Refer to caption: /html/2402.02890/assets/x2.png

Figure: Figure 3: Examples of a path for the traversal procedure. The task is 5-dimensional, so indices 5, 6, and 7 (green boxes) as well as their parents (purple circles) are never visited. Ranks of all links, except for those leading to inactive indices, are equal to 3.
Refer to caption: /html/2402.02890/assets/diagrams/tree_big_v2.png

The values of the upper $v^{\text{\scriptsize\tt up}}_{l,j}$ and down $v^{\text{\scriptsize\tt down}}_{l,j}$ indices change dynamically and the manner and sequence of their change is the subject of this study.
These values $v^{\text{\scriptsize\tt up}}_{l,j}$ and $v^{\text{\scriptsize\tt down}}_{l,j}$ represent sets of size equal to the rank, associated with the corresponding node:
$\left|v^{\text{\scriptsize\tt up}}_{l,j}\right|=\left|v^{\text{\scriptsize\tt down}}_{l,j}\right|={r}^{(l)}_{j}$.
Each element of this set is a vector with values of indices stored in the corresponding ($i^{\text{\scriptsize\tt up}}_{l,j}$ or $i^{\text{\scriptsize\tt down}}_{l,j}$) index set.
The main goal of the iterative search for index values (the detailed implementation of which will be described below) is to find the submatrix of maximum volume at the intersection of the given indices.
Finding a submatrix of maximal volume serves two purposes: first, we can more accurately reconstruct the original matrix using it, and second, we expect that this matrix has elements close to maximal in modulo.
Let us elaborate on the construction of this matrix.

Let $\mathcal{Y}^{U}(I)$ for the given index be the *unfolding* matrix of the $d$-way tensor $\mathcal{Y}$ in the given index $I=({k}_{1},\,{k}_{2},\,\ldots,\,{k}_{n})$, $1\leq n\leq d$, if for all its elements holds

$$ $\mathcal{Y}^{U}(I)[\overline{i_{k_{1}}\cdots i_{k_{n}}},\,\overline{i_{p_{1}}\cdots i_{p_{d-n}}}]=\mathcal{Y}[i_{1},\,i_{2},\,\ldots,\,i_{d}],\\ \{{p}_{1},\,{p}_{2},\,\ldots,\,{p}_{d-n}\}=\{1,\,2,\,\ldots,\,d\}\setminus I.$ $$

By a line on a group of indices, we mean a multi-index composed of the given indices, i. e. the position of the corresponding sequence of indices in the list of all possible values.
We do not fix a particular sorting type of this sequence (lexicographic order can be taken) since the rearrangement
does not affect the rank of the matrix or the property of its submatrix of maximal volume.

For a non-leaf node ${D}^{(l)}_{j}$, the node up indices $i^{\text{\scriptsize\tt down}}_{l,j}$ and up and down indices values $v^{\text{\scriptsize\tt up}}_{l,j}$ and  $v^{\text{\scriptsize\tt down}}_{l,j}$ we can construct the unfolding $\mathcal{Y}^{U}_{l,\,j}$ as
$Y_{l,\,j}=\mathcal{Y}^{U}(i^{\text{\scriptsize\tt up}}_{l,j}).$ If we consider submatrix $Y_{l,\,j}[v^{\text{\scriptsize\tt up}}_{l,j},\,v^{\text{\scriptsize\tt down}}_{l,j}]\in\mathbb{R}^{r\times r}$
of this matrix based on the values $v^{\text{\scriptsize\tt up}}_{l,j}$ and  $v^{\text{\scriptsize\tt down}}_{l,j}$,
where $r={r}^{(l)}_{j}$ is the corresponding rank, when our goal is to make a volume of the matrix $Y_{l,\,j}[v^{\text{\scriptsize\tt up}}_{l,j},\,v^{\text{\scriptsize\tt down}}_{l,j}]$
as large as possible by choosing indices values $v^{\text{\scriptsize\tt up}}_{l,j}$ and $v^{\text{\scriptsize\tt down}}_{l,j}$.
Recall that the *volume* of any (tall) matrix $A$ is defined as

$$ $\operatorname{vol}A=\sqrt{\det A^{T}A},\quad A\in\mathbb{R}^{n\times m},\quad n\geq m.$ $$

### 3.2 Index Values Update Algorithm

Figure: Algorithm 1 Indices values update algorithm.

While our method is running, we update all indices (both up and down), using the same procedure, as presented in Algorithm [1](#alg1).
Here the function $zeros$ reserves the specified number of elements for a vector, matrix, etc.
The $\operatorname{QRP}$ function returns a QR decomposition with permutations (i. e., the elements on the diagonal of $R$ do not decrease; we use the implementation from the Python package scipy).
The operation of the algorithm can be briefly described as follows.
We construct a tall matrix, whose rows correspond to the tensor product of index values, which are conditionally called “incoming” and columns to “outgoing” ones.
Then, using the $\operatorname{MaxVol}$ procedure, we select rows from this matrix so that the submatrix corresponding to them is of quasi-maximal volume, and the index values corresponding to these rows are returned.

#### Algorithm Input Indices.

The “incoming” $i_{1}$, $i_{2}$, and “outgoing” $i$ input indices in our algorithm depend on the indices that are updated at each step (see [Figure 2](#S3.F2)).
Namely, if we update the upper indices values $v^{\text{\scriptsize\tt up}}_{l,j}$ for some node ${D}^{(l)}_{j}$, then “incoming” indices are the upper indices of children ${D}^{(l+1)}_{{j}_{1}}$ and ${D}^{(l+1)}_{{j}_{2}}$ of this node: $i_{1}\leftarrow i^{\text{\scriptsize\tt up}}_{l+1,j_{1}}$, $i_{2}\leftarrow i^{\text{\scriptsize\tt up}}_{l+1,j_{2}}$.
The “output” indices are down indices of the node ${D}^{(l)}_{j}$:
$i\leftarrow i^{\text{\scriptsize\tt down}}_{l,j}$ (see [Figure 2(a)](#S3.F2.sf1)).
If, in turn, we update the values of down indices for the link that connects parent ${D}^{(l)}_{j}$ and child ${D}^{(l+1)}_{{j}_{1}}$, then for “incoming” indices we have: $i_{1}\leftarrow i^{\text{\scriptsize\tt down}}_{l,j}$, $i_{2}\leftarrow i^{\text{\scriptsize\tt up}}_{l+1,j_{2}}$, where $j_{2}$ is the number of another child ${D}^{(l+1)}_{{j}_{2}}$ of the node ${D}^{(l)}_{j}$ (which differs from the original child ${D}^{(l+1)}_{{j}_{1}}$).
The upper indices of the given link are the “output” ones: $i\leftarrow i^{\text{\scriptsize\tt up}}_{l+1,j_{1}}$ (see [Figure 2(b)](#S3.F2.sf2)).

#### Transformation of the Tensor Values.

The point-wise transformation is needed when we search for the minimum.
In this case, we can transform tensor values by any monotonic decreasing function $T$.
In our experiments, we use the following adaptive (i. e. its parameters dependent on the given data batch) transformation $T$

$$ $\displaystyle T(x)=\exp\bigl{(}-(x-x_{0})/\sigma\bigr{)},$ $\displaystyle x_{0}=\operatorname{mean}(x),\quad\sigma=\operatorname{std}(x),$ $$

where $\operatorname{mean}(x)$ and $\operatorname{std}(x)$ are the sample mean and sample variance of the set of numbers respectively.
When searching for the maximum value, we do a similar transformation (note that transformation can be avoided in this case; however, we apply it for greater stability of the method)

$$ $T(x)=\exp\bigl{(}(x-x_{0})/\sigma\bigr{)}.$ $$

#### MaxVol Procedure.

The $\operatorname{MaxVol}$ procedure in our algorithm is the so-called *rectangular* maximum volume search method (Mikhalev & Oseledets, ).
Note that it can return not only square matrices but also rectangular matrices, making a decision on the number of returned rows based on a heuristic procedure based on the possible increase in volume when adding a candidate row and the given tuning parameters.
In our numerical experiments, we allowed to expand output index set by at most $\Delta r=1$ element, so the ranks grew by at most 1 per pass.

### 3.3 Traversal Procedure

When updating indices, we walk sequentially to the neighboring (linked) node, going back only if we reach a leaf node.
As each visited node, we increment its visit counter by one, whereas at the beginning all counters were reset to zero.
When we pass through an edge, we update only one set of index values at a time: if we go from parent to child, we update down indices values; if we go from child to parent, we update upper indices values.

To decide which of the two nodes to take the next step to (in case there are two options), we count the average number of visits in each part of the tree that separates each of the two paths.
Namely, we cut the edge that was traveled last, and we cut the edges connecting the current node to the two candidate nodes.
Since by definition there are no loops in the tree and each edge is a cut edge, we get three components of connectedness.
Then we calculate the average number of visits (the sum of the number of visits on all nodes divided by the number of nodes) in each of the two connectivity components and go where the number is smaller.
If the average number of visits is close, namely, they differ by no more than a given $\alpha$ value, then we go to the random side.
Please, see [Figure 3](#S3.F3) for example of path, where each number in a list inside a node (blue circle) represents the number of steps when updates occur in this node.

### 3.4 Cores Building

After all indices are found by the search procedure described above, we can build all cores based on these indices.
First, consider a leaf node ${D}^{(L)}_{j}$, and let its down indices are $i^{\text{\scriptsize\tt down}}_{L,j}$ and the values of these indices are $v^{\text{\scriptsize\tt down}}_{L,j}$ (recall, that $\{j\}\cup i^{\text{\scriptsize\tt down}}_{L,j}=\{1,\,2,\,\ldots,\,d\}$).
Then the core, associated with this node is calculated as follows.
First, we form a matrix $V$ of values of the BB using these indices

$$ $\begin{split}V[i,\,k]=f(I_{ik})\;\text{ with}\quad I_{ik}[j]=i,\;\;\\ I_{ik}[i^{\text{\scriptsize\tt down}}_{L,j}]=v^{\text{\scriptsize\tt down}}_{L,j}[k],\;\;\forall 1\leq i\leq N_{j},\;1\leq k\leq{r}^{(L)}_{j},\end{split}$ $$

and then we let the core ${\mathcal{G}}^{(L)}_{j}$ be the transposed factor $Q$ of the QR-decomposition of this matrix $V$

$$ ${\mathcal{G}}^{(L)}_{j}=Q^{\text{T}},\text{ where }\{Q,\,R\}=\operatorname{QR}(V).$ $$

For the non-leaf and non-root node ${D}^{(l)}_{j}$ we perform a similar procedure.
Let $i^{\text{\scriptsize\tt down}}_{l,j}$ and $v^{\text{\scriptsize\tt down}}_{l,j}$ be its down indices and its down indices value, respectively.
Let $i^{\text{\scriptsize\tt up}}_{l+1,j_{c}}$ and $v^{\text{\scriptsize\tt up}}_{l+1,j_{c}}$ be upper indices and upper indices value, respectively, for the $c$th child of this node, where $c=1,2$ (recall, that $i^{\text{\scriptsize\tt up}}_{l+1,j_{1}}\cup i^{\text{\scriptsize\tt up}}_{l+1,j_{2}}\cup i^{\text{\scriptsize\tt down}}_{l,j}=\{1,\,2,\,\ldots,\,d\}$).
Then we first build the matrix $V$

$$ $V[\overline{in},\,k]=f(I_{ink})\;\text{ with}\\ I_{ink}[i^{\text{\scriptsize\tt up}}_{l+1,j_{1}}]=v^{\text{\scriptsize\tt up}}_{l+1,j_{1}}[i],\;\;I_{ink}[i^{\text{\scriptsize\tt up}}_{l+1,j_{2}}]=v^{\text{\scriptsize\tt up}}_{l+1,j_{2}}[n],\\ I_{ink}[i^{\text{\scriptsize\tt down}}_{l,j}]=v^{\text{\scriptsize\tt down}}_{l,j}[k],\\ \forall 1\leq i\leq{r}^{(l+1)}_{j_{1}},\;1\leq n\leq{r}^{(l+1)}_{j_{2}},\;1\leq k\leq{r}^{(l)}_{j},$ $$

and then we let the values of the core ${\mathcal{G}}^{(l)}_{j}$ be the “reshaped” values of the factor $Q$ of the QR-decomposition of $V$

$$ ${\mathcal{G}}^{(l)}_{j}[i,\,k,\,n]=Q[\overline{in},\,k],\text{ where }\{Q,\,R\}=\operatorname{QR}(V).$ $$

Finally, for the root node ${D}^{(1)}_{1}$, we let the values of the assigned core be the values of the given BB in the corresponding points.
Namely, let $i^{\text{\scriptsize\tt up}}_{2,j_{c}}$ and $v^{\text{\scriptsize\tt up}}_{2,j_{c}}$ be upper indices and upper indices value for the $c$th child of the root node, $c=1,2$.
Then for all $1\leq i\leq{r}^{(2)}_{j_{1}},\;1\leq n\leq{r}^{(2)}_{j_{2}}$ we have

$$ ${\mathcal{G}}^{(1)}_{1}[i,\,1,\,n]=Q[i,\,n],\;\;Q[i,\,n]=f(I_{in})\\ \text{ with}\quad I_{in}[i^{\text{\scriptsize\tt up}}_{2,j_{1}}]=v^{\text{\scriptsize\tt up}}_{2,j_{1}}[i],\;I_{in}[i^{\text{\scriptsize\tt up}}_{2,j_{2}}]=v^{\text{\scriptsize\tt up}}_{2,j_{2}}[n].$ $$

Note, that due to this procedure, the obtained cores are orthogonalized and, therefore, their maximum modulo values are moderated.

## 4 Related Work

In many practical situations, the problem-specific target function is not differentiable, too complex, or its gradients are not helpful due to the non-convex nature of the problem, and it has to be treated as a black box (BB).
In this case, two important problems naturally arise: approximation (Bhosekar & Ierapetritou, ) and optimization (Alarie et al., ).

The approximation carried out in the offline phase allows us to build a surrogate (simplified) model of the BB, which can then be used in the online phase to quickly calculate its values and various characteristics.
In the multidimensional case, it becomes difficult to construct a surrogate model, and low-rank tensor approximations are often the most effective.
Several recent works (Kapushev et al., ; Ahmadi-Asl et al., ; Chertkov et al., ) proposed various new algorithms based on the TT-decomposition for approximating high-dimensional functions.
If we have access to the BB and can perform dynamic queries, then the powerful TT-cross method (Oseledets & Tyrtyshnikov, ) is often used, and if only a training dataset is available, then the TT-ALS method (Holtz et al., ) is preferred.
In this work, we consider the case of adaptive queries to the BB, so we select the TT-cross method as the main baseline for the approximation problem.

Gradients are not available for the BB, so only gradient-free methods can be used for the optimization problem. Particle Swarm Optimization (PSO) (Kennedy & Eberhart, ) and Simultaneous Perturbation Stochastic Approximation (SPSA) (Maryak & Chin, ) are rather useful methods in this case.
There is also a large variety of other heuristic methods for finding the global extremum.
Recently, the TT-decomposition has been actively used for black-box optimization, since it turns out to be more effective than standard approaches in the multidimensional case.
An iterative method TTOpt based on the maximum volume approach is proposed in the work (Sozykin et al., ).
The authors applied this approach to the problem of optimizing the weights of neural networks in the framework of reinforcement learning problems in (Sozykin et al., ) and to the QUBO problem in (Nikitin et al., ).
A similar optimization approach was also considered in (Selvanayagam et al., ) and (Shetty et al., ).
One more promising algorithm, Optima-TT, which is based on the probabilistic sampling from the TT-tensor, was proposed in recent work (Chertkov et al., ).
We also note the work (Soley et al., ), where an optimization method based on the iterative power algorithm in terms of the quantized version of the TT-decomposition is proposed.
As a result, we consider classical PSO and SPSA methods as well as the TTOpt method as baselines for the optimization problem.

## 5 Numerical Experiments

**Table 1: Benchmark functions for performance analysis of the proposed method.**
| Function | Bounds | Analytical formula |
| --- | --- | --- |
| Alpine | [$-10$, $10$] | $\sum_{i=1}^{d}|x_{i}\sin{x_{i}}+0.1x_{i}|$ |
| Chung | [$-10$, $10$] | $\left(\sum_{i=1}^{d}x_{i}^{2}\right)^{2}$ |
| Dixon | [$-10$, $10$] | $\textsf{f}(\bm{x})=(x_{1}-1)^{2}+\sum_{i=2}^{d}i\cdot\left(2x_{i}^{2}-x_{i-1}\right)^{2}$ |
| Griewank | [$-100$, $100$] | $\sum_{i=1}^{d}\frac{x_{i}^{2}}{4000}-\prod_{i=1}^{d}\cos{\left(\frac{x_{i}}{\sqrt{i}}\right)}+1$ |
| Pathological | [$-100$, $100$] | $\sum_{i=1}^{d-1}\left(0.5+\frac{\sin^{2}\sqrt{100x_{i}^{2}+x_{i+1}^{2}}-0.5}{1+0.001\left(x_{i}^{2}-2x_{i}x_{i+1}+x_{i+1}^{2}\right)^{2}}\right)$ |
| Pinter | [$-10$, $10$] | $\sum_{i=1}^{d}\left(ix_{i}^{2}+20i\sin^{2}{A_{i}}+i\log_{10}{(1+iB_{i}^{2})}\right)$,<br>where<br>$A_{i}=x_{i-1}\sin{x_{i}}+\sin{x_{i+1}}$,<br>$B_{i}=x_{i-1}^{2}-2x_{i}+3x_{i+1}-\cos{x_{i}}+1)$<br>with $x_{0}=x_{d}$ and $x_{d+1}=x_{1}$ |
| Qing | [0 0, $500$] | $\textsf{f}(\bm{x})=\sum_{i=1}^{d}\left(x_{i}^{2}-i\right)^{2}$ |
| Rastrigin | [$-5.12$, $5.12$] | $A\cdot d+\sum_{i=1}^{d}\left(x_{i}^{2}-A\cdot\cos{(2\pi x_{i})}\right),$<br>where $A=10$ |
| Schaffer | [$-100$, $100$] | $\sum_{i=1}^{d-1}\left(0.5+\frac{\sin^{2}{\left(\sqrt{x_{i}^{2}+x_{i+1}^{2}}\right)}-0.5}{\left(1+0.001\cdot(x_{i}^{2}+x_{i+1}^{2})\right)^{2}}\right)$ |
| Schwefel | [0 0, $500$] | $-\frac{1}{d}\sum_{i=1}^{d}x_{i}\cdot\sin{(\sqrt{|x_{i}|})}$ |
| Sphere | [$-5.12$, $5.12$] | $\sum_{i=1}^{d}x_{i}^{2}$ |
| Squares | [$-10$, $10$] | $\sum_{i=1}^{d}ix_{i}^{2}$ |
| Trigonometric | [0 0, $\pi$] | $\sum_{i=1}^{d}\left(d-\sum_{j=1}^{d}\cos{x_{j}}+i(1-\cos{x_{i}}-\sin{x_{i}})\right)^{2}$ |
| Wavy | [$-\pi$, $\pi$] | $1-\frac{1}{d}\sum_{i=1}^{d}\cos{(kx_{i})}\cdot e^{-\frac{x_{i}^{2}}{2}}$,<br>where $k=10$ |

**Table 2: Approximation relative error for the HTBB and TT-cross applied to all considered $d=256$-dimensional benchmarks. The reported values are averaged over $10$ independent runs.**
| Benchmark | HTBB | TT-cross |
| --- | --- | --- |
| Alpine | 2.83e-15 | 1.73e-02 |
| Chung | 7.87e-03 | 2.86e-02 |
| Dixon | 5.65e-03 | 1.00e-01 |
| Griewank | 2.83e-15 | 1.43e-02 |
| Pathological | 3.92e-02 | 1.08e-01 |
| Pinter | 1.23e-02 | 1.47e-02 |
| Qing | 3.67e-02 | 4.87e-02 |
| Rastrigin | 1.01e-14 | 1.47e-02 |
| Schaffer | 1.87e-02 | 1.88e-02 |
| Schwefel | 3.39e-14 | 6.31e-01 |
| Sphere | 1.20e-14 | 1.44e-02 |
| Squares | 1.07e-14 | 1.77e-02 |
| Trigonometric | 2.76e-02 | 4.82e-02 |
| Wavy | 8.56e-05 | 2.46e-03 |

**Table 3: Approximation relative error for the HTBB applied to all considered $512$ and $1024$-dimensional benchmarks. The reported values are averaged over $5$ independent runs.**
| Benchmark | $d=512$ | $d=1024$ |
| --- | --- | --- |
| Alpine | 4.92e-15 | 3.81e-04 |
| Chung | 7.86e-03 | 7.64e-03 |
| Dixon | 3.75e-03 | 2.83e-03 |
| Griewank | 1.37e-14 | 3.16e-14 |
| Pathological | 3.80e-02 | 3.76e-02 |
| Pinter | 8.80e-03 | 8.38e-03 |
| Qing | 1.85e-02 | 1.60e-02 |
| Rastrigin | 1.63e-14 | 1.02e-04 |
| Schaffer | 1.94e-02 | 1.52e-02 |
| Schwefel | 2.59e-13 | 1.23e-13 |
| Sphere | 1.16e-14 | 4.58e-14 |
| Squares | 1.08e-14 | 2.38e-14 |
| Trigonometric | 2.74e-02 | 2.38e-02 |
| Wavy | 1.18e-04 | 3.38e-04 |

**Table 4: Minimization results for the HTBB, TTOpt, One+One, SPSA, and PSO applied to all considered $256$-dimensional benchmarks. The reported values are averaged over $10$ independent runs.**
| Benchmark | HTBB | TTOpt | One+One | SPSA | PSO |
| --- | --- | --- | --- | --- | --- |
| Alpine | 6.75e+01 | 4.48e+02 | 3.66e+02 | 3.99e+02 | 4.76e+02 |
| Chung | 1.45e+06 | 7.74e+07 | 1.48e+06 | 1.54e+06 | 6.98e+07 |
| Dixon | 1.89e+06 | 1.99e+08 | 2.33e+06 | 3.20e+06 | 2.68e+08 |
| Griewank | 3.11e+01 | 2.21e+02 | 3.19e+01 | 3.12e+01 | 2.09e+02 |
| Pathological | 6.97e+01 | 1.02e+02 | 1.14e+02 | 9.32e+01 | 1.06e+02 |
| Pinter | 5.17e+05 | 1.19e+06 | 5.67e+05 | 5.94e+05 | 1.51e+06 |
| Qing | 4.99e+06 | 2.98e+12 | 8.47e+10 | 1.26e+12 | 1.76e+12 |
| Rastrigin | 9.19e+02 | 3.61e+03 | 1.09e+03 | 9.31e+02 | 3.72e+03 |
| Schaffer | 9.87e+01 | 1.15e+02 | 1.06e+02 | 1.02e+02 | 1.20e+02 |
| Schwefel | -3.85e+02 | -1.77e+02 | -1.92e+02 | -1.89e+02 | -1.38e+02 |
| Sphere | 3.16e+02 | 2.30e+03 | 3.24e+02 | 3.17e+02 | 2.19e+03 |
| Squares | 1.55e+05 | 7.37e+05 | 1.57e+05 | 1.55e+05 | 1.02e+06 |
| Trigonometric | 8.72e+04 | 9.30e+06 | 2.62e+05 | 1.77e+07 | 1.01e+07 |
| Wavy | 3.19e-01 | 6.21e-01 | 3.64e-01 | 3.22e-01 | 6.36e-01 |

To demonstrate the effectiveness of the proposed HTBB approach, we select $14$ popular $256$-dimensional benchmarks (Jamil & Yang, ; Vanaret et al., ; Dieterich & Hartke, ), which correspond to analytical functions with complex landscape and are described in detail in [Table 1](#S5.T1).
For each benchmark, we fix the input dimension at $256$ and consider the approximation and optimization problem in the black-box settings for the tensor that arises when the corresponding function is discretized on a Chebyshev grid with $8$ nodes in each dimension.
In all cases, we limited the budget (the number of requests to the BB) to $10^{4}$, and the HT-rank was taken to be $2$.

### 5.1 Multidimensional approximation

For each $256$-dimensional benchmark we perform the approximation with the proposed HTBB method and compare it with the TT-cross method,(^5^55 We used the implementation of the TT-cross method from [https://github.com/AndreiChertkov/teneva](https://github.com/AndreiChertkov/teneva).) constrained by the same budget ($10^{4}$ requests to BB).
The relative $L2$ errors on test sets of $10^{4}$ random points which were generated for each benchmark are reported in [Table 2](#S5.T2) (the computations were repeated $10$ times for both methods and the averaged results are presented).
Also in [Figure 4](#S5.F4) we provide a graphical comparison of the results for two benchmarks for the case of different values of the problem dimension ($5$, $10$, $50$, $100$, $200$).
As follows from the presented results, for all problems our method turns out to be more accurate than the baseline, and in some cases its accuracy turns out to be many orders of magnitude higher.
For the case of higher dimensions for the considered problem classes, running the TT-cross method leads to failures in software implementation due to instability, while our approach remains stable and gives high accuracy, as follows from values reported in [Table 3](#S5.T3) for dimensions $512$ and $1024$.

### 5.2 Multidimensional optimization

For each $256$-dimensional benchmark we perform the optimization (namely the search for a global minimum) with the proposed HTBB method.
We consider as baselines the tensor-based optimization method TTOpt(^6^66 We used the implementation of the TTOpt method from [https://github.com/AndreiChertkov/ttopt](https://github.com/AndreiChertkov/ttopt).) and three popular gradient-free optimization algorithms from the nevergrad framework (Bennet et al., ):(^7^77 See [https://github.com/facebookresearch/nevergrad](https://github.com/facebookresearch/nevergrad).) One+One, SPSA, and PSO.
The limit on the number of requests to the objective function was fixed at the value $10^{4}$.
The calculations were repeated $10$ times for all methods and the averaged results are presented in [Table 4](#S5.T4).
Also in [Figure 5](#S5.F5) we show the convergence plots for two benchmarks.
As follows from the reported values, HTBB, in contrast to alternative approaches, gives a consistently top result for all model problems.

Figure: Figure 4: Approximation results for Alpine and Dixon functions for cases of dimensions $5$, $10$, $50$, $100$, and $200$. For both methods, we plot the relative error of the solution averaged over $10$ runs with a solid line and fill in the area between the worst and best result with the same color.
Refer to caption: /html/2402.02890/assets/figures/FuncAlpine_appr.png

Figure: Figure 5: Minimization results for Alpine and Dixon functions. For each of the optimizers, we plot the value of the solution averaged over $10$ runs with a solid line and fill in the area between the worst and best result with the same color.
Refer to caption: /html/2402.02890/assets/figures/FuncAlpine_opti.png

## 6 Conclusions

In this work, we presented a new method HTBB for simultaneously solving the problem of multidimensional approximation and gradient-free optimization for functions given in the form of a black box.
Our approach is based on the low-rank hierarchical Tucker decomposition, which makes it especially effective in the multidimensional case.
The key features of the presented work are a) using the MaxVol algorithm which allows efficiently finding the required indices
and b) using the sequential traversal of cores,
allowing to move to one of the neighboring nodes and making it more efficient to find indexes that need updating.

The HTBB method can be applied to a wide class of practically significant problems,
including optimal control and various machine learning applications.
As future work, we point out the possibility of a rather simple extension on the HT-structure of the algorithms that now exist for the TT-decomposition: rounding, orthogonalization, search for the maximum element by the top-k-like methods, etc.

## Peer Reviews

confidence: 3
contribution: 2
initial_rating: 3
presentation: 2
questions: - typo on page 5, line 47: $i^{down}_{l,j}$
- typo on page 6, line 323: "As each visited node"
soundness: 3
strengths: - The proposed algorithm achieved better numerical performances than the prior works.
- The authors did a good job explaining the details of the proposed algorithm despite the complicated notations, which are common in tensor-related work.
summary: This paper develops new black-box approximation and gradient-free optimization methods for multidimensional functions based on Hierarchical Tucker (HT) decomposition. Numerical experiments were conducted to illustrate the improvement over the prior method that was developed based on Tensor Train (TT) decomposition.
weaknesses: - A key concern is the lack of theoretical guarantees. Compared to the existing algorithms like TT-Cross and TTOpt, there is no clear improvement in terms of time complexity, approximation guarantees, or optimization guarantees with respect to the dimension $d$ and rank $r$. The algorithm seems too heuristic to me and the contribution is not enough for an ICLR publication.
- In addition, I have reservations about the practicality of the proposed "black-box" approach. In the numerical experiments, all functions have well-behaved analytical expressions, which may not reflect real-world scenarios. It remains unclear whether the approach would perform well in applications where functions may lack favorable properties, such as continuity or convexity.
- The presentation of the paper, particularly in terms of English language usage, could be improved.

---

confidence: 4
contribution: 2
initial_rating: 3
presentation: 1
questions: - For the 14 benchmark instances, [Jamil-Yang, J. of Math. Modeling and Num. Opt. 2013] has 100+
  tensors in their dataset. How did you choose this subset of 14 across multiple works?
- "gradient-free" methods are mentioned regularly throughout the paper. Can you give an example
  of an effective method that you compare against that is gradient-based?
- Are the sizes (i.e., the number of optimizable parameters) in HTBB and TT-cross equivalent in Table 1?
- In Figure 4, why are the SPSA and OPO methods better up until $10^3$ requests?
- Note: The [Goreinov et al., 2010] citation should include the journal.
soundness: 2
strengths: - The main technical details in Section 4 are sophisticated and well organized.
  That said, they could benefit from more context and a better setup in the
  preceeding sections.
- Experiments are run against well defined benchmarks tensors from previous
  work, e.g., [Jamil-Yang, J. of Math. Modeling and Num. Opt. 2013]. These
  tensors are discretized continuous functions with known ground truth.
- Table 3 and Figure 4 in the experiments suggest that HTBB has a lot of
  potential, e.g., there are clear improvements in the $10^4$ number-of-requests
  regime.
summary: This work proposes a method called Hierarchiacl Tucker Black Box (HTBB) based
on low-rank hierarchical Tucker decomposition and the MaxVol subroutine for
finding good submatrices [Goreinov et al., Matrix Methods: Theory, Algorithms
and Applications 2010]. The approach is inspired by the algorithms:
* TT-cross [Oseledets-Eugene, Linear Algebra and Its Applications, 2010]
* TTOpt [Sozykin et al., NeurIPS 2022]

with the goal of extending these works
to hierarchical Tucker tree structures. The main techincal contribution of
this work is Algorithm 2, which decomposes dimensions of the tensor to
construct a good low-rank hierarchical Tucker tree structure. Finally, this
work compares HTBB with several other algorithms across a benchmark of
high-dimensional tensors representing discretized continuous functions.
weaknesses: - First and foremost, the presentation quality of this work makes it
  challenging to (1) have a precise understanding of what problem is being
  solved, and (2) the difference between your approach and other existing works.
  The good ideas in this work are getting lost, largely because there is not a
  crisp and standalone mathematical problem definition (relies too much on
  pointing to related work).
- Related works cover many different method types (e.g., particle swarm
  optimization).  How many of these solve the same problem and could be
  compared against in the experiments?

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: 1. The introduction of the equation in Section 3 appears abrupt. It is recommended that the authors include a figure to provide further illustration and clarity.

2. What would be the performance of Algorithm 2 if QR decomposition were replaced with singular value decomposition (SVD)?

3. Given the existence of a canonical form for the hierarchical Tucker decomposition, could the authors provide a comparison between the cross-approximation method and an SVD-based method in the context of hierarchical Tucker?

4. In the first experiment, in addition to the TT cross-approximation, the Tucker cross-approximation should be considered, as hierarchical Tucker is an extension of the Tucker decomposition.

5. The experimental results should include performance comparison with regard to different tensor orders, as this would provide insight into the scalability of the proposed method.

6. All experimental results should be normalized for comparability.
soundness: 3
strengths: The authors propose a low-rank hierarchical Tucker cross-approximation method, utilizing the MaxVol index selection procedure. The paper is well-structured, with a comprehensive introduction. In my assessment, this work constitutes a notable contribution to the field.
summary: In this paper, the authors propose a low-rank hierarchical Tucker cross-approximation method, utilizing the MaxVol index selection procedure. The paper is well-structured, with a comprehensive introduction. In my assessment, this work constitutes a notable contribution to the field.
weaknesses: Please refer to the questions below.

---

confidence: 2
contribution: 2
initial_rating: 5
presentation: 2
questions: Q1) I believe the paper would be considerably improved by taking into accound the point W1 stated above. In particular, I suggest the following:
- Briefly recall the rationale and the mechanics of TT-Cross and TTOpt.
- Write Algorithm 1 in a more explicit fashion, clearly stating its inputs and outputs, using a more formal notation.

Q2) In Section 3, the authors state that by HT they "mean a tensor tree that is not necessarily balanced". Thus, the topology of the HT structure used in the numerical experiments should be specified. Is it always taken to be simply a balanced tree? 

Q3) Related to the previous question, how is this topology chosen in practice? Is it simply by trial-and-error or are there some heuristics for doing so? 

Q4) Intuitively, the higher flexibility of the proposed scheme in comparison with the TT-based methods should lead to smaller approximation errors, provided that the HT structure is correctly chosen. Hence:
- Is it possible to derive any quantitative comparison between the TT and HT approaches?
- More generally, can you give any approximation guarantees of your method for some reasonable classes of functions, say as function of the number of point evaluations? 

Q5) Regarding the point W3 raised above, is the novelty here all in the way where HTBB selects the next node that should be visited, the rest of it being similar to TT-Cross and TTOpt?

Q6) Some details and relevant discussion are lacking in the numerical experiments. Namely, what are the TT ranks used in the experiments of sections 6.1 and 6.2? What are the total number of parameters in both models? If the TT ranks are increased (but the HT ones are held fixed), is the TT model able to gradually bridge the performance gap, and if so, at which additional cost?
soundness: 3
strengths: S1) The paper generalizes tensor-tain-based function approximation algorithms to more general hierarchical Tucker structures.

S2) Its simulation results convincingly illustrate the benefits of the proposed approach in terms of function approximation and maximization performance, presumably thanks to the higher degree of flexibility of HT in comparison with TT.
summary: This paper is concerned with the approximation and maximization (or minimization) of a real-valued function defined on a subset of $\mathbb{R}^d$. The function is considered as being a "black box", meaning its properties are unknown, but one can evaluate it over a predefined discretization grid, implicitly defining a tensor (of point function evaluations). Specifically, the authors propose extensions of the TT-Cross and TTOpt algorithms, which rely on the so-called tensor train (TT) tensor format, to a more general hierarchical Tucker (HT) format that is specified by a tree topology which indicates how the cores (factors) involved in the low-rank approximation are connected. The main technical difficulty seems to be how to navigate that tree for selecting a subset of indices used for building the approximation.
weaknesses: W1) Despite the proposed algorithms being potentially useful for practitioners, in my view the contribution of this paper is fairly incremental, as it extends some well-known algorithms to a different tensor format, but without introducing new ideas, insights or relevant theoretical results.

W2) The paper is fairly hard to read for non-experts, which I believe will be the case for the majority of the ICLR audience. In that regard:
- A brief description of the rationale behind the TT-Cross and TTOpt should be given, so that the readers understand what the authors are trying to implement. This will make the paper minimally self-contained, since the extension of these algorithms is the main goal here. 
- Related do that point, is not quite clear what is the output of Algorithm 1. It seems to be the cores that will define the approximation, but at the same time this algorithm is expected to yield the max (or min) of the function, as announced in the paper.

W3) The paper should also clearly state exactly where the novelty is in terms of the introduced algorithms. 

W4) A careful revision is needed for improving English usage. Examples:
- "possible maximum rank increment" -> maximum possible rank increment
- "then we go to the random side" -> then we randomly choose one side to go to.
- "see Figure 3 for example of path" -> for an example of a path
- "For the core of the root node it is hold"  (unclear what is meant)

Meta-review: This work proposes to approximate a black-box function of discrete inputs with hierarchical Tucker decomposition. A new heuristic algorithm is presented, but no theoretical guarantee is provided. Numerical experiments are conducted on some well-behaving functions, therefore not particularly convincing about its universal approximation ability.

Readability could be improved.

Decision: Reject