Several parameters entering the constitutive laws of metal matrix composites are not amenable to direct measurement. Their determination is therefore demanded to indirect identification procedures, based on the comparison between the output of some experimental works and the results of the iterative simulation of the performed tests. Largely repetitive numerical analyses can also support the optimized design of composite microstructures, the understanding of the actual stress-transfer mechanisms, and virtual testing. Surrogate analytical models can replace effectively non-linear finite element (FE) computations in parametric studies and identification procedures, reducing execution times and costs without compromising the accuracy of the results. In this context, a frequently used approach consists of the interpolation, by means of radial basis functions (RBFs), of data filtered by proper orthogonal decomposition (POD). This methodology is often used to reproduce the smooth output of different mechanical systems. This work is rather focused on the homogenized response of damaging metal matrix composites subjected to strain localization phenomena, and explores the accuracy achievable in these situations by the POD-RBF approximation of more traditional FE analyses. In all considered cases, the POD-RBF results are accurate, and able to distinguish between apparently similar situations. The approach is flexible, and the performance of the surrogate models can be tailored to the requirements of each application. In particular, various analytical approximations can be introduced to support the design of new microstructures and material couplings, and to understand the role of any material and/or geometry imperfections resulting from the production processes.