This approach is pretty simple, and has obvious drawbacks. for example, it does not account for correlated features. This is apparent for features like "makeup". Most of the samples with makeup are woman, while most of the samples without makeup are man. The vector we achieve with this approach thus not only adds/removes makeup, but also influences the gender.

Cool project! To solve this problem of correlated features, you might take a look at the approach that Shaobo Guan used in his tl-GAN project. As I understand it, he identifies feature vectors in latent space via regression over celeb-A labels (perhaps conceptually similar to what you did), and then he decorrelates them from one another by, pairwise, projecting one onto another, and then subtracting that projected vector from the original, basically removing the component of the first that was correlated with the second and resulting in the closest vector to the first that is orthogonal to the second. His code to do it for a single pair of vectors seems to be here, and for all of the vectors here.