Nonlinear parametric inverse problems appear in many applications. In this work, we focus on diffuse optical tomography (DOT) in medical image reconstruction, in which we aim to recover an unknown image of interest (defined by a set of parameters), such as cancerous tissue in a given medium, using a mathematical (forward) model. The forward model in DOT is a diffusion-absorption model for the photon flux. The main computational bottleneck in such inverse problems is the repeated evaluation of the large-scale forward model. For DOT, this corresponds to solving large linear systems for each source and frequency at each optimization step. In addition, as Newton-type methods are the method of choice for these problems, one needs to solve additional linear systems with the adjoint to efficiently compute derivative information. As rapid advances in technology allow for large numbers of sources and detectors, these problems become computationally prohibitively expensive. In the past, the use of reduced order models (ROM) has been proposed to drastically reduce the size of the linear systems solved in each optimization step, while still solving the inverse problem accurately. However, as the number of sources and detectors increases, even the construction of the ROM bases incurs a substantial cost. Interpolatory model reduction, when performed in the setting of full matrix-interpolation to match the full parameter gradient matrix, requires the solution of large linear systems for all sources and frequencies as well as for all detectors and frequencies for each interpolation point in parameter space to build a candidate basis for the ROM projection space. As this candidate basis numerically has low rank, this construction is followed by a rank-revealing factorization that typically reduces the number of vectors in the candidate basis substantially. Since only the ROM projection space matters, not its basis, we propose to use randomization to approximate this basis efficiently and drastically reduce the number of large linear solves for constructing the global ROM basis for DOT. We also provide a detailed analysis for the low-rank structure of the candidate basis for our problem of interest. Even though we focus on the DOT problem here, the ideas presented are relevant to many other large scale inverse problems and optimization problems.