Linear discriminant analysis (LDA) is a commonly used method for dimensionality reduction. Despite its successes, it has limitations under some situations, including the small sample size problem, the homoscedasticity assumption that different classes have the same Gaussian distribution, and its inability to produce probabilistic output and handle missing data. In this paper, we propose a semi-supervised and heteroscedastic extension of probabilistic LDA, called S$^2$HPLDA, which aims at overcoming all these limitations under a common principled framework. Moreover, we apply automatic relevance determination to determine the required dimensionality of the low-dimensional space for dimensionality reduction. We empirically compare our method with several related probabilistic subspace methods on some face and object databases. Very promising results are obtained from the experiments showing the effectiveness of our proposed method.