Our paper presents a novel high dimensional probability density estimation technique using any dimensionality reduction method. Our method first performs subspace reduction using any matrix factorization algorithm and estimates the density in the low-dimensional space using sample-point variable bandwidth kernel density estimation. Subsequently, the high dimensional density is approximated from the low dimensional density parameters. The reconstruction error due to dimensionality reduction process is also modeled in a principled and efficient manner to obtain the high dimensional density estimate. We show the effectiveness of our technique by using two popular dimensionality reduction tools, principal component analysis and non-negative matrix factorization. This technique is applied to AT&T, Yale, Pointing'04 and CMU-PIE face recognition datasets and improved performance compared to other dimensionality reduction and density estimation algorithms is obtained.