Nonlinear local Lyapunov vectors (NLLVs), theoretically inherited from the linear Lyapunov vectors (LVs) in a nonlinear framework, are developed to indicate orthogonal directions in phase space with different perturbation growth rates. Practically NLLVs generates flow dependent perturbations with full nonlinear models as in the breeding method, but regularly conducts the Gram–Schmidt reorthonormalization processes on the perturbations. The advantages of the sampling of the unstable subspace with the mutually orthogonal NLLVs instead of the most unstable direction with the dependent bred vectors (BVs) are clarified by using their respective unstable mode to predict the structure of forecast error growth. Additionally, the NLLVs have overall higher but correlated local dimensions compared to the BVs which may be beneficial for the former to increase the ensemble spread and capture the instabilities as well. The NLLV method is used to generate initial perturbations and implement ensemble forecasts in nonlinear models (the Lorenz63 model, Lorenz96 model, a barotropic quasi-geostrophic model and the Zebiak-Cane model) to explore the validity of the NLLV method. The performance of the NLLV method is compared comprehensively and systematically with other methods such as the ensemble transform Kalman filter (ETKF) method, the BV and the random perturbation (Monte Carlo) methods. Overall, the ensemble spread of NLLVs is more consistent with the errors of the ensemble mean, which indicates the better performance of NLLVs in simulating the evolution of analysis errors. In addition, the NLLVs perform significantly better than the BVs in terms of reliability and the random perturbations in resolution. The NLLV scheme has slightly higher ensemble forecast skill than the ETKF scheme. The NLLV scheme has significantly simpler algorithm and less computation time than the ETKF.