This paper discusses the extension of Probabilistic Independent Component Analysis (PICA) to multi-subject and multi-session functional MRI (fMRI) data, leading to a three-way decomposition that captures temporal, spatial, and subject-dependent variations. The technique, called Tensorial PICA (tensor-PICA), is derived from Parallel Factor Analysis (PARAFAC) and is compared with PARAFAC using simulated and real fMRI data. Tensor-PICA is shown to be more accurate, reduce cross-talk between estimated sources, and provide robustness against deviations from model assumptions and overfitting. It also offers computational speed advantages. On real fMRI data, tensor-PICA extracts plausible activation maps, time courses, and session/subject modes, providing a rich description of additional processes such as image artifacts or secondary activation patterns. This approach enhances the interpretation and optimization of group fMRI studies beyond what can be achieved with model-based analysis techniques.This paper discusses the extension of Probabilistic Independent Component Analysis (PICA) to multi-subject and multi-session functional MRI (fMRI) data, leading to a three-way decomposition that captures temporal, spatial, and subject-dependent variations. The technique, called Tensorial PICA (tensor-PICA), is derived from Parallel Factor Analysis (PARAFAC) and is compared with PARAFAC using simulated and real fMRI data. Tensor-PICA is shown to be more accurate, reduce cross-talk between estimated sources, and provide robustness against deviations from model assumptions and overfitting. It also offers computational speed advantages. On real fMRI data, tensor-PICA extracts plausible activation maps, time courses, and session/subject modes, providing a rich description of additional processes such as image artifacts or secondary activation patterns. This approach enhances the interpretation and optimization of group fMRI studies beyond what can be achieved with model-based analysis techniques.