No problem, glad that functionality is being used and seems to be working (I wrote it mostly for internal usage in the package but I think it is useful more generally).
Getting a global error of around 6e-3 using a cutoff of 1e-5 doesn't sounds so unreasonable, since the cutoff is used at every bond of the MPS (so the total error would be something like the sum of the truncation error at each bond). Depending on the size of your MPS, the global error measured by `norm(A - prod(psi_optimised))` can be larger than the cutoff you use, as you were seeing.
Also note that the `isapprox` function has a few keyword arguments for setting the tolerance for the approximation, for example `isapprox(A, prod(psi); rtol = 1e-5)` (where `rtol` sets a relative tolerance scaled by the norms of the input tensors). A reasonable choice for `rtol` may be something like `isapprox(A, prod(psi); rtol = cutoff * N)` where `cutoff` is the cutoff you set and `N` is the system size. You would have to play around with this, of course. Doing floating point comparisons is always a bit tricky.