This work presents a theoretical model of boson sampling with optical feedback, in which a subset of the interferometer’s output modes is looped back into the inputs and fresh photons are injected at successive time steps. We develop complementary methods for computing finite-time output distributions and the state retained in the loop, including spatiotemporal unfolding, direct evolution of the reduced density matrix, a Kraus-operator description, and a correlation-tensor formalism. Using the optical transformation law of characteristic functions, we prove that the loop approaches a unique stationary state whenever the loop-to-loop block (U_{LL}) has spectral radius smaller than one; the exceptional unit-modulus case is shown to correspond to a perfectly nonleaking loop mode. The methods include photon loss and variable photon number, and the numerical results illustrate the stabilization timescale, stationary photon statistics, and reconstruction of the stationary density matrix from correlation tensors. Finally, we formulate stationary-distribution boson sampling, in which fresh photons interfere with the stationary optical memory, and identify the complexity of scalable recurrent instances as a direction for future research.