Genetic algorithm (GA) is combined with finite-element method for the first time as an alternative method to invert gravity anomaly data for reconstructing the 3D density structure in the subsurface. The computational efficiency is significantly improved by storing the coefficient matrix and using it in all the forward calculations, and by dividing the region of interest into many sub-regions and applying the GA to the sub-regions in parallel. Central Taiwan - a geologically complex region - is used as an example to demonstrate the utility of the method. A crustal block, 120 x 150 km2 in area and 34 km in thickness, is represented by a finite-element model of 76,500 cubic elements, each 2x2x2 km3 in size. An initial density model is reconstructed from the regional 3D tomographic seismic velocity, using an empirical relation between velocity and density. The difference between the calculated and the observed gravity anomaly (i.e., the residual anomaly) shows an elongated minimum of large magnitude that extends along the axis of the Taiwan mountain belt. Among the interpretative models tested, the best model shows a crustal root beneath the axis of the Central Ranges and a density contrast of 400 or 500 kg/m3 across the Moho. Both predictions appear to be supported by independent seismological and laboratory evidences.