Network function virtualization (NFV) is an emerging design paradigm that replaces physical middlebox devices with software modules running on general purpose commodity servers. While gradually transitioning to NFV, Internet service providers face the problem of where to introduce NFV in order to make the most benefit of that; here, we measure the benefit by the amount of traffic that can be serviced through the NFV. This problem is non-trivial as it is composed of two challenging subproblems: 1) placement of nodes to support virtual network functions (referred to as VNF-nodes); and 2) allocation of the VNF-nodes resources to network flows; the two subproblems need to be considered jointly to satisfy the objective of serving the maximum amount of traffic. This problem has been studied recently but for the one-dimensional setting, where all network flows require one network function, which requires a unit of resource to process a unit of flow. In this work, we extend to the multi-dimensional setting, where flows can require multiple network functions, which can also require a different amount of each resource to process a unit of flow. The multi-dimensional setting introduces new challenges in addition to those of the onedimensional setting (e.g., NP-hardness and non-submodularity) and also makes the resource allocation a multi-dimensional generalization of the generalized assignment problem with assignment restrictions. To address these difficulties, we propose a novel two-level relaxation method and utilize the primal-dual technique to design two approximation algorithms that achieve an approximation ratio of (e−1)(Z−1) 2e 2 Z(kR) 1/(Z−1) and (e−1)(Z−1) 2e(Z−1+eZR 1/(Z−1) ) , where k (resp. R) is the number of VNF-nodes (resp. resources), and Z is a measure of the available resource compared to flow demand. Finally, we perform extensive trace-driven simulations to show the effectiveness of the proposed algorithms.