Accurate approximations of the change of system's output and its statistics with respect to the input are highly desired in computational dynamics. Ruelle's linear response theory provides breakthrough mathematical machinery for computing the sensitivity of chaotic dynamical systems, which enables a better understanding of chaotic phenomena. In this paper, we propose an algorithm for sensitivity analysis of discrete chaos with an arbitrary number of positive Lyapunov exponents. We combine the concept of perturbation space-splitting regularizing Ruelle's original expression together with measure-based parameterization of the expanding subspace. We use these tools to rigorously derive trajectory-following recursive relations that exponentially converge, and construct a memory-efficient Monte Carlo scheme for derivatives of the output statistics. Thanks to the regularization and lack of simplifying assumptions on the behavior of the system, our method is immune to the common problems of other popular systems such as the exploding tangent solutions and unphysicality of shadowing directions. We provide a ready-to-use algorithm, analyze its complexity, and demonstrate several numerical examples of sensitivity computation of physically-inspired low-dimensional systems.