Use the additional output err to automatically obtain the fit in terms of sum squared error and use the input Options(1) to set the convergence criterion to 1e-8.
Estimate the model several times and compare fit values:
Options=[];
Options(1)=1e-7;
% We need to start from different places every time Options(2)=2;
[Factors1,it1,err1] = parafac(X,5,Options);
[Factors2,it2,err2] = parafac(X,5,Options);
Note, however, that PARAFAC by default uses direct trilinear decomposition for initialization of three-way arrays (and SVD initialization if there are missing values).
This means that repeating the analysis several times will give the exact same result since the initial parameters are identical.
Set the initialization method to random orthogonal loadings(Options(2)=2;) to avoid this