Heterosis describes the phenomenon whereby a hybrid population has higher fitness than an inbred population, and has previously been explained by either Mendelian dominance or overdominance, where it is generally assumed that one gene controls one trait. However, recent studies have demonstrated that genes interact through a complex gene regulatory network (GRN). Furthermore, phenotypic variance due to noise is reportedly lower for heterozygotes, whereas the origin of such variance-related heterosis remains elusive. Therefore, a theoretical analysis linking heterosis to GRN evolution and stochastic gene expression dynamics is required. Here, we investigate heterosis related to fitness and phenotypic variance in a system with interacting genes, by numerically evolving diploid GRNs. According to the results, the heterozygote population exhibited higher fitness than the homozygote population, that is, fitness-related heterosis resulting from evolution. In addition, the heterozygote population expressed lower noise-related phenotypic variance in expression levels than the homozygous population, implying that the heterozygote population is more robust to noise. Furthermore, the distribution of the ratio of heterozygote phenotypic variance to homozygote phenotypic variance exhibited quantitative agreement with previous experimental results. By applying dominance and overdominance to the gene expression pattern rather than only a single gene expression, we confirmed the correlation between heterosis and overdominance. We explain our results by proposing that the convex high-fitness region is evolutionarily shaped in the genetic space to gain noise robustness under genetic mixing through sexual reproduction.