If your matrices happen to be low-rank, you might want to consider a simultaneous low-rank approximation of your matrices (as in https://doi.org/10.1109/CVPR.2006.112), i.e. you compute $$ \min_{A,M_i,B} \sum_{i=1}^p || S_i - A M_i B ||_F^2, $$ where $M_i$ are potentially smaller than $S_i$. You could then compute $$ S \approx A \big( \sum_{i=1}^p M_i \big) B.$$

I am however not aware of any result on the error of $S$. I would suspect that the error is equal to the sum of the individual low-rank approximation errors.

If your matrices happen to be low-rank, you might want to consider a simultaneous low-rank approximation of your matrices (as in Inoue and Urahama - Equivalence of Non-Iterative Algorithms for Simultaneous Low Rank Approximations of Matrices), i.e. you compute $$ \min_{A,M_i,B} \sum_{i=1}^p \lVert S_i - A M_i B \rVert_F^2, $$ where $M_i$ are potentially smaller than $S_i$. You could then compute $$ S \approx A \bigl( \sum_{i=1}^p M_i \bigr) B.$$

I am however not aware of any result on the error of $S$. I would suspect that the error is equal to the sum of the individual low-rank approximation errors.