Sunny Katyara, Fanny Ficuciello, Darwin G Caldwell, Fei Chen, Bruno Siciliano
11institutetext: Department of Advanced Robotics, Istituto Italiano di
Tecnologia, Via Morego 30, 16163, Genova, Italy
22institutetext: Interdepartmental Center for Advances in Robotic Surgery,
University of Naples Federico II, Naples 80125, Italy

# Reproducible Pruning System on Dynamic Natural Plants for Field Agricultural
Robots

Sunny Katyara 1122 Fanny Ficuciello 22 Darwin G. Caldwell 11 Fei Chen 11 Bruno
Siciliano 22

###### Abstract

Pruning is the art of cutting unwanted and unhealthy plant branches and is one
of the difficult tasks in the field robotics. It becomes even more complex
when the plant branches are moving. Moreover, the reproducibility of robot
pruning skills is another challenge to deal with due to the heterogeneous
nature of vines in the vineyard. This research proposes a multi-modal
framework to deal with the dynamic vines with the aim of sim2real skill
transfer. The 3D models of vines are constructed in blender engine and
rendered in simulated environment as a need for training the network. The
Natural Admittance Controller (NAC) is applied to deal with the dynamics of
vines. It uses force feedback and compensates the friction effects while
maintaining the passivity of system. The faster R-CNN trained on 3D vine
models, is used to detect the spurs and then the statistical pattern
recognition algorithm using K-means clustering is applied to find the
effective pruning points. The proposed framework is tested in simulated and
real environments.

###### keywords:

Dynamic plants, Sim2Real transfer, Interaction controller, Deep learning,
Pattern recognition

This research is supported in part by the project “Grape Vine Perception and
Winter Pruning Automation” funded by joint lab of Istituto Italiano di
Tecnologia and Università Cattolica del Sacro Cuore, and the project
“Improving Reproducibility in Learning Physical Manipulation Skills with
Simulators Using Realistic Variations” funded by EU H2020 ERA-Net Chist-Era
program

##  1 Introduction

Field robotics is an important research domain especially in agricultural
applications as promising solution to increase the yield and reduce the
environmental effects. Field robotics uses advanced technologies to give more
autonomy to mobile robots to help or do the tasks that are hazardous and
dangerous to humans. Now-a-days, most of the agricultural robots are primarily
being used for monitoring and quality control such as aerial inspection but
can also be used for heavy tasks like harvesting, seeding, pruning etc. As an
emerging research domain for agricultural robotics, robotic pruning system is
receiving much attention in the community these days. Pruning is an important
agricultural activity of removing the branches, limbs and spurs from the
natural plants. Pruning becomes effective if it only removes the unwanted,
unhealthy and poorly positioned branches with least effects on the healthy
parts. The effective pruning makes the difference to the quality and yield of
plants in the field [1]. Pruning which requires qualified and experienced
professionals, can eventually be performed by intelligent and autonomous
robotic system. The stand-alone mobile robot equipped with vision and machine
learning is used to travel through the agricultural field and to identify,
manipulate and finally cut the branches. The general idea of pruning is
illustrated in Fig. 1. As a representative application for commercial plants,
the grape vine winter pruning is chosen for analysis. It can be seen from Fig
1a that grape vine pruning is more difficult as compared to other plants,
because of complexity in manipulating, locating and deciding the potential
cutting points.

Figure 1: Pruning schema (a) represents general pruning scenarios (b)
illustrates the robot pruning system.

Manipulation skills of robot require active perception and interactive
learning technologies. Active perception defines the behaviour of agents i.e,
robots in the environment on the basis of collected sensor data in order to
increase its autonomy and performance. However, robot learning is an iterative
process and thus the robot learner tries eventually to learn the specific set
of skills and generalize them to new conditions [2]. Since, in agricultural
applications, the reproducibility of manipulations skills is another important
aspect because most of the plants even in the same field have different
physical characteristics. Therefore, a simulator with real time rendering
allows developing large data base and testing the learned skills under the new
conditions. Moreover, the sim2real interface allows teaching and interacting
with the robots in their virtual environments (simulated) and enables them to
learn different manipulation skills in the multi-structured scenarios. In
order to create the virtual environments for robots to learn, the 3D modeling
of objects is one of the key steps towards the object manipulation.

However, manipulation in field robotics has not yet achieved maturity because
of uncertainties associated with objects and environments. Most of the machine
learning techniques consider fully known, discrete, stationary, deterministic
and discrete environments for robot manipulation skill learning and transfer
[3]. But in order to deal with dynamic, unknown, stochastic and continuous
environments, more advanced learning approach with promising long-term
solution has to be studied, e.g., Reinforcement Learning (RL) [4]. Due to the
number of issues such as natural environmental constraints, limited data sets,
need for repetition of experiments, limited resources etc the sim2real
approach with highly efficient rendering capabilities are chosen to teach the
specific set of skills to robots in simulated environments.

A lot of research has already been done on manipulating rigid body objects but
the manipulation of soft and moving objects is still an open question to
answer [5][6]. There are many hindrances in automating the manipulation of
soft objects and the most common among them is the estimation of their
deformation properties because even they are not same among the objects of
similar groups. Moreover, in order to deal with moving objects, especially
plants, the interaction controllers are needed [7]. The most popular
interaction controllers are impedance and admittance controllers, which offer
gentle and fast response. Impedance controller under the large variations of
system parameters such as compliance, stiffness, damping and inertia affects
the passivity and stability of system [8]. However, the admittance controller,
more precisely Natural Admittance Controller (NAC) proposed by [9] uses the
force feedback to overcome the effects of friction during the interaction
without affecting the passivity of system. In our case, NAC is one of the
feasible solution to deal with moving plants whose dynamics are unknown during
the interaction and thus requires the passivity of system to be maintained.
The complex task of pruning has several steps to be followed but the two most
important among them are; the perception of agricultural field and
identification of spurs on the plants to be pruned and the second is related
to motion planning of robot, with the intention of making cuts.

To the best of our knowledge, most of work on robotic pruning has been done on
stationary branches and with least consideration on reproducibility of robotic
pruning skills. The contributions of this research are; (1) proposing multi-
modular framework which combines deep learning (Faster R-CNN) , computer
vision (graph morphometry) and robot control (natural admittance controller)
for dynamic and reproducible pruning skills (2) Designing 3D models of the
grape vines similar to real ones in blender for training the network in the
simulator (3) Teaching and testing pruning skills on mobile manipulator in
simulated environment initially and then evaluate them on real system. For the
sake of convenience, the grape vine winter pruning is taken as an
representative case study, shown in Fig. 1(b).

##  2 Related Works

A systematic literature review on the challenges and solutions for effective
pruning, sensing and automation technologies is investigated in [10]. It
points out that the effective pruning relies on the set of key technologies:
pruning methods, 3D modelling of dynamic branches, training of robot
perception, and handling the dynamics of process. Even though the study is
limited to apple tree pruning but its discussion inspires for general natural
plant pruning. Another design of robot system for harvesting apple orchard is
presented in [11]. This system uses a redundant robot manipulator equipped
with global camera, which can effectively harvest the desired area at a least
error rate using an open loop controller. However, the reliability of system
need to be improved with grasping status force feedback. A robot system is
introduced in [12] to address the grape vine pruning. The robot system is
equipped with multiple cameras for online 3D reconstruction of the vines.
Having this geometrical information of vines, the robot plans a collision free
trajectory to apply reasonable cuts. However, the research does not take into
account the dynamics of vines and their setup uses a tractor pulling a house
with robot system inside the vineyard, limiting its wide adoption to different
complex conditions. A dual-arm mobile manipulator for pruning and fruit
picking is discussed in [13]. It deals with the non-holonomic constraint
introduced by the wheeled platform using recursive Gibbs-Appell method. The
system is tested under symmetrical and unsymmetrical conditions during the
pruning and fruit picking. The system considered the dynamics of rigid objects
but is not tested on moving and soft object manipulation. The effectiveness of
algorithm is measured in terms of simplicity, computational timing and
repeatability. The idea of sim2real skill transfer for the simple manipulation
tasks i.e, object pushing is suggested in [14]. The idea is based on
randomizing the dynamics of simulator to teach the recurrent policies to robot
so that it can adapt to different environments in real world. However, this
approach can be extended to the manipulation of complex tasks i.e, pruning
with the incorporation of deep learning and vision system.

Figure 2: Block diagram of proposed methodology.

##  3 Research Methodology

For the two important sub-tasks discussed in previous section, initially a
convolutional neural network (ResNet-101 R-CNN) trained on modeled 3D vines,
is used to detect the possible spurs on the analyzed branches from the set of
RGB images and then the graph morphometry together with K-means is applied to
find the Potential Pruning Points (PPP) on the vines. Once the pruning points
are determined, the next task is to command the robot towards the respective
branches on the vines with desired trajectories and hence the quadratic
programming is being used for such path planning [15]. The block diagram of
proposed methodology is shown in Fig. 2 and associated kinematic model in Fig.
3.

Figure 3: Kinematic model of robot pruning system.

###  3.1 Modeling of 3D vines and robot platform

The 3D models of vines are constructed in blender engine and then rendered
into simulated environment. The 2D images of vines captured from agricultural
field are used as reference for creating 3D models. The 2D images are imported
into blender as background in vector graphics form and resized upto 76.3 cm to
give them more realistic shape. In the edit mode, the circular extruder with
3.121 cm diameter is selected and the path is traced along the 2D image but
with precise rotation, resizing and extension to take into account all the
spurs, canes and buds on the vine branches. Once the extrude is done, the
proportional editor is activated to give the tilt to the branches according to
their real shape. Moreover, the sub-division surface modifier is applied to
shape the irregularities on the vines. In order to give real looking texture,
UV mapping is done with unwrapping on the selected texture image to follow the
lines on real vines and cubic projection is applied to cover the whole surface
area of designed vine. Since, the vines are given multi-textures similar to
realistic ones and the buds are carefully being placed on alternatives sides
of the canes. For adding the softness to vines, specially the jiggle effect
when wind blows, the dynamics of vines are changed from rigid body to cloth
physics. The pinning is activated and the weighted paint is selected to add
jiggle to the branches and then a vertex group is created and scene is
activated. In the add brush setting, the paint is applied on the parts in wire
frame mode which are not needed to jiggle until they turn reddish. The final
result for all the chosen three vines, designed and rendered in blender are
shown in Fig 4 (a, b, c).

In this study, the mobile manipulator, named “Rolling Panda”, is consist of
7-DOF robot arm from Franka Emika, a two wheel mobile platform from Neobotix,
two Intel Realse D435i RGBD cameras attached to the end-effector, a cutting
tool with pruning shear attached to the end-effector, and a SICK laser scanner
for 2D SLAM.

Figure 4: 3D models of three different grape vines, (a) illustrates the
simpler vine with five spurs, (b) is the vine with six spurs and complex
geometry, (c) represents the more complex vine with interconnected shoots on
several individual spurs.

###  3.2 Estimation of Potential Pruning Points (PPP)

Before estimating the pruning points, it is required to find the dedicated
spur regions so as to narrow down the search space for pruning algorithm.
Identification of possible spurs on the vines is the crucial step in finding
the pruning region. The Faster R-CNN is used to mathematically characterize
the spurs [16]. It calculates the centroid of each region and assigns bounding
box around to ensure the average variance does not go beyond the centroid
boundaries. Such information about the centroid of spur is used to find its
pose. R-CNN model with ResNet-101 maps the spurs on the vines with mean
precision of 0.679 and at a learning rate of 0.0005. The architecture of
faster R-CNN network used for spur detection is shown in Fig. 5.

Figure 5: Architecture of faster R-CNN for spur detection.

The faster R-CNN uses convolution layers to extract the desired feature map
and the anchors are used to define the possible regions from the feature map.
Total 9 anchors are used with three aspect ratios 1:1, 1:2 and 2:1 which
stride over the feature map skipping 16 pixels at each time. The possible
regions predicted by anchors are fed to the regional proposal network (RPN).
The RPN is a small neural networks that predicts the region of the object of
interest. It actually predicts whether there is object (spur) or not and also
frames the bounding box around those spurs. In order to predict the labels of
anchors, the output from RPN is passed through softmax or regressor activation
functions. Since, RPN results in CNN feature maps of different sizes and hence
Region of Interest (RoI) pooling is applied to harmonize the given maps to the
same size. The faster R-CNN uses original R-CNN network to perform RoI
pooling. It takes the feature vector predicted by RPN, normalize it and passes
through two fully connected convolution layers with ReLU activation and thus
predicts the class and bounding boxes of object. In order to train and test
our faster R-CNN network, 1210 and 245 images of modeled grape vines are used
respectively. The dimensions of feature map generated are $25{\times}25$ and
the size of anchor is 9, therefore there are total 2916 potential anchors. The
test time per image is 0.2 secs and 8359 bounding boxes are detected in the
regions of interest with a minimum threshold of 0.7 and which are then finally
passed to classifier layer as an output.

Once the spurs are detected, the next step is to find the Potential Pruning
Points (PPP) which mainly depend upon the geometric information of vines and
their intrinsic properties. A predictive HOG (Histogram of Gradients) is
applied to find out the normalized patches in the spur regions. The regions
are subdivided into cells on the basis of gradient intensity and histogram of
gradient directions is computed as shown in Fig. 6. The pipeline adopted to
estimate the pruning points is defined in Fig . 7.

Figure 6: predictive HOG applied to spur region, (a) represents the spur
region of interest, (b) illustrates the binary information of region , (c)
denotes the predictive histogram of selected spur region.

Since, the algorithm results into a flattening matrix with dense information
about the intensity gradients. Therefore,the Mean Predictive HOG (mpHOG) is
calculated under a Bayesian Linear Regression for inference matrix
distribution which is normalized by a probabilistic model evidence. On the
basis of mpHoG descriptor, a statistical-based pattern recognition pipeline is
applied to identify the geometric properties of vines and is merged with
K-means clustering to classify the regions in terms of similarity of their
intrinsic characteristics. Since, the K-means is a unsupervised machine
learning technique used to find the centroid of a given data set and locates
the average of data around it [17]. The algorithm continues to minimize the
distance between consecutive centroids until it does not move any more. The
algorithm is formulated as minimization problem and is given by Eq. 1

$\bar{F}=\sum_{m=1}^{k}\sum_{n=1}^{j}{||U_{m}-V_{n}||^{2}}$ (1)

The objective function minimizes the distance between $U_{m}$ points and
$V_{n}$ cluster centre. The structural graph based morphometry uses the graph
kernel to identify the similarity between the obtained clusters. The
classifier with local information locates the centre of data in each
neighborhood, thus provides the information that’s maps the series of pruning
points in the given regions. Hence, the cutting pose based on the determined
pruning points is obtained by using $\bar{X_{p}}={R}{t}$. Where R and t
denotes the rotation and translation of pruning region with the given set of
pruning points. According to Rodrigues rotation formula, with the known unit
rotation axis vector r($r_{x}$ , $r_{y}$, $r_{z}$) and the angle of rotation
Ø, the R is expressed by Eq. 2

$R(\phi,r)={(cos{\phi})}{I}+{(1-cos{\phi})}{rr}^{t}+{sin{\phi}}{[r]_{x}}$ (2)

I is the identity matrix of $3\times 3$ and t is the skew symmetric matrix.
Let $p_{1}$ and $p_{2}$ be the two pruning points in the region, they should
satisfy the constraint defined by the Eq. 3

${p_{2}}^{T}{X_{p}}{p_{1}}=0$ (3)

The homogeneous vectors of the pruning points $p_{1}$ and $p_{2}$ are always
being pre-multiplied by the inverse of intrinsic camera matrix Mc.

Figure 7: Pipeline for estimation of pruning points.

###  3.3 Motion Planner

Motion planners are used to navigate the robot to a goal through the sequence
of valid configurations. The configurations of robot subjected to constraints
are used to avoid collisions with the obstacles and maximize the
manipulability of system. The motion planner formulates the task as a cost
function to be minimized on the basis of constraints solved by using the
equality and inequality solvers and is expressed by the Eq. 4.

${{X}\otimes{Y}}=\\{ax-y\in{R}^{n}|x\in X,y\in Y,a\in R\\}$ (4)

Where X, Y and R represent goals, constraints and robot configuration space
respectively. In order to define the pruning task and to set the priorities on
it, a stack is formed. In our research, the Stack of Task (SoT) is formulated
by the Eq. 5[18], with $L_{jointposition},$ and $L_{jointvelocity}$
representing constraints on robot system joint positions and velocities
respectively

$\begin{bmatrix}\\{T_{base}\\}^{World}\\\
\\{T_{Arm}\\}^{Base}+\\{T_{Endeffector}\\}^{Arm}\\\
\\{T_{Postural}\\}^{world}\end{bmatrix}$ (5)

s.t

$L_{jointposition}+L_{jointvelocity}$

T represents task and is generally defined by
${T_{i}}=({J_{i}},{\dot{e_{i}}})$

$J_{i}$ \- Jacobian of associated links:
$\dot{e_{i}}$ \- Error in joint velocities

The tasks are given hard and soft priorities and the constraints are defined
by the quadratic programming. The task has hard priority if it does not
deteriorate the performance execution of first task while the soft priorities
are defined at same task levels and they may influence each other on the basis
of their associated weights. In Eq. 5, the postural task is set at minimum
priority as compared to base, arm and end-effector levels. The tasks related
to the arm and end-effector are being solved without affecting the
configuration of the mobile base.

In order to deal with such kind of motion planning, the motion planner using
Quadratic Programming (QP) optimization interface is used[19]. The general
task using quadratic programming in joint space is described by the Eq. 6 and
is subjected to $R\dot{q}\leq{C}$ constraints

${\dot{q_{i}}}=argmin({J_{i}}{q_{i}},{\dot{e_{i}}})$ (6)

###  3.4 Natural Admittance Controller and dynamics of vines

When dealing with vibrating and moving objects, the friction during the
interaction is always ambiguous and passivity of the system is also lost.
Natural Admittance Controller (NAC) is the method to analyze the system under
the force feedback that results in the loss of passivity due to reduction of
inertia during the interaction [20]. Since, NAC declares the designated values
of inertias close to the actual physical values and thus ensures to maintain
the passivity by eliminating the friction effects. From the differential
kinematics, the relationship between the joint and operational space
velocities is given by ${\dot{X_{i}}}={J_{robot}}{\dot{q_{i}}}$. Where
$\dot{X_{i}}$ is $7\times 1$ vector of velocities of franka arm in operational
space and $\dot{q_{i}}$ are its associated joint velocities and $J_{robot}$ is
the Jacobian of franka arm.

The forced applied by the end-effector on the plants at the goal position is
determined by Eq. 7

${F_{e}}={K}({P_{g}}-{P_{base}})-{B}{\dot{X_{i}}}$ (7)

Where, K and B define the stiffness and damping of panda arm respectively,
$P_{g}$ and $P_{base}$ are the goal points and the base points respectively.
The total interaction force is the sum of end-effector force and feedback
force from the F/T sensor for plants ($F_{net}=F_{e}+F_{sensor}$). The
equivalent joint torque thus is computed by ${\tau_{i}}={J_{robot}}{F_{net}}$,
with obvious meaning of variables. The values of joint accelerations needed to
reach the goal are given by Eq. 8

${\ddot{q}_{desired}}={H}^{-1}{\tau_{i}}-{B_{H}}{\dot{q}_{actual}}$ (8)

Where, H is the matrix of inertias and $B_{H}$ is the desired damping. With
the use of Euler’s integration rule, the joint velocities are computed.

With the control loop running at $100Hz$, the desired joint velocities are
sent to command the robot to interact with stationary and moving vines. The
filter is used at force-feedback to suppress the noisy data and also the anti-
wind up criteria is applied to avoid the saturation of controller and system
parameters.

The motions of the vines are modeled with three virtual joints as PRR
kinematic chain in simulated environment, which can translate in and rotate
around $y$ direction and also generates rotation around $z$-axis. The dynamics
of plant is defined by ${F_{plant}}={m}{\dot{v}}+{b}{v}$

Where, ${F_{plant}}$ is the force applied by plant and its sensed by F/T
sensor as ${F_{sensor}}$, m represents the mass of plant, $\dot{v}$ is its
acceleration and b is the friction. The desired admittance of the system in
Laplace domain is defined by Eq. 9

${C_{desired}(S)}=\frac{1}{{m}{s}+{B}+{K/s}}$ (9)

Since, the mass of vine is not altered but the friction b is compensated to
get the desired behavior with the feedback gain given by Eq. 10 and the effect
of friction further diminishes with an increase of velocity loop gain $g_{v}$.
It has been found that the system under NAC is passive, no matter how large
the value of $g_{v}$ is increased for minimizing the friction effects.

${g_{f}}=\frac{({b}-{g_{v}}-{B}){s}-{K}}{{m}{s}^{2}+{B}{s}+{K}}$ (10)

##  4 Results And Discussion

Figure 8: Robot pruning scenarios in simulated and real environments, (a and
g) represent the mobile manipulator navigating the field, (b and h) illustrate
the pose of mobile manipulator reaching the designated vine, (c and i), (d and
j), (e and k), (f and l) report the first, second, third and forth cutting
actions respectively.

The results from the spur detection algorithm trained for almost 6 hours are
shown in Fig. 9(a, b). The bounding boxes represent the possible spurs on the
vines with confidence interval of more than $70\%$ in Fig. 9(b). Figures 9(c,
d) illustrate the Potential Pruning Points (PPP) obtained with the application
of proposed algorithm. Fig 9(c) represents the morphometry data of vine. The
big circle in Fig 9(d) on the regions indicates the classification of spurs,
either to prune or not and yellow dots in each region indicate the range of
pruning points detected. The pink dots denote the effective PPP on the vine to
perform the required cuts.

The cutting pose shown by orange line in the Fig. 9(d) between pruning points
is sent to the motion planner to generate the desired trajectories for the
robot platform, executed within 60 secs approximately. The application of
proposed framework on robot pruning scenario in simulated and real
environments is reported in Fig. 8. It can be seen from the Fig. 8 that the
results in the simulated environment Fig. 8(a-f) are successfully reproduced
on the real robot Fig. 8(g-l) with minor deviations in pruning time and
pruning accuracy. Total two spurs, 8 potential pruning points and four cutting
poses are detected within 6 secs in both the environments successfully. Since,
during the experiments on the real robots it is found that it skips one cut on
the designated vine position while the robot in simulated environment performs
all the four cuts more accurately and in a less amount of time. It is obvious,
because the changing real time conditions i.e, lighting, temperature, sensor
fidelity i.e, F/T information, friction of vines etc are however steady in the
simulated vineyard scenario. However, the domain randomization plays an
important role in bridging the gap between sim2real skill transfer but still
the $100\%$ reproduciblity of manipulation skills on real system suffers
mainly due to inaccuracy in modeling of system dynamics and varying
environmental conditions as what happens in our case too.

Figure 9: Detection of spurs and pruning points on vines. in (a) the chosen
vine for pruning is shown, in (b) the possible spurs on the selected vine are
found, in (c) the geometrical data of vine within the spur region is shown (d)
the location of pruning points, region and cutting pose are illustrated.

The associated joint position trajectories for the base and arm are shown in
Fig. 10 (a, b) respectively. It can be seen that the pose of robot arm is not
affected when the mobile base is navigated to a desired vine under 40 sec time
interval and then the robot arm is actuated to reach the cutting region
without violating the current configuration of mobile base. This case
validates the definition of hard constraints between the mobile base and robot
arm used in our framework.

Figure 10: Whole body motion control of robot platform. (a) shows the
trajectories followed by base to reach the desired location, in (b) the
trajectories to reach selected vine for the franka arm are shown, (c)
illustrates the pruning action by cutting blades attached to end-effector.

However, the soft constraints are defined between the cutting tool and robot
arm in order to facilitate the adaptive translation and orientation to reach
the desired cutting pose on the vines. The vines are initially stationary and
the cutting tool after reaching the pose of PPP, performs the required cut and
is shown in Fig. 10(c). The tool centre point (TCP) is defined on the cutting
tool at 5 mm away from tip of the blades. It can be seen that the blades on
the cutting tool moves $0.075$ m close to each other to execute the pruning
action. Next, the vines are given motion in different directions i.e,
$y$-translation, $y$-rotation, $z$-rotation and finally combination of all
(mixed). The motion of plants within all defined directions is shown in Fig.
11(a).

Figure 11: Dynamics of vine branches. (a) denotes designated motions of vines,
(b) shows the force profile of vines in static and moving conditions.

Due to the interaction between the robot arm and vine, the force is exerted by
both on each other. The force exerted by end-effector on plant is constant
regardless of the motion of vines. The force profile of vine in stationary and
dynamic conditions is shown in Fig. 11(b). Since the vines are placed in
Y-direction and when they touch the cutting tool during the motion equipped
with F/T sensor, it experiences a spike of impact as shown in Fig. 11(b). The
force applied by the end-effector $F_{e}$ is composed of two terms, one is the
force of spring between the arm and vine and the other is due to the damper
(inertia) of arm acting between reference and actual velocities of the end-
effector.

In case of dynamic pruning, the F/T information is crucial for effective
pruning. When cutting tool and moving vine interact, the force is experienced
by F/T sensor and is sent to the planner via NAC as feedback to execute a cut
on the vine at the instant of interaction. Since, the cutting action defines
the number of attempts made to reach the final assigned end goal. During
stationary pruning, the cutting tool mostly performs all the desired cuts in
one trial but for moving vines, several attempts are being made to finally
reach the required cut.

The various pruning tests under different scenarios are listed in Table 1. The
pruning accuracy is measured in terms of the cuts made by blades at the
assigned location with a correct position and orientation [9]. Since, the
value of vine force is averaged over its time to reach the pruning points. The
more complex pruning scenario is obviously a mixed one, in which the vines
have given random motions in all the defined directions and hence the cutting
tool after sensing force tries several times to reach the final cut.

Table 1: Different Setups of Pruning Scenarios. Sr. No. |  Motion |  Vine Force (N) |  Pruning time (sec) |  Pruning accuracy (%) |  Pruning actions (times)
---|---|---|---|---|---
1 |  Stationary |  2.31494 |  18 |  100 |  1
2 |  Y-trans |  6.12634 |  31 |  97.2 |  3
3 |  Y-rot |  9.79863 |  24 |  98.7 |  2
4 |  Z-rot |  10.02314 |  28 |  96.4 |  3
5 |  Mixed |  19.36213 |  49 |  95.1 |  5

##  5 CONCLUSIONS

This research proposed a multi-modal framework for automating the pruning of
dynamic vines. For the reproducibility of pruning skills under different
scenarios with the aim of sim2real skill transfer, the 3D models of vines were
constructed and rendered in blender. The RGB images taken with Intel Real-
Sense camera mounted on end-effector were used for spurs detection using
Faster RCNN at a threshold of 0.7. The robot arm then moved close to a
bounding box at the distance of $8.35cm$ between the cutting tool and vine and
then within the vicinity of spurs, the pruning points were located using
statistical pattern recognition algorithm based on the graph morphometry. The
cutting poses of these points were constructed to command the robot. The vines
were initially stationary and then moved in random directions. The proposed
framework was trained and tested on designed 3D models of vines in simulators
and then the pruning skills were also repeated on the real system in
laboratory environment. NAC was applied to deal with dynamic nature of the
vines. It was found that the pruning timing was compromised while dealing with
moving vines but the accuracy was still maintained within the $\pm 5\%$
tolerance limits.

The current research has been carried out in the controlled laboratory
environment with limited sim2real skill transfer. The next goal is to test the
learned skills in real vineyard on the field and also to increase the autonomy
of robot platform by using deep SLAM under distinct environmental conditions.

## References

  * [1] Adhikari, S., and Kandel, T. P. (2015). Effect of time and level of pruning on vegetative growth, flowering, yield, and quality of guava. International journal of fruit science, 15(3), 290-301.
  * [2] Osentoski, S., Jay, G., Crick, C., Pitzer, B., DuHadway, C., and Jenkins, O. C. (2011). Robots as web services: Reproducible experimentation and application development using rosjs. In 2011 IEEE International Conference on Robotics and Automation (pp. 6078-6083).
  * [3] Russell, S., and Norvig, P. (2002). Artificial intelligence: a modern approach; Pearson Education Limited.
  * [4] Matas, J., James, S., and Davison, A. J. (2018). Sim-to-real reinforcement learning for deformable object manipulation. arXiv preprint arXiv:1806.07851.
  * [5] Wang, Z., and Hirai, S. (2011). Modeling and estimation of rheological properties of food products for manufacturing simulations. Journal of food engineering, 102(2), 136-144.
  * [6] Katyara, S., Ficuciello, F., Caldwell, D., Siciliano, B., and Chen, F. (2020). Leveraging Kernelized Synergies on Shared Subspace for Precision Grasp and Dexterous Manipulation. arXiv preprint arXiv:2008.11574.
  * [7] Biagiotti, L., Melchiorri, C., and Vassura, G. (2001). Position/force control of an arm/gripper system for space manipulation. In 2001 IEEE/ASME International Conference on Advanced Intelligent Mechatronics. Proceedings (Cat. No. 01TH8556) (Vol. 2, pp. 1175-1180)
  * [8] Colgate, J. E., and Hogan, N. (1987). On the stability of a manipulator interacting with its environment. In Proc. 25th Allerton Conf. Communication Control and Computing.
  * [9] Newman, W. S. (1992). Stability and performance limits of interaction controllers. Journal of dynamic systems, measurement, and control, 114(4):563-570.
  * [10] He, L., and Schupp, J. (2018). Sensing and automation in pruning of apple trees: A review. Agronomy, 8(10), 211.
  * [11] Silwal, A., Davidson, J. R., Karkee, M., Mo, C., Zhang, Q., and Lewis, K. (2017). Design, integration, and field evaluation of a robotic apple harvester. Journal of Field Robotics, 34(6), 1140-1159.
  * [12] Botterill, T., Paulin, S., Green, R., Williams, S., Lin, J., Saxton, V., … and Corbett‐Davies, S. (2017). A robot system for pruning grape vines. Journal of Field Robotics, 34(6), 1100-1122.
  * [13] Korayem, M. H., Shafei, A. M., and Seidi, E. (2014). Symbolic derivation of governing equations for dual-arm mobile manipulators used in fruit-picking and the pruning of tall trees. Computers and electronics in agriculture, 105, 95-102.
  * [14] Peng, X. B., Andrychowicz, M., Zaremba, W., and Abbeel, P. (2018). Sim-to-real transfer of robotic control with dynamics randomization. In 2018 IEEE international conference on robotics and automation (ICRA) (pp. 1-8).
  * [15] Hoffman, E. M., Rocchi, A., Laurenzi, A., and Tsagarakis, N. G. (2017). Robot control for dummies: Insights and examples using opensot. In 2017 IEEE-RAS 17th International Conference on Humanoid Robotics (Humanoids) (pp. 736-741).
  * [16] Faster, R. C. N. N. (2015). Towards real-time object detection with region proposal networks. Advances in neural information processing systems, pages 91-99.
  * [17] Na, S., Xumin, L., and Yong, G. (2010). Research on k-means clustering algorithm: An improved k-means clustering algorithm. In 2010 Third International Symposium on intelligent information technology and security informatics (pp. 63-67).
  * [18] Mansard, N., Stasse, O., Evrard, P., and Kheddar, A. (2009). A versatile generalized inverted kinematics implementation for collaborative working humanoid robots: The stack of tasks. In 2009 International Conference on Advanced Robotics (pp. 1-6).
  * [19] Rocchi, A., Hoffman, E. M., Caldwell, D. G., and Tsagarakis, N. G. (2015). Opensot: a whole-body control library for the compliant humanoid robot coman. In 2015 IEEE International Conference on Robotics and Automation (ICRA) (pp. 6248-6253)
  * [20] Dohring, M., and Newman, W. (2003). The passivity of natural admittance control implementations. In 2003 IEEE International Conference on Robotics and Automation (Cat. No. 03CH37422) (Vol. 3, pp. 3710-3715).


11footnotetext: HSE University22footnotetext: Faculty of Computational
Mathematics and Cybernetics, Moscow State University33footnotetext: Moscow
Center for Fundamental and Applied Mathematics44footnotetext: Research
Computing Center, Moscow State University55footnotetext: Marchuk Institute of
Numerical Mathematics of Russian Academy of Sciences

# Tensor-Train Numerical Integration of Multivariate Functions with
Singularities

Lev I. Vysotsky1,2,3
<EMAIL_ADDRESS>&Alexander V. Smirnov 4,3
<EMAIL_ADDRESS>&Eugene E. Tyrtyshnikov 5,4,3
<EMAIL_ADDRESS>
###### Abstract

Numerical integration is a classical problem emerging in many fields of
science. Multivariate integration cannot be approached with classical methods
due to the exponential growth of the number of quadrature nodes. We propose a
method to overcome this problem. Tensor-train decomposition of a tensor
approximating the integrand is constructed and used to evaluate a multivariate
quadrature formula. We show how to deal with singularities in the integration
domain and conduct theoretical analysis of the integration accuracy. The
reference open-source implementation is provided.

_K_ eywords multidimensional integration $\cdot$ multivariate integration
$\cdot$ numerical integration $\cdot$ sampling point $\cdot$ tensor train
$\cdot$ TT-format

##  1 Introduction

Numerical integration is a classical problem that emerges in many fields of
science. This paper was inspired by a physical application: a numerical
approach to the evaluation of Feynman integrals based on so-called sector
decomposition (see, for example, [1]).

The classical naive approach (see e.g. [2]) to multivariate integration meets
the so-called curse of dimensionality. Let us for simplicity consider the case
of integrating over a unitary hypercube of dimension $d$. In the classical
approach to multivariate integration one fixes $n$ nodes on each of the
coordinate axes and sums the values of the integrand in each of the $n^{d}$
points with weights prescribed by the used quadrature rule (i.e. Gaussian or
Simpson’s). However whatever rule is chosen, the $n^{d}$ grows exponentially
with $d$, so the classical approaches turn out to be completely inefficient in
high-dimensional case.

There exist different approaches to the multivariate integration avoiding the
curse of dimensionality, both deterministic and nondeterministic. The
nondeterministic methods are usually called Monte-Carlo algorithms.
Historically the first method of this class is the Vegas algorithm originally
presented in [3, 4]. It uses importance sampling for variance reduction. One
might pose a questions why there is a place for randomness in a task when a
deterministic result is required. The idea is that the algorithms randomly
disseminating sampling points have almost a hundred percent chance of hitting
the “special” regions of integration domain. On the contrary, a deterministic
improperly chosen method might skip such regions thus leaning to an incorrect
result. Hence while using random point sampling those algorithms are supposed
to end in almost the same result each time with the difference being
negligible and smaller than the error estimate.

There are other known nondeterministic algorithms such as Suave combining the
advantages of two popular methods: importance sampling by Vegas and subregion
sampling in a manner similar to Miser [5]. By dividing into subregions, Suave
manages to a certain extent to get around Vegas’ difficulty to adapt its
weight function to structures not aligned with the coordinate axes. However,
the disadvantage of this method lies it the obligation of the user to provide
some extra information of the function structure which is not possible in a
general case.

There are also deterministic methods known such as Cuhre using a cubature rule
for subregion estimation in a globally adaptive subdivision scheme [7].
However, such approaches still suffer from the curse of dimensionality and
might compete with the Monte-Carlo class only at lower dimensions such as 2–4,
while in the physical problem that inspired this work one might expect up to
10 dimensions or even more.

The development of the Monte-Carlo approach led to emergence of the so-called
Quasi Monte Carlo (QMC) integration methods [6, 8]. The review of them can be
found in [9], where one of them was applied for the evaluation of Feynman
integrals. For a description of a modern implementation of the QMC algorithm
see [10].

A completely new approach to multivariate integration based on the tensor-
train decomposition was originally suggested in [11]. The idea is to use low-
parametric approximation (tensor-train, or TT, decomposition) to the
multivariate function to overcome the curse of dimensionality of the naive
multivariate integration. The proposed decomposition can be constructed using
the function values computed in moderate number of points. It also allows to
efficiently calculate the multidimensional weighted sum implied by quadrature
rule. However, in [11] the authors provided only a “proof-of-concept” of TT
integration and applied it to some model functions. Recently, in [12] the
method was used to compute high-dimensional (up to 1024 dimensions) integrals
and was shown to outperform the stochastic methods. Moreover, the method seems
applicable to very high-dimensional integrals (up to hundreds of thousands of
dimensions) like those found e.g. in condensed matter physics, as the
computational complexity grows only linearly with the number of dimensions.

In this paper we provide a complete description of the tensor-train-based
multivariate integration of functions with singularities. A configurable and
practically applicable algorithm is presented and its reference open-source
implementation is provided. We conduct the theoretical error analysis and show
how to deal with some types of singularities in the integration domain. The
algorithm was applied to the computation of a Feynman integral and shown to
perform better than the Monte-Carlo-based method.

The rest of the paper is organized as follows. In Section 2 we introduce
concepts and facts necessary for the formulation of the method. Section 3.1
contains the bird’s eye view of the integration algorithm. In Sections 3.2 and
3.3 the details of TT cross approximation approach are provided along with
some supporting arguments. In Section 3.4 we consider the problem of choosing
the approximation tolerance given the limit on the number of function
evaluations. Our approach to accurate integration of functions with
singularities is described in Section 3.5. The theoretical analysis of the
integration accuracy is provided in Section 3.6. Finally, the discussion of
the implementation details and usage of the library can be found in Section 4.
There we also provide a comparison of the TT and Monte-Carlo-based methods
applied to a computation of a Feynman integral.

##  2 Definitions and preliminary facts

Let $\mathbb{C}^{n_{1}\times\cdots\times n_{d}}$ denote the set of
$d$-dimensional tensors (i.e. arrays) over field $\mathbb{C}$ with mode sizes
$n_{1},\dots,n_{d}$. For a given tensor
$A\in\mathbb{C}^{n_{1}\times\cdots\times n_{d}}$ denote $A(i_{1},\dots,i_{d})$
its element with indices $i_{1},\dots,i_{d}$. We will use zero-based
indexation, i.e. $0\leq i_{k}\leq n_{k}-1$.

If in indexation expression $A(i_{1},\dots,i_{d})$ we replace variables
$i_{k_{1}},\dots,i_{k_{s}}$ with colons we obtain an $s$-dimensional tensor
$B\in\mathbb{C}^{n_{k_{1}},\dots,n_{k_{s}}}$ with elements
$B(i_{k_{1}},\dots,i_{k_{s}})=A(i_{1},\dots,i_{d}).$ The tensor $B$ is thus
determined by “non-colon” variables. As an example, consider a 3-dimensional
tensor $A\in\mathbb{R}^{2\times 3\times 4}$ with elements $A(i,j,k)=i+j+k$.
Then $B:=A(:,:,2)$ is a $2\times 3$ matrix with elements $B(i,j)=i+j+2$, and
$v:=A(1,1,:)$ is a vector of length $4$ with elements $v(k)=k+2$. The notation
is adopted from Matlab and Python programming languages. Other useful
adoptions are the $\mathsf{vec}$ and $\mathsf{reshape}$ functions. The former
one makes a “long” vector:

$\mathsf{vec}(A)\in\mathbb{C}^{n_{1}\cdots
n_{d}},~{}\mathsf{vec}(A)(i_{1}+i_{2}n_{1}+\dots+i_{d}n_{1}\dots
n_{d-1})=A(i_{1},\dots,i_{d}),$

The $\mathsf{reshape}$ function takes a tensor
$A\in\mathbb{C}^{n_{1}\times\dots\times n_{d}}$ and new shape
$m_{1}\times\dots\times m_{D}$ and returns tensor
$B\in\mathbb{C}^{m_{1}\times\dots\times m_{D}}$ such that
$\mathsf{vec}(A)=\mathsf{vec}(B)$ (naturally, $n_{1}\cdots n_{d}$ must be
equal to $m_{1}\cdots m_{D}$). For each $k=1,\dots,d-1$ the $k$-th unfolding
matrix $A_{k}$ of a given tensor $A$ is defined as

$A_{k}=\mathsf{reshape}(A,(n_{1}\cdots n_{k})\times(n_{k+1}\cdots n_{d})).$

We say that a tensor $A\in\mathbb{C}^{n_{1}\times\cdots\times n_{d}}$ admits
_tensor-train (TT) decomposition_ [13] if

$A(i_{1},\dots,i_{d})=\sum_{\alpha_{0}=0}^{r_{0}-1}\dots\sum_{\alpha_{d}=0}^{r_{d}-1}G_{1}(\alpha_{0},i_{1},\alpha_{1})G_{2}(\alpha_{1},i_{2},\alpha_{2})\dots
G_{d}(\alpha_{d-1},i_{d},\alpha_{d})$ (1)

for some natural numbers $r_{1},\dots,r_{d-1}$ and 3-dimensional tensors
$G_{k}\in\mathbb{C}^{r_{k-1}\times r_{k}\times r_{k+1}}$. Note that we force
$r_{0}$ and $r_{d}$ to be equal to 1, thus in equation (1) variables
$\alpha_{0}$ and $\alpha_{d}$ are “dummy” ones and tensors $G_{1}$ and $G_{d}$
are actually 2-dimensional. Tensors $G_{k}$ are called the _cores_ and numbers
$r_{k}$ the _compression ranks_ (or TT-ranks) of decomposition (1). An
equivalent way to define tensor-train decomposition is

$A(i_{1},\dots,i_{d})=G_{1}(:,i_{1},:)\dots G_{d}(:,i_{d},:).$ (2)

Note that each multiplier $G_{k}(:,i_{k},:)$ above is a matrix (one index is
fixed while the two other are “free”) of size $r_{k-1}\times r_{k}$. The
product is thus correctly defined, moreover, as $r_{0}=r_{d}=1$ it has size
$1\times 1$, so it is indeed a number from $\mathbb{C}$. It can be easily
checked that definitions (1) and (2) are equivalent. The latter is best suited
for computations as it exploits fast BLAS level-3 operation of matrix
multiplication. The fact that $A$ admits TT-decomposition with cores
$G_{1},\dots,G_{d}$ is denoted $A=\mathsf{TT}(G_{1},\dots,G_{d})$.

Obviously, the storage required for a given TT-representation of a tensor is
proportional to $\sum_{k}n_{k}r_{k-1}r_{k}$ or if all mode sizes are bounded
from above by $n$ and compression ranks by $r$, the storage can be estimated
as $O(dnr^{2})$.

For a given tensor $A\in\mathbb{C}^{n_{1}\times\cdots\times n_{d}}$ we can (at
least theoretically) find the smallest compression ranks among all the
possible TT-decompositions of it. These compression ranks are equal to
corresponding ranks of unfolding matrices: $r_{k}=\mathrm{rank}\,A_{k}$. The
existence of a TT-decomposition with compression ranks equal to
$\mathrm{rank}\,A_{k}$ is proven in [13]. The non-existence of TT-
decompositions with smaller compression ranks follows from the fact that any
TT-decompositions with TT-ranks $r_{k}$ yields a skeleton decomposition of
unfolding matrix $A_{k}$ with $r_{k}$ summands.

In practice, however, we are interested in low rank TT-approximations, i.e.
tensors $A^{\prime}$ close to $A$ in some norm that admit TT-decomposition
with low ranks. In [11] it was proven that if unfolding matrices are
approximated by rank-$r_{k}$ matrices $B^{(k)}$:
$\|A_{k}-B^{(k)}\|_{F}\leq\varepsilon_{k},$ $k=1,\dots,d-1$, then there exists
a tensor $A^{\prime}$ with compression ranks $r_{k}$ satisfying

$\|A^{\prime}-A\|_{F}^{2}\leq\sum_{k=1}^{d-1}\varepsilon_{k}^{2}$ (3)

(here $\|A\|_{F}$ denotes Frobenius norm of a tensor, i.e. square root of
$\sum_{i_{1},\dots,i_{d}}|A(i_{1},\dots,i_{d})|^{2}$).

If tensors $A$ and $B$ are represented in TT-format (i.e. some TT-
decomposition is known for each of them), several basic operations may be
performed directly using such representation [13]. Addition, Hadamard
(element-wise) product, dot product and multidimensional contraction are the
examples of such operations. The dot product of tensors $\langle A,B\rangle$
is defined as $\langle A,B\rangle=\sum_{i_{i},\dots,i_{d}}(\overline{A}\circ
B)(i_{1},\dots,i_{d}),$ where $\overline{A}$ is tensor whose elements are
complex conjugates of elements of $A$:
$\overline{A}(i_{1},\dots,i_{d})=\overline{A(i_{1},\dots,i_{d})}$ and $A\circ
B$ denotes element-wise product of $A$ and $B$. Multidimensional contraction
of tensors $A$ with columns $u_{k}\in\mathbb{C}^{n_{k}}$ is defined as
$W=\langle\overline{A},U\rangle$, where $U$ is a tensor with compression ranks
equal to 1 and cores $G_{k}(0,:,0)=u_{k}$.

##  3 The integration algorithm

###  3.1 High-level description

Suppose we are given a function $f:\mathbb{R}^{d}\to\mathbb{C}$ and need to
compute (approximately) the integral

$I(f)=\int_{0}^{1}\dots\int_{0}^{1}f(x_{1},\dots,x_{d})dx_{1}\cdots dx_{d}.$
(4)

We assume that some nodes $\\{x_{k}^{(i_{k})}\\}_{i_{k}=0}^{n_{k}-1}$,
$x_{k}^{(i_{k})}\in\mathbb{R}$, and weights
$\\{w_{k}^{(i_{k})}\\}_{i_{k}=0}^{n_{k}-1}$, $w_{k}^{(i_{k})}\in\mathbb{R}$,
are fixed for each axis $x_{k}$. Consider the $d$-dimensional grid $X$:

$X:=\\{x_{1}^{(i_{1})}\\}_{i_{1}=0}^{n_{1}-1}\times\dots\times\\{x_{1}^{(i_{d})}\\}_{i_{d}=0}^{n_{d}-1}$

and weight tensor $W\in\mathbb{C}^{n_{1}\times\cdots\times n_{d}}$ with
elements $W(i_{1},\dots,i_{d})=w_{1}^{(i_{1})}\cdots w_{d}^{(i_{d})}.$ We will
approximate integral (4) with the sum

$S_{X,W}(f):=\sum_{i_{1}=0}^{n_{1}-1}\dots\sum_{i_{d}=0}^{n_{d}-1}f\left(x_{1}^{(i_{1})},\dots,x_{d}^{(i_{d})}\right)W(i_{1},\dots,i_{d}).$
(5)

If we denote $A_{X,f}$ the tensor from $\mathbb{C}^{n_{1}\times\cdots\times
n_{d}}$ with elements $f\big{(}x_{1}^{(i_{1})},\dots,x_{d}^{(i_{d})}\big{)}$,
the introduced sum can be written as
$S_{X,W}(f)=\langle\overline{A_{X,f}},W\rangle$. But even if weights
$w_{k}^{(i_{k})}$ are chosen wisely and the sum $S_{X,W}(f)$ accurately
approximates $I(f)$, this approach is not free from the “curse of
dimensionality” as direct computation of $S_{X,W}(f)$ requires $n_{1}\cdots
n_{d}$ evaluations of $f(x_{1},\dots,x_{d})$. Let us assume, however, that
$A_{X,f}$ admits an approximation $A^{\prime}$ with moderate TT-ranks $r_{k}$.
Then the value $S^{\prime}=\langle\overline{A^{\prime}},W\rangle$ can be
computed as multidimensional contraction (as tensor $W$ has compression ranks
equal to 1) with complexity $O\left(\sum_{k}n_{k}r_{k}r_{k-1}\right)$ which is
polynomial in mode sizes $n_{k}$ and compression ranks $r_{k}$.

The main difficulty consists in obtaining a good approximation of $A_{X,f}$
without computing all (or even a significant ratio of) its entries. For many
functions arising from real-world applications this is an achievable target.

Algorithm 1 presents the pseudocode for constructing a TT-approximation of a
given black-box (i.e. with specified procedure of computing an element with
given indices) tensor. It is described in detail in Section 3.3. The workhorse
of the method are the subroutines TTCrossLeftToRight (see Algorithm 2) and
TTCrossRightToLeft (pseudocode is not presented due to its extreme
similarity). After each pass we obtain a (hopefully) improved approximation to
$A$. During one pass for each dimension $k$ a carefully selected submatrix of
the unfolding matrix $A_{k}$ is approximated by a cross consisting of $r_{k}$
rows and $r_{k}$ columns. For this we use the approach which pseudocode is
presented in Algorithm 3 and which is described in detail in Section 3.2). The
obtained cross is used to construct the $k$-th TT core of the approximation.
The lines 1 and 1 of Algorithm 1 are responsible for choosing the
approximation accuracy according to the user-specified limit on the number of
tensor element evaluations (for details refer to Section 3.4).

Input: (black-box) tensor $A\in\mathbb{C}^{n_{1}\times\dots\times n_{d}}$,
evaluation count limit $N_{\mathrm{lim}}$

Output: TT cores $G_{1},\dots,G_{d}$ of approximation $A^{\prime}$ of $A$

1 Randomly choose sets of column multi-indices
$\widehat{\mathcal{J}}_{1},\dots,\widehat{\mathcal{J}}_{d-1}$;

$G_{1},\dots,G_{d},\mathcal{I}_{1},\dots,\mathcal{I}_{d-1}\leftarrow\mathtt{TTCrossLeftToRight}(A,\widehat{\mathcal{J}}_{1},\dots,\widehat{\mathcal{J}}_{d-1},\varepsilon_{\mathrm{test}})$
// measure $N_{\mathrm{test}}$

2 $\varepsilon_{\mathrm{eff}}\leftarrow$ effective tolerance (see Section 3.4)
;

3 repeat

4
$\widehat{\mathcal{I}}_{1},\dots,\widehat{\mathcal{I}}_{d-1}\leftarrow\mathcal{I}_{1},\dots,\mathcal{I}_{d-1}$
extended with random multi-indices;

5
$H_{1},\dots,H_{d},\mathcal{J}_{1},\dots,\mathcal{J}_{d-1}\leftarrow\mathtt{TTCrossRightToLeft}(A_{X,f},\widehat{\mathcal{I}}_{1},\dots,\widehat{\mathcal{I}}_{d-1},\varepsilon_{\mathrm{eff}})$;

6
$\widehat{\mathcal{J}}_{1},\dots,\widehat{\mathcal{J}}_{d-1}\leftarrow\mathcal{J}_{1},\dots,\mathcal{J}_{d-1}$
extended with random multi-indices;

7
$G_{1},\dots,G_{d},\mathcal{I}_{1},\dots,\mathcal{I}_{d-1}\leftarrow\mathtt{TTCrossLeftToRight}(A_{X,f},\widehat{\mathcal{J}}_{1},\dots,\widehat{\mathcal{J}}_{d-1},\varepsilon_{\mathrm{eff}})$;

8

9until  _
$\Big{\|}\mathsf{TT}(H_{1},\dots,H_{d})-\mathsf{TT}(G_{1},\dots,G_{d})\Big{\|}_{F}<\varepsilon\Big{\|}\mathsf{TT}(G_{1},\dots,G_{d})\Big{\|}_{F}$
_;

10return $G_{1},\dots,G_{d}$;

Algorithm 1 TTCross constructs a TT approximation of a tensor (Section 3.3)

Input: (black-box) tensor $A\in\mathbb{C}^{n_{1}\times\dots\times n_{d}}$,
sets of column multi-indices
$\widehat{\mathcal{J}}_{1},\dots,\widehat{\mathcal{J}}_{d-1}$, approximation
tolerance $\varepsilon$

Output: TT-approximation of $A$ with cores $G_{1},\dots,G_{d-1}$, sets of row
multi-indices $\mathcal{I}_{1},\dots,\mathcal{I}_{d-1}$

$\mathcal{I}_{0}\leftarrow\\{()\\}$ // empty multi-index

1 $r_{0}\leftarrow 1$;

2 for  _$k\leftarrow 1,\dots,d-1$_ do

3
$\widehat{\mathcal{I}}_{k}\leftarrow\big{\\{}(i_{1},\dots,i_{k})~{}:~{}(i_{1},\dots,i_{k-1})\in\mathcal{I}_{k-1},i_{k}\in\\{0,\dots,n_{k}-1\\}\big{\\}}$;

4
$\mathcal{I}_{k},\mathcal{J}_{k}\leftarrow\mathtt{MatrixCrossApproximate}\big{(}A(\widehat{I}_{k},\widehat{\mathcal{J}}_{k}),\varepsilon/\sqrt{d-1}\big{)}$;

5 $C_{k}\leftarrow
A(\widehat{\mathcal{I}}_{k},\mathcal{J}_{k})\big{(}A(\mathcal{I}_{k},\mathcal{J}_{k}\big{)}^{+}_{\tau}$;

6 $r_{k}\leftarrow|\mathcal{I}_{k}|$;

7 Store $G_{k}\leftarrow\mathsf{reshape}(C_{k},r_{k-1}\times n_{k}\times
r_{k})$;

8

9$G_{d}\leftarrow A(\mathcal{I}_{d-1},:)$;

10 return $G_{1},\dots,G_{d}$, $\mathcal{I}_{1},\dots,\mathcal{I}_{d-1}$;

Algorithm 2 TTCrossLeftToRight performs left-to-right pass of TT-cross
approximation (Section 3.3)

Input: (black-box) matrix $M\in\mathbb{C}^{m\times n}$, approximation
tolerance $\widehat{\varepsilon}$

Output: Row and column indices $\mathcal{I},\mathcal{J}$ such that
$M(:,\mathcal{J})M(\mathcal{I},\mathcal{J})^{-1}M(\mathcal{I},:)$ approximates
$M$

1 $M^{\prime},j_{1},r\leftarrow M,0,0$;

2 repeat

3 $r\leftarrow r+1$;

4 $i_{r}\leftarrow\arg\max_{i}|M^{\prime}(i,j_{r})|$;

5 $j_{r+1}\leftarrow\arg\max_{j}|M^{\prime}(i_{r},j)|$;

6 $\mu\leftarrow|M^{\prime}(i_{r},j_{r})|$;

7 $M^{\prime}\leftarrow
M^{\prime}-M^{\prime}(:,j_{r})(M^{\prime}(i_{r},j_{r}))^{-1}M^{\prime}(i_{r},:)$;

8

9until  _ $\mu\sqrt{(m-r)(n-r)}<\widehat{\varepsilon}\|M-M^{\prime}\|_{F}$ _;

10return $\\{i_{1},\dots,i_{r}\\},\\{j_{1},\dots,j_{r}\\}$;

Algorithm 3 MatrixCrossApproximate computes a cross approximation of a matrix
(Section 3.2)

###  3.2 Cross approximation of matrices

In this section we describe methods to approximate a matrix when we are
allowed to access only a small fraction of its elements. There is a well-known
result that any matrix $M\in\mathbb{C}^{m\times n}$ of rank $r$ satisfies

$M=M(:,\mathcal{J})M(\mathcal{I},\mathcal{J})^{-1}M(\mathcal{I},:)$

where $\mathcal{I}$ and $\mathcal{J}$ are any size-$r$ row and column index
sets such that submatrix at their intersection is nonsingular. A more subtle
result is that a matrix with approximate rank $r$ can be sometimes
approximated with matrix $M(:,\mathcal{J})GM(\mathcal{I},:)$ for some $r\times
r$ matrix $G$. The most intriguing problem here is the choice of $r$ rows and
columns. It turns out that it makes sense to choose the submatrix with maximum
volume (modulus of the determinant) among all the $r\times r$ submatrices: in
[18] it was proven that if $M(\mathcal{I},\mathcal{J})$ is such a submatrix
and moreover is nonsingular, then

$\left\|M-M(:,\mathcal{J})M(\mathcal{I},\mathcal{J})^{-1}M(\mathcal{I},:)\right\|_{C}\leq(r+1)\sigma_{r+1}(M).$
(6)

Here $\|M\|_{C}$ denotes Chebyshev norm ($\|M\|_{C}=\max_{i,j}|M(i,j)|$) and
$\sigma_{r+1}(M)$ the $(r+1)$-th singular value of $M$. The inversion of the
matrix $M(\mathcal{I},\mathcal{J})$ can, however, be numerically unstable if
the latter is ill-conditioned. Actually, the pseudoinverse can be used instead
of the inverse: in [19] it was proven that if again the submatrix
$M(\mathcal{I},\mathcal{J})$ has maximum volume and is nonsingular then

$\left\|M-M(:,\mathcal{J})\big{(}M(\mathcal{I},\mathcal{J})\big{)}^{+}_{\tau}M(\mathcal{I},:)\right\|_{2}\leq
c(\sigma_{r+1}(M)+\tau)\sqrt{1+r(\max\\{m,n\\}-r)},$

where $c$ and $\tau$ are positive constants and $G^{+}_{\tau}$ denotes
$\tau$-pseudoinverse of $G$ (if a matrix $B$ has singular decomposition
$B=U\Sigma V^{*}$, and $\Sigma_{\tau}$ coincides with $\Sigma$ except the
elements smaller than $\tau$ are set to zero, then $\tau$-pseudoinverse of $B$
is the pseudoinverse of $U\Sigma_{\tau}V^{*}$).

The important practical question is how to find the submatrix of maximal
volume. Alas, this problem is NP-hard [20], so we usually are looking for a
submatrix with “sufficiently” large volume. One approach was proposed in [21]
and it is used in the current implementation. The idea is to construct
approximations to $M$ with increasing ranks. We start with a random column
index $j_{1}$, choose the row with index $i_{1}$ satisfying
$i_{1}=\arg\max_{i}|M(i,j_{1})|$ and use
$M_{1}=M(:,j_{1})(M(i_{1},j_{1}))^{-1}M(i_{1},:)$ as the first approximation.
These steps are repeated for the residual matrix $M-M_{1}$, for the next
residual etc., until the latest residual $M-M_{1}-\dots-M_{r}$ is small enough
or the limit for $r$ is reached. The indices
$\mathcal{I}=\\{i_{1},\dots,i_{r}\\}$ and
$\mathcal{J}=\\{j_{1},\dots,j_{r}\\}$ often yield a good submatrix.

The Algorithm 3 presents a more rigorous description of the method. There are
two nontrivial moments in it. First, the line 3 looks like we need to compute
and store the whole $m\times n$ matrix. In fact, $M^{\prime}$ is a black-box
matrix and its elements are computed basing on the corresponding elements of
$M$ and rank-1 matrices $M_{k}$ (not present explicitly in pseudocode).
Second, the stopping criterion in line 3 may also look weird. Actually, the
left part of it is an estimate of $\|M^{\prime}\|_{F}$ as there are not more
than $(m-r)(n-r)$ nonzero elements in this matrix and its largest in modulus
element does not (hopefully) strongly exceed the value of $\mu$. The value
$\|M-M^{\prime}\|_{F}$ in the right part of the comparison is our current
rank-$r$ approximation of $M$ and can be computed efficiently. Indeed, if we
write down the rank-1 matrices $M_{k}$ as $M_{k}=u_{k}v_{k}^{*}$ and denote
$U:=[u_{1},\dots,u_{r}]$ and $V=[v_{1},\dots,v_{r}]$, then

$\|M-M^{\prime}\|_{F}^{2}=\|M_{1}+\dots+M_{r}\|_{F}^{2}=\|UV^{*}\|_{F}^{2}=\mathrm{Tr}(VU^{*}UV^{*})=\mathrm{Tr}(U^{*}UV^{*}V).$

Matrices $U^{*}U$ and $V^{*}V$ have size $r\times r$ and can be computed in
time $O((m+n)r^{2})$.

###  3.3 TT cross approximation

In this section we give the details of algorithms 1 and 2 for computing the TT
approximation of a black-box tensor. Suppose we are given $d-1$ sets of multi-
indices

$\displaystyle\widehat{\mathcal{J}}_{1}=\left\\{\left(i_{2}^{(1)},\dots,i_{d}^{(1)}\right),\dots,\left(i_{2}^{(R_{1})},\dots,i_{d}^{(R_{1})}\right)\right\\},\cdots,\widehat{\mathcal{J}}_{d-1}=\left\\{i_{d}^{(1)},\dots,i_{d}^{(R_{d-1})}\right\\}.$

Let us consider the $n_{1}\times R_{1}$ matrix
$A(:,\widehat{\mathcal{J}}_{1})$ and choose its submatrix
$\widehat{A}_{1}=A(\mathcal{I}_{1},\mathcal{J}_{1})\in\mathbb{C}^{r_{1}\times
r_{1}}$ of large volume as described in Section 3.2. Also store the
corresponding columns as $C_{1}:=A(:,\mathcal{J}_{1})$. Next, let us take the
subtensor $A(\mathcal{I}_{1},:,\dots,:)$ of size $r_{1}\times
n_{2}\times\dots\times n_{d}$ and consider the new tensor

$B_{1}:=\mathsf{reshape}\left(A(\mathcal{I}_{1},:,\dots,:),r_{1}n_{2}\times
n_{3}\times\dots\times n_{d}\right).$

We continue in the same vein with the tensor $B_{1}$ and index sets
$\widehat{\mathcal{J}}_{2},\dots,\widehat{\mathcal{J}}_{d-1}$. Eventually for
each $k=1,\dots,d-1$ the matrices $\widehat{A}_{k}\in\mathbb{C}^{r_{k}\times
r_{k}}$ and $C_{k}\in\mathbb{C}^{r_{k-1}n_{k}\times r_{k}}$ will be
constructed. On the last step we work with 2-dimensional tensor $B_{d-2}$ of
size $r_{d-2}n_{d-1}\times n_{d}$. Now along with matrices $C_{d-1}$ and
$\widehat{A}_{d-1}$ we store the row matrix
$C_{d}:=B_{d-2}(\mathcal{I}_{d-1},:)$. In the paper [11] it is proven that if
the numbers $r_{1},\dots,r_{d-1}$ happen to coincide with the TT-ranks of $A$
and all of $\widehat{A}_{k}$ are nonsingular, then $A$ admits TT-decomposition
with cores

$\displaystyle G_{1}$ $\displaystyle=C_{1}\widehat{A}_{1}^{-1},$
$\displaystyle G_{k}$
$\displaystyle=\mathsf{reshape}(C_{k}\widehat{A}_{k}^{-1},r_{k-1}\times
n_{k}\times r_{k}),~{}~{}k=2,\dots,d-1,$ $\displaystyle G_{d}$
$\displaystyle=C_{d}.$

However, we are usually interested in low-rank approximations so our $r_{k}$
are likely to be smaller than actual TT-ranks of $A$. Moreover, as discussed
in Section 3.2, it is more robust to use $\tau$-pseudoinverse of
$\widehat{A}_{k}$ instead of $\widehat{A}_{k}^{-1}$. Thus we construct only an
approximate TT-decomposition and beside that its quality depends on the choice
of index sets $\widehat{\mathcal{J}}_{k}$. So to find a good approximation we
need to repeat the whole procedure several times. We start with arbitrary sets
$\widehat{\mathcal{J}}_{k}$ and after the execution of the above algorithm we
obtain multi-index sets of rows, namely
$\mathcal{I}_{1},\dots,\mathcal{I}_{d-1}$, $|\mathcal{I}_{k}|=r_{k}$. These
sets are extended with random multi-indices to get the sets
$\widehat{\mathcal{I}}_{1},\dots,\widehat{\mathcal{I}}_{d-1}$,
$|\widehat{\mathcal{I}}_{k}|=R_{k}$. Now we run the same algorithm “from right
to left”, i.e. we start with $A(\widehat{\mathcal{I}}_{d-1},:)$ and choose a
submatrix of large volume in it. The we proceed for $k=d-2,\dots,1$, finally
obtaining new approximation to $A$.

Assume that after the $j$-th pass of the above algorithm we obtain the TT-
decomposition of the approximation $A^{(j)}$. How can we estimate its
closeness to $A$ or at least the convergence of the process? In fact, we can
efficiently compute $\|A^{(j)}-A^{(j-1)}\|_{F}$ to compare it with prescribed
tolerance:

$\left\|A^{(j)}-A^{(j-1)}\right\|_{F}^{2}=\left\langle
A^{(j)},A^{(j)}\right\rangle+\left\langle
A^{(j-1)},A^{(j-1)}\right\rangle-2\mathrm{Re}\left\langle
A^{(j)},A^{(j-1)}\right\rangle.$

Each dot product can be computed efficiently using the TT-decompositions of
$A^{(j-1)}$ and $A^{(j)}$ as described in [13].

Note that for the sake of more straightforward implementation in the
pseudocode of Algorithm 2 we use row multi-indices $\mathcal{I}_{k}$ instead
of tensors $B_{k}$. Also notice that in line 2 we call subroutine
MatrixCrossApproximate with $\widehat{\varepsilon}=\varepsilon/\sqrt{d-1}$.
The idea is that if we use $\widehat{\varepsilon}$-approximations for the
unfolding matrices, the inequality (3) lets us hope that the overall error
will be bounded by $\widehat{\varepsilon}\sqrt{d-1}=\varepsilon$.

###  3.4 Choosing approximation tolerance

This section describes our solution to the problem of choosing the
approximation tolerance in Algorithm 1. On the one hand, this choice can be
made by the user (a rough estimate of the desired integration error can be
used; for more formal analysis refer to Section 3.6). On the other hand,
sometimes the user is limited by computational resources available for the
evaluation of an integral. In such a situation the Monte-Carlo methods are
more flexible as one can straightforwardly provide the limit on the number of
function evaluations. In the TT integration approach, however, all the
function evaluations are performed inside the MatrixCrossApproximate
subroutine. If at some moment we find out that the limit is reached, we must
abort immediately thus not obtaining a high-quality approximation of the
current matrix. Moreover, if this happens (in one of the executions of line 2
of Algorithm 2), we cannot proceed with the current TT-cross pass and are
doomed to use the previous approximation. And what if it was the first pass?

As you can see, the “hard limit” approach does not work for our method.
Instead, we propose another strategy. It is based on the observation that for
many real-life tensors the singular values of the unfolding matrices decay
exponentially, i.e.

$\sigma_{\ell}(A_{k})\approx\alpha_{k}e^{-\beta_{k}\ell},~{}\ell=1,\dots,L_{k},$

where $\alpha_{k},\beta_{k}$ are some positive constants and $L_{k}$ is the
smallest size of $A_{k}$. If the rank-$r_{k}$ approximation to $A_{k}$ found
in MatrixCrossApproximate is close enough to the optimal rank-$r_{k}$
approximation, we can write

$\widehat{\varepsilon}^{2}\approx\sum_{\ell=r_{k}+1}^{L_{k}}(\sigma_{\ell}(A_{k}))^{2}\approx\frac{\alpha_{k}^{2}}{2\beta_{k}}e^{-2\beta_{k}r_{k}}\Rightarrow
r_{k}\approx-\beta_{k}\ln\widehat{\varepsilon}+\gamma_{k}.$

As we have mentioned, all tensor element evaluations happen in the
MatrixCrossApproximate subroutine. The matrix passed to the subroutine on the
$k$-th step has size $r_{k-1}n_{k}\times R_{k}$. All the elements of a cross
of width $r_{k}$ (i.e. exactly $r_{k-1}n_{k}r_{k}+(R_{k}-r_{k})r_{k}$ ones)
are evaluated. If we assume for simplicity that $R_{k}-r_{k}\ll r_{k-1}n_{k}$,
we can write for the total number of element evaluations $N$:

$N\approx\sum_{k=1}^{d-1}n_{k}r_{k-1}r_{k}\approx(\ln\widehat{\varepsilon})^{2}\Big{(}\sum_{k=1}^{d-1}n_{k}\beta_{k-1}\beta_{k}\Big{)}+\sum_{k=1}^{d-1}\gamma_{k}.$

So if we ignore the sum of $\gamma_{k}$, we can assume that the number of
evaluations is proportional to the squared logarithm of the tolerance.

Now, our heuristic works as follows (lines 1 and 1 of Algorithm 1). First, we
run the Algorithm 2 with some “test” tolerance $\varepsilon_{\mathrm{test}}$
(e.g. $0.01$). Let us assume it has required $N_{\mathrm{test}}$ evaluations
of tensor elements. If the “soft limit” specified by the user is
$N_{\mathrm{lim}}$ we will run the following passes with “effective tolerance”
$\varepsilon_{\mathrm{eff}}$ which is computed basing on the formula

$\frac{N_{\mathrm{lim}}-N_{\mathrm{test}}}{N_{\mathrm{test}}}=\frac{(\ln\varepsilon_{\mathrm{eff}})^{2}}{(\ln\varepsilon_{\mathrm{test}})^{2}}.$

###  3.5 Functions with singularities

It is well known that the numerical integration of functions with
singularities (i.e. unbounded in some point of integration domain) give large
errors. There are several approaches to solve this problem (see e.g. [24]).
One of the most powerful is the variable substitution. Let us substitute $x$
with some monotonic function $g(t)$ (here $g^{-1}(x)$ denotes the inverse
function):

$\int_{0}^{1}f(x)dx=\int_{g^{-1}(0)}^{g^{-1}(1)}f(g(t))g^{\prime}(t)dt.$

After the transformation any quadrature rule (e.g. Gaussian) can be applied.
Note, however, that there is no need to apply it to the function
$f(g(t))g^{\prime}(t)$. Instead, the quadrature rule may be transformed
accordingly, i.e. instead of the nodes $t_{i}$ and weights $w_{i}$ we can use
the nodes $g(t_{i})$ and weights $w_{i}g^{\prime}(t_{i})$ and apply the new
rule directly to the function $f$. The advantage of this approach is most
clearly seen when it used inside the TT integration and saves quite a lot of
extra computations.

The appropriate function $g(t)$ is chosen basing on the location and type of
the singularity. As an example, consider logarithmic singularity at zero:
$f(x)=O(\ln x)$ as $x\to 0$ (i.e. $f(x)$ approaches negative infinity
asymptotically not faster than $\ln x$). In this case the power transformation
$x=t^{p}$, $p>1$ can be used. Indeed, if we apply the transformed integral is

$\int_{0}^{1}f(t^{p})pt^{p-1}dt=p\int_{0}^{1}O(\ln(t)t^{p-1})dt$

and now the integrand is regular (i.e. does not have singularities) as for
$\alpha>0$ the function $\ln(t)t^{\alpha}$ tends to $0$ as $t\to 0$. Other
useful substitutions include $\tanh$-$\sinh$ [25]
($x=\tanh(\frac{1}{2}\pi\sinh t)$) and erf [24]
($x=\frac{1}{2}(1-\mathrm{erf}(t))$) ones.

The substitution approach can be easily adopted for TT integration technique
if there is only one singularity in the multivariate integration domain
$[0,1]^{d}$. For simplicity let us assume that it is located at zero. Then for
each $k=1,\dots,d$ we choose a quadrature rule with nodes $x_{k}^{(i_{k})}$
and weights $w_{k}^{(i_{k})}$ such that under the corresponding transforms
$x_{k}=g_{k}(t_{k})$ the new integrand

$f(x_{1},\dots,x_{k-1},g_{k}(t_{k}),x_{k+1},\dots,x_{d})g^{\prime}_{k}(t_{k})$

viewed as a function of $t_{k}$ is regular for any
$x_{1},\dots,x_{k-1},x_{k},\dots,x_{d}$.

###  3.6 Error analysis

The integration algorithm introduces two kinds of errors: the approximation
error (tensor $A^{\prime}$ with hopefully moderate TT-ranks differs from
$A_{X,f}$) and the numerical integration error ($S_{X,W}(f)$ differs from
$I(f)$). In this section we estimate these two type of errors separately and
eventually come to the bound of $|I(f)-S^{\prime}|$.

For $k=1,\dots,d-1$ we will denote by $I_{k}(x_{k+1},\dots,x_{d})$ the
function

$I_{k}(x_{k+1},\dots,x_{d})=\int_{0}^{1}\cdots\int_{0}^{1}f(x_{1},\dots,x_{d})dx_{1}\cdots
dx_{k}.$

Naturally, we define $I_{0}(x_{1},\dots,x_{d})\equiv f(x_{1},\dots,x_{d})$ and
$I_{d}\equiv I(f)$.

######  Lemma 1.

Let us assume that for all $k=1,\dots,d$ two properties hold:

  1. 1.

the quadrature rule with nodes $x_{k}^{(i_{k})}$ and weights $w_{k}^{(i_{k})}$
has nonnegative weights and is exact for constant functions, i.e.

$w_{k}^{(i_{k})}\geq 0,~{}~{}1=\int_{0}^{1}1\cdot
dx=\sum_{i_{k}=0}^{n_{k}-1}w_{k}^{(i_{k})};$ (7)

  2. 2.

for all $x_{k+1},\dots,x_{d}\in[0,1]$ we have

$\left|\int_{0}^{1}I_{k-1}(x_{k},\dots,x_{d})dx_{k}-\sum_{i_{k}=0}^{n_{k}-1}I_{k-1}\left(x_{k}^{(i_{k})},x_{k+1},\dots,x_{d}\right)w_{k}^{(i_{k})}\right|\leq\varepsilon_{\mathrm{int}}.$
(8)

Then

$\left|I(f)-S_{X,W}(f)\right|\leq d\varepsilon_{\mathrm{int}}.$ (9)

###### Proof.

We will use induction to prove the following inequality for all $k=1,\dots,d$
and all $x_{k+1},\dots,x_{d}\in[0,1]$:

$\left|I_{k}(x_{k+1},\dots,x_{d})-\sum_{i_{1},\dots,i_{k}}f\left(x_{1}^{(i_{1})},\dots,x_{k}^{(i_{k})},x_{k+1},\dots,x_{d}\right)w_{1}^{(i_{1})}\dots
w_{k}^{(i_{k})}\right|\leq k\varepsilon_{\mathrm{int}}.$ (10)

The base case corresponds to $k=1$ and follows immediately from (8) as by
definition $I_{0}(x_{1},\dots,x_{d})\equiv f(x_{1},\dots,x_{d})$.

Now proceed to the induction step. To estimate the considered difference we
add and subtract the same value and then use the triangle inequality:

$\displaystyle\left|I_{k}(x_{k+1},\dots,x_{d})-\sum_{i_{1},\dots,i_{k}}f\left(x_{1}^{(i_{1})},\dots,x_{k}^{(i_{k})},x_{k+1},\dots,x_{d}\right)w_{1}^{(i_{1})}\dots
w_{k}^{(i_{k})}\right|$ $\displaystyle\leq$
$\displaystyle\leq\left|I_{k}(x_{k+1},\dots,x_{d})-\sum_{i_{k}=0}^{n_{k}-1}I_{k-1}\left(x_{k}^{(i_{k})},x_{k+1},\dots,x_{d}\right)w_{k}^{(i_{k})}\right|$
$\displaystyle+$
$\displaystyle+\left|\sum_{i_{k}=0}^{n_{k}-1}I_{k-1}\left(x_{k}^{(i_{k})},x_{k+1},\dots,x_{d}\right)w_{k}^{(i_{k})}-\sum_{i_{1},\dots,i_{k}}f\left(x_{1}^{(i_{1})},\dots,x_{k}^{(i_{k})},x_{k+1},\dots,x_{d}\right)w_{1}^{(i_{1})}\dots
w_{k}^{(i_{k})}\right|.$

The first summand is exactly the left part of (8) and thus not greater than
$\varepsilon_{\mathrm{int}}$. The second summand can be estimated as follows:

$\displaystyle
M:=\left|\sum_{i_{k}=0}^{n_{k}-1}I_{k-1}\left(x_{k}^{(i_{k})},x_{k+1},\dots,x_{d}\right)w_{k}^{(i_{k})}-\sum_{i_{1},\dots,i_{k}}f\left(x_{1}^{(i_{1})},\dots,x_{k}^{(i_{k})},x_{k+1},\dots,x_{d}\right)w_{1}^{(i_{1})}\dots
w_{k}^{(i_{k})}\right|$ $\displaystyle=$
$\displaystyle=\left|\sum_{i_{k}=0}^{n_{k}-1}\left(I_{k-1}\left(x_{k}^{(i_{k})},x_{k+1},\dots,x_{d}\right)-\sum_{i_{1},\dots,i_{k-1}}f\left(x_{1}^{(i_{1})},\dots,x_{k}^{(i_{k})},x_{k+1},\dots,x_{d}\right)w_{1}^{(i_{1})}\dots
w_{k-1}^{(i_{k-1})}\right)w_{k}^{(i_{k})}\right|$ $\displaystyle\leq$
$\displaystyle\leq\sum_{i_{k}=0}^{n_{k}-1}\left|I_{k-1}\left(x_{k}^{(i_{k})},x_{k+1},\dots,x_{d}\right)-\sum_{i_{1},\dots,i_{k-1}}f\left(x_{1}^{(i_{1})},\dots,x_{k}^{(i_{k})},x_{k+1},\dots,x_{d}\right)w_{1}^{(i_{1})}\dots
w_{k-1}^{(i_{k-1})}\right|w_{k}^{(i_{k})}.$

By induction hypothesis and the property (7) of the quadrature rule the
inequality chain can be continued:

$\displaystyle
M\leq\sum_{i_{k}=0}^{n_{k}-1}(k-1)\varepsilon_{\mathrm{int}}w_{k}^{(i_{k})}\leq(k-1)\varepsilon_{\mathrm{int}}\sum_{i_{k}=0}^{n_{k}-1}w_{k}^{(i_{k})}=(k-1)\varepsilon_{\mathrm{int}}.$

The induction step is thus proven. The conclusion of the lemma follows from
(10) for $k=d$. ∎

Let us discuss the applicability of the proven lemma to some classes of
functions. How restrictive is the property (8)? We will consider two examples.
First, if we use the $s$-point Gaussian quadrature rule and the function
$f(x_{1},\dots,x_{n})$ is a polynomial where each variable $x_{k}$ has degree
not greater than $2s-1$, then the quadrature rule gives the exact result:
$I(f)=S_{X,W}(f)$. Indeed, for every monomial $x_{1}^{p_{1}}\cdots
x_{d}^{p_{d}}$

$I_{k-1}(x_{k},\dots,x_{d})=\frac{1}{(1+p_{1})\cdots(1+p_{k-1})}x_{k}^{p_{k}}\cdots
x_{d}^{p_{d}}$

and

$\displaystyle\int_{0}^{1}I_{k-1}(x_{k},\dots,x_{d})dx_{k}=\frac{1}{(1+p_{1})\cdots(1+p_{k-1})}x_{k+1}^{p_{k+1}}\cdots
x_{d}^{p_{d}}\int_{0}^{1}x_{k}^{p_{k}}dx_{k}=$
$\displaystyle=\frac{1}{(1+p_{1})\cdots(1+p_{k-1})}x_{k+1}^{p_{k+1}}\cdots
x_{d}^{p_{d}}\sum_{i_{k}=0}^{n_{k}-1}\left(x_{k}^{(i_{k})}\right)^{p_{k}}w_{k}^{(i_{k})}=\sum_{i_{k}=0}^{n_{k}-1}I_{k-1}\left(x_{k}^{(i_{k})},x_{k+1},\dots,x_{d}\right)w_{k}^{(i_{k})},$

where the second equality follows from the fact that the Gaussian quadrature
rule is exact for polynomials of degree not greater than $2s-1$ (see e.g. [16,
Thm. 16.5.1]). By linearity the equality can be extended to the integral of
the polynomial $f(x_{1},\dots,x_{n})$.

In the second example we again use the $s$-point Gaussian quadrature rule and
the function $f$ has $2s$ continuous partial derivatives w.r.t. each variable
with a known bound:

$\left|\frac{\partial^{2s}}{\partial
x_{k}^{2s}}f(x_{1},\dots,x_{d})\right|\leq D.$

In this case the following inequality holds:

$\left|I(f)-S_{X,W}(f)\right|\leq\frac{\left(s!\right)^{4}}{(2s+1)\big{(}\left(2s\right)!\big{)}^{3}}\cdot
Dd.$

Indeed, by Leibniz integral rule (see e.g. [17]) we can move differentiation
under the integral sign and obtain the following estimate:

$\displaystyle\left|\frac{\partial^{2s}}{\partial
x_{k}^{2s}}\int_{0}^{1}\cdots\int_{0}^{1}f(x_{1},\dots,x_{d})dx_{1}\cdots
dx_{k-1}\right|$
$\displaystyle=\left|\int_{0}^{1}\cdots\int_{0}^{1}\frac{\partial^{2s}}{\partial
x_{k}^{2s}}f(x_{1},\dots,x_{d})dx_{1}\cdots dx_{k-1}\right|\leq$
$\displaystyle\leq\int_{0}^{1}\cdots\int_{0}^{1}\left|\frac{\partial^{2s}}{\partial
x_{k}^{2s}}f(x_{1},\dots,x_{d})\right|dx_{1}\cdots dx_{k-1}\leq D.$

It is known (e.g. see [15, §5.2]) that for some $\xi\in(0,1)$ the integration
error $\varepsilon_{\mathrm{int}}$ from (8) can be bounded by

$\frac{\left(s!\right)^{4}}{(2s+1)\big{(}\left(2s\right)!\big{)}^{3}}\left|\frac{\partial^{2s}}{\partial
x_{k}^{2s}}I_{k-1}(\xi,x_{k+1},\dots,x_{d})\right|\leq\frac{\left(s!\right)^{4}}{(2s+1)\big{(}\left(2s\right)!\big{)}^{3}}\cdot
D.$

The required bound for $|I(f)-S_{X,W}(f)|$ follows immediately.

It can be seen from the above examples that the property (8) is not extremely
restrictive and can sometimes be derived from the corresponding one-
dimensional results.

Now we can prove the main theorem.

######  Theorem 1.

If the conditions of lemma 1 hold and tensor
$A^{\prime}\in\mathbb{C}^{n_{1}\times\cdots\times n_{d}}$ satisfies the
inequality

$\left\|A_{X,f}-A^{\prime}\right\|_{F}\leq\varepsilon_{\mathrm{appr}},$

then

$\left|I(f)-\langle\overline{A^{\prime}},W\rangle\right|\leq
d\varepsilon_{\mathrm{int}}+\varepsilon_{\mathrm{appr}}.$

###### Proof.

We start by using the triangle property of modulus:

$\displaystyle\left|I(f)-\langle\overline{A^{\prime}},W\rangle\right|$
$\displaystyle\leq\left|I(f)-\langle\overline{A_{X,f}},W\rangle\right|+\left|\langle\overline{A_{X,f}},W\rangle-\langle\overline{A^{\prime}},W\rangle\right|=$
$\displaystyle=\left|I(f)-\langle\overline{A_{X,f}},W\rangle\right|+\left|\langle\overline{A_{X,f}-A^{\prime}},W\rangle\right|.$
(11)

The first term has been estimated in lemma 1 by $d\varepsilon_{\mathrm{int}}$.
The second summand can be bounded using the Cauchy–Bunyakovsky–Schwarz
inequality:

$\left|\langle\overline{A_{X,f}-A^{\prime}},W\rangle\right|\leq\sqrt{\left|\langle\overline{A_{X,f}-A^{\prime}},\overline{A_{X,f}-A^{\prime}}\rangle\right|}\cdot\sqrt{\left|\langle
W,W\rangle\right|}=\|A_{X,f}-A^{\prime}\|_{F}\|W\|_{F}.$

It suffices to estimate the value of $\|W\|_{F}$. Notice that by (7) all the
quadrature weights satisfy $0\leq w_{k}^{(i_{k})}\leq 1$, thus we can write

$\displaystyle\|W\|_{F}^{2}$
$\displaystyle=\sum_{i_{1},\dots,i_{d}}\prod_{k=1}^{d}\left(w_{k}^{(i_{k})}\right)^{2}=\prod_{k=1}^{d}\sum_{i_{k}=0}^{n_{k}-1}\left(w_{k}^{(i_{k})}\right)^{2}\leq\prod_{k=1}^{d}\sum_{i_{k}=0}^{n_{k}-1}w_{k}^{(i_{k})}=\prod_{k=1}^{d}1=1.$

∎

##  4 Implementation and experiments

We have implemented the described algorithm in C++. The source code is open
and distributed under the MIT license. The repository resides at
https://bitbucket.org/vysotskylev/c-tt-library. The implementation began as a
fork of tensor-train library by Dmitry Zheltkov
(https://bitbucket.org/zheltkov/c-tt-library). The only serious dependency is
OpenBLAS library (https://www.openblas.net) used for matrix operations.

Several examples of integration invocation are provided in
tests/test_integration.cpp. The key function is TTIntegrate. It takes as
parameters the integrand, the domain and configuration. The integrand is a
functor with parameters int count, double* points, double* result. The
integrand must be evaluated in count points with coordinates passed in points
array (which must be of size count * dimensionality) and the results must be
written to result array. The integration domain is supposed to be an axis-
aligned parallelepiped which is specified through the box parameter. The
config parameter allows the user to choose the limit on the number of function
evaluations, prescribed accuracy, type of quadrature etc. For the most recent
information on integration parameters and invocation refer to README.md file
in the repository.

By default the library uses multi-threading for matrix operations. However,
the calls to the functor are guaranteed to be sequential. The user is
encouraged to use parallel computations to evaluate the function in provided
points. The library strives to pass to the functor batches large enough to
keep multiple CPU cores (or even a GPU) busy. In the applications where
multiple integrals are to be computed the parallelization is even more
straightforward.

As a model example for integration consider the following integral with a
singularity at zero:

$\int_{0}^{1}\dots\int_{0}^{1}\ln(x_{1}\cdots x_{d})dx_{1}\cdots dx_{d}.$

The analytical answer is easily computed to be $-d$. We compare the proposed
algorithm with Monte-Carlo type method called Vegas implemented in the Cuba
library [26]. The plot of relative accuracy depending on the number of
variables is shown at Figure 1. The power-3 transform and Gauss-Legendre
quadrature were used for the TT integration. Both algortihms were limited by 1
million function evaluations.

Figure 1: Comparison of TT and Vegas integration on the model example

We added the new TT integrator as one of the integrators used in a development
version of the FIESTA [14] program used to evaluate Feynman integrals
numerically. Currently there exists a problem to numerically verify some
results for Feynman integrals111They appear as ingredients in quantum
chromodynamics (QCD) calculations and in supersymmetric $N=4$ Yang-Mills
theory within perturbation theory. These Feynman integrals are found in many
important functions of these theories, in particular, various anomalous
dimensions and form factors. obtained analytically. Most of them were checked
with the Vegas integrator, but we also used one of master integrals as a
cross-check for the TT integrator proposed in this paper. In the sector
decomposition approach this Feynman integral is represented as a sum of 1208
multivariate integrals of dimensions 5 or 6. These integrals were computed
numerically using TT and Vegas integration. The analytical result obtained
independently is a Laurent polynomial of variable $\varepsilon$ so we could
compare the accuracy of numerical integration methods. The formula for the
polynomial is

$\frac{1}{24}\varepsilon^{-3}+\frac{9}{16}\varepsilon^{-2}+\Big{(}\frac{451}{96}+\frac{5\pi^{2}}{72}\Big{)}\varepsilon^{-1}+\frac{18165+540\pi^{2}-1040\cdot\zeta(3)}{576}$

(here $\zeta(s)$ denotes Riemann Zeta function). The integrators (TT and
Vegas) were allowed to evaluate the integrand in 1 million points. For TT
integration we used Gauss-Legendre quadrature with 13 nodes on each axis and
performed power transform (see Section 3.5) of form $x=t^{3}$ as the integrand
has a logarithmic singularity at zero. The relative errors for each term are
summarized in Table 1.

| $\varepsilon^{-3}$ | $\varepsilon^{-2}$ | $\varepsilon^{-1}$ | $\varepsilon^{0}$
---|---|---|---|---
TT | $27.3\times 10^{-6}$ | $0.3\times 10^{-6}$ | $0.06\times 10^{-6}$ | $0.88\times 10^{-6}$
Vegas | $188.6\times 10^{-6}$ | $7.1\times 10^{-6}$ | $0.11\times 10^{-6}$ | $3.96\times 10^{-6}$
Table 1: Relative errors of numerical integration for a Feynman integral

The results are quite promising and show that in some circumstances the new
integrator can outperform the Vegas algorithm. We are working on a new release
of the FIESTA program where the new TT integrator might become one of the key
components.

##  5 Acknowledgements

The work is supported by Russian Ministry of Science and Higher Education,
agreement No. 075-15-2019-1621.

## References

  * [1] C. Bogner and S. Weinzierl, Resolution of singularities for multi-loop integrals, Comput.Phys.Commun., 178, 596–610 (2008). https://doi.org/10.1016/j.cpc.2007.11.012
  * [2] R. Piessens, E. de Doncker, C. Uberhuber and D. Kahaner, Quadpack: a subroutine package for automatic integration (Springer-Verlag, 1983). https://doi.org/10.1007/978-3-642-61786-7
  * [3] G P. Lepage, A new algorithm for adaptive multidimensional integration, Journal of Computational Physics, 27(2), 192–203 (1978). https://doi.org/10.1016/0021-9991(78)90004-9
  * [4] G. P. Lepage, VEGAS: An adaptive multi-dimensional integration routine, Technical Report CLNS-80/447, Newman Laboratory of Nuclear Studies, Cornell University, Ithaca, NY, 1980.
  * [5] W. H. Press and G. R. Farrar, Recursive stratified sampling for multidimensional Monte Carlo integration, Comp. Phys., 4(2), 190–195 (1990). https://doi.org/10.1063/1.4822899
  * [6] R. Schürer, Adaptive Quasi-Monte Carlo Integration Based on MISER and VEGAS, In: Niederreiter H. (eds) Monte Carlo and Quasi-Monte Carlo Methods 2002. (Springer, Berlin, Heidelberg. 2004). https://doi.org/10.1007/978-3-642-18743-8_25
  * [7] J. Berntsen, T. Espelid and A. Genz, An adaptive algorithm for the approximate calculation of multiple integrals, ACM Trans. Math. Software, 17, 437–451 (1991). https://doi.org/10.1145/210232.210233
  * [8] J. Dick, F. Y. Kuo and I. H. Sloan, High-dimensional integration: The quasi-monte carlo way, Acta Numerica, 22, 133–288 (2013).
  * [9] S. Borowka, G. Heinrich, S. Jahn, S. P. Jones, M. Kerner, J. Schlenk and T. Zirke, pySecDec: a toolbox for the numerical evaluation of multi-scale integrals, Comput. Phys. Commun., 222 313–326 (2018). https://doi.org/10.1016/j.cpc.2017.09.015
  * [10] S. Borowka, G. Heinrich, S. Jahn, S. P. Jones, M. Kerner and J. Schlenk, A GPU compatible quasi-Monte Carlo integrator interfaced to pySecDec, Comput. Phys. Commun., 240, 120–137 (2019). https://doi.org/10.1016/j.cpc.2019.02.015.
  * [11] I. V. Oseledets and E. E. Tyrtyshnikov, TT-cross approximation for multidimensional arrays, Linear Algebra and Its Applications, 432(1), 70–88 (2010). https://doi.org/10.1016/j.laa.2009.07.024
  * [12] S. Dolgov and D. Savostyanov, Parallel cross interpolation for high-precision calculation of high-dimensional integrals. Computer Physics Communications, 246, 106869 (2020). https://doi.org/10.1016/j.cpc.2019.106869
  * [13] I. V. Oseledets, Tensor-Train Decomposition, SIAM Journal on Scientific Computing, 33(5), 2295–2317 (2011). https://doi.org/10.1137/090752286
  * [14] A. V. Smirnov, FIESTA 4: optimized Feynman integral calculations with GPU support, Comp. Phys. Comm, 204, 189–199 (2016). https://doi.org/10.1016/j.cpc.2016.03.013
  * [15] D. Kahaner, C. Moler, S. Nash, Numerical methods and software (Prentice–Hall, Inc., USA, 1989).
  * [16] E. E. Tyrtyshnikov, A brief introduction to numerical analysis (Springer Science & Business Media, 2012).
  * [17] M. H. Protter and B. Charles Jr, Intermediate calculus (Springer Science & Business Media, 2012).
  * [18] S. A. Goreinov and E. E. Tyrtyshnikov, The maximal-volume concept in approximation by low-rank matrices, Contemporary Mathematics, 208, 47–52 (2001)
  * [19] S. A. Goreinov, N. L. Zamarashkin and E.E Tyrtyshnikov, Pseudo-skeleton approximations by matrices of maximal volume, Mathematical Notes, 62(4), 515–519 (1997). https://doi.org/10.1007/bf02358985
  * [20] A. Civril and M. Magdon-Ismail, On selecting a maximum volume sub-matrix of a matrix and related problems. Theoretical Computer Science, 410(47–49), 4801–4811 (2009). https://doi.org/10.1016/j.tcs.2009.06.018
  * [21] M. Bebendorf, Approximation of boundary element matrices, Numerische Mathematik, 86, 565–589 (2000). https://doi.org/10.1007/PL00005410
  * [22] E. E. Tyrtyshnikov, Incomplete Cross Approximation in the Mosaic-Skeleton Method, Computing (Vienna/New York), 64(4), 367–380 (2000). https://doi.org/10.1007/s006070070031
  * [23] S. A. Goreinov, I. V. Oseledets, D. V. Savostyanov, E. E Tyrtyshnikov and N. L. Zamarashkin, How to Find a Good Submatrix, Matrix Methods: Theory, Algorithms and Applications, 247–256 (2010). https://doi.org/10.1142/9789812836021_0015
  * [24] S. M. Kirkup, J. Yazdani and G. Papazafeiropoulos, Quadrature Rules for Functions with a Mid-Point Logarithmic Singularity in the Boundary Element Method Based on the $x=t^{p}$ Substitution, American Journal of Computational Mathematics, 09(04), 282–301 (2019). https://doi.org/10.4236/ajcm.2019.94021
  * [25] M. Mori and M. Sugihara, The double-exponential transformation in numerical analysis. Journal of Computational and Applied Mathematics, 287–296 (2001). https://doi.org/10.1016/S0377-0427(00)00501-X
  * [26] T. Hahn, Cuba - a library for multidimensional numerical integration, Comput. Phys. Commun., 168 78–95 (2005) https://doi.org/10.1016/j.cpc.2005.01.010