This paper develops the Non-Intrusive Least Squares Shadowing (NILSS) method, which computes the sensitivity for long-time averaged objectives in chaotic dynamical systems. In NILSS, we represent a tangent solution by a linear combination of one inhomogeneous tangent solution and several homogeneous tangent solutions. Next, we solve a least squares problem using this representation; thus, the resulting solution can be used for computing sensitivities. NILSS is easy to implement with existing solvers. In addition, for chaotic systems with many degrees of freedom but few unstable modes, NILSS has a low computational cost. NILSS is applied to two chaotic PDE systems: the Lorenz 63 system and a CFD simulation of flow over a backward-facing step. In both cases, the sensitivities computed by NILSS reflect the trends in the long-time averaged objectives of dynamical systems.where f (u, s) : R m × R → R m is a smooth function, u is the state, and s is the parameter. The initial condition u 0 is a smooth function of φ. A solution u(t) is called the primal solution.In this paper, The objective is a long-time averaged quantity. To define it, we first let J(u, s) : R m × R → R be a continuous function that represents the instantaneous objective function. The objective is obtained by averaging J over a infinitely long trajectory:J T depends on s, φ, and T , while J ∞ is determined only by s and u 0 . Here we make the assumption of ergodicity [3], which means that u 0 , hence φ does not affect J ∞ . As a result, J ∞ only depends on s. The purpose of this paper is to develop an algorithm that computes the sensitivity d J ∞ /ds. The sensitivity can help scientists and engineers design products [4,5], control processes and systems [6,7], solve inverse problems [8], estimate simulation errors [9,10,11], assimilate measurement data [12,13] and quantify uncertainties [14]. When a dynamical system is chaotic, computing a meaningful d J ∞ /ds is challenging, since in general:That is, if we fix u 0 (φ), the process of T → ∞ does not commute with differentiation with respect to s [15]. As a result, the transient method, which employs the conventional tangent method with a fixed u 0 , does not converge to the correct sensitivity for chaotic systems. In fact, the transient method diverges most of the time [15]. Many sensitivity analysis methods have been developed to compute d J ∞ /ds. The conventional methods include the finite difference and transient method. The ensemble method, developed by Lea et al. [16,17], computes the sensitivity by averaging results from the transient method over an ensemble of trajectories. Another recent approach is based on the fluctuation dissipation theorem (FDT), as seen in [18,19,20,21,22,23].In this research study, we consider the Least Squares Shadowing (LSS) approach, developed by Wang, Hu and Blonigan [24,14]. LSS computes a bounded shift of a trajectory under an infinitesimal parameter change, which is called the LSS solution. The LSS solution can then be used to compute the derivative d J ∞ /ds. LSS...