  Recently, a tensor-on-tensor (ToT) regression model has been proposed to
generalize tensor recovery, encompassing scenarios like scalar-on-tensor
regression and tensor-on-vector regression. However, the exponential growth in
tensor complexity poses challenges for storage and computation in ToT
regression. To overcome this hurdle, tensor decompositions have been
introduced, with the tensor train (TT)-based ToT model proving efficient in
practice due to reduced memory requirements, enhanced computational efficiency,
and decreased sampling complexity. Despite these practical benefits, a
disparity exists between theoretical analysis and real-world performance. In
this paper, we delve into the theoretical and algorithmic aspects of the
TT-based ToT regression model. Assuming the regression operator satisfies the
restricted isometry property (RIP), we conduct an error analysis for the
solution to a constrained least-squares optimization problem. This analysis
includes upper error bound and minimax lower bound, revealing that such error
bounds polynomially depend on the order $N+M$. To efficiently find solutions
meeting such error bounds, we propose two optimization algorithms: the
iterative hard thresholding (IHT) algorithm (employing gradient descent with
TT-singular value decomposition (TT-SVD)) and the factorization approach using
the Riemannian gradient descent (RGD) algorithm. When RIP is satisfied,
spectral initialization facilitates proper initialization, and we establish the
linear convergence rate of both IHT and RGD.
