  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.
