We apply the vectorized Non-negative Matrix Factorization (NMF) method to post-processing of direct imaging data for exoplanetary systems such as circumstellar disks. NMF is an iterative approach, which first creates a non-orthogonal and non-negative basis of components using given reference images, then models a target with the components. The constructed model is then rescaled with a factor to compensate for the contribution from a disk. We compare NMF with existing methods (classical reference differential imaging method, and the Karhunen-Lo\`eve image projection algorithm) using synthetic circumstellar disks, and demonstrate the superiority of NMF: with no need for prior selection of references, NMF can detect fainter circumstellar disks, better preserve low order disk morphology, and does not require forward modeling. As an application to a well-known disk example, we process the archival Hubble Space Telescope (HST) STIS coronagraphic observations of HD~181327 with different methods and compare them. NMF is able to extract some circumstellar material inside the primary ring for the first time. In the appendix, we mathematically investigate the stability of NMF components during iteration, and the linearity of NMF modeling.