Scaling tensor approximation by symmetric tensor decomposition with SciPy's L-BFGS-B

I am trying to approximate a symmetric tensor of which the values are in the range of [1e-7,1e-4], by a symmetric tensor decomposition of lower rank. For this I am using the L-BFGS-B method in SciPy’s optimization package.

The relative errors (based on the Frobenius norm) I am obtaining are quite big (greater than 1). After some research I have come to the conclusion that there is need for scaling, as my tensor is ‘poorly scaled’. When I multiply all values of X with 1e7, I do obtain a smaller relative error (in the range of [1e-4,1e-3]), but only when the modes of X have a small dimension in comparison with the chosen rank of the symmetric tensor decomposition.

Since I am not that familiar with the concept of scaling in numerical optimization, is there a better way to solve this scaling problem, so that I can obtain a small relative error even when the dimension of the modes of X is large in comparison with the chosen rank?

I am also doing some research about the existence of a low rank approximation, as this may not even exist. While this may be the case, reproducing the first experiment of Sherman & Kolda (2020), does not give me the same order of magnitude of the relative errors. This means that there should exist some improvements to be made to my implementation, of which I think only the scaling aspect can be improved.