The synthesis of high-resolution (HR) hyperspectral image (HSI) by fusing a low-resolution HSI with a corresponding HR multispectral image has emerged as a prevalent HSI super-resolution (HSR) scheme. Recent researches have revealed that tensor analysis is an emerging tool for HSR. However, most off-the-shelf tensor-based HSR algorithms tend to encounter challenges in rank determination and modeling capacity. To address these issues, we construct nonlocal patch tensors (NPTs) and characterize low-rank structures with coupled Bayesian tensor factorization. It is worth emphasizing that the intrinsic global spectral correlation and nonlocal spatial similarity can be simultaneously explored under the proposed model. Moreover, benefiting from the technique of automatic relevance determination, we propose a hierarchical probabilistic framework based on Canonical Polyadic (CP) factorization, which incorporates a sparsity-inducing prior over the underlying factor matrices. We further develop an effective expectation-maximization-type optimization scheme for framework estimation. In contrast to existing works, the proposed model can infer the latent CP rank of NPT adaptively without tuning parameters. Extensive experiments on synthesized and real datasets illustrate the intrinsic capability of our model in rank determination as well as its superiority in fusion performance.