Pseudo-hermitian matrices are matrices hermitian with respect to an indefinite metric. They can be thought of as the truncation of pseudo-hermitian operators, defined over some Krein space, together with the associated metric, to a finite dimensional subspace. As such, they can be used, in the usual spirit of random matrix theory, to model chaotic or disordered PT-symmetric quantum systems, or their gain-loss-balanced classical analogs, in the phase of broken PT-symmetry. The eigenvalues of pseudo-hermitian matrices are either real, or come in complex-conjugate pairs. In this paper we introduce a family of pseudo-hermitian random matrix models, depending parametrically on their metric. We apply the diagrammatic method to obtain its averaged resolvent and density of eigenvalues as explicit functions of the metric, in the limit of large matrix size $N$. As a concrete example, which is essentially an ensemble of elements of the non-compact unitary Lie algebra, we choose a particularly simple set of metrics, and compute the resulting resolvent and density of eigenvalues in closed form. The spectrum consists of a finite fraction of complex eigenvalues, which occupy uniformly two two-dimensional blobs, symmetric with respect to the real axis, as well as the complimentary fraction of real eigenvalues, condensed in a finite segment, with a known non-uniform density. The numbers of complex and real eigenvalues depend on the signature of the metric, that is, the numbers of its positive and negative eigenvalues. We have also carried thorough numerical analysis of the model for these particular metrics. Our numerical results converge rapidly towards the asymptotic analytical large-$N$ expressions.