We design an importance sampling scheme for backward stochastic differential equations (BSDEs) that minimizes the conditional variance occurring in least-squares Monte Carlo (LSMC) algorithms. The Radon-Nikodym derivative depends on the solution of BSDE, and therefore it is computed adaptively within the LSMC procedure. To allow robust error estimates with respect to the unknown change of measure, we properly randomize the initial value of the forward process. We introduce novel methods to analyze the error: firstly, we establish norm stability results due to the random initialization; secondly, we develop refined concentration-of-measure techniques to capture the variance of reduction. Our theoretical results are supported by numerical experiments.