Prescribing arbitrary stress-strain behavior in soft metamaterials is of broad interest for various applications, such as soft robotics, energy dissipation, shock absorption, soft wearables, and biomedical implants. In this regard, the state-of-the-art approaches include inverse design methods using topology optimization[1,2] and machine learning[3] that finds the optimized material distribution within a specified design domain. However, these methods often struggle to achieve highly nonlinear target responses that involve multiple snap-through or snap-back behavior. These challenges stem from the high-dimensional, non-convex nature of the inverse problem and the limited design resolution imposed by practical computational constraints. In this study, we propose a two-step inverse optimization framework that overcomes these limitations by controlling the collective response of a finite assembly of heterogeneous unit designs, each exhibiting comparatively simpler mechanical behavior. In the first step, the unit responses are parameterized in a reduced dimensional space and optimized to achieve a prescribed complex assembly-level response. In the second step, topology optimization is employed to determine the corresponding material distributions for the optimized unit responses. Using this approach, we demonstrate realization of complex stress-strain assembly responses featuring multiple plateau, snap-buckling, and multi-step energy dissipation. The proposed design framework is further validated through fabrication and experimental testing of selected optimized designs. Overall, this method provides a practical and scalable route for designing soft metamaterials with heterogeneous units and arbitrarily prescribed nonlinear mechanical responses.