To characterize circumstellar systems in high contrast imaging, the fundamental step is to construct a best point spread function (PSF) template for the non-circumstellar signals (i.e., star light and speckles) and separate it from the observation. With existing PSF construction methods, the circumstellar signals (e.g., planets, circumstellar disks) are unavoidably altered by over-fitting and/or self-subtraction, making forward modeling a necessity to recover these signals. We present a forward modeling--free solution to these problems with data imputation using sequential non-negative matrix factorization (DI-sNMF). DI-sNMF first converts this signal separation problem to a "missing data" problem in statistics by flagging the regions which host circumstellar signals as missing data, then attributes PSF signals to these regions. We mathematically prove it to have negligible alteration to circumstellar signals when the imputation region is relatively small, which thus enables precise measurement for these circumstellar objects. We apply it to simulated point source and circumstellar disk observations to demonstrate its proper recovery of them. We apply it to Gemini Planet Imager (GPI) K1-band observations of the debris disk surrounding HR 4796A, finding a tentative trend that the dust is more forward scattering as the wavelength increases. We expect DI-sNMF to be applicable to other general scenarios where the separation of signals is needed.