This paper develops a general theory of optimal income taxation with multiple dimensions of agent heterogeneity. The main technical hurdle in developing this theory is the possibility that individuals have multiple optimal incomes. Using a perturbation approach, optimal tax formulas are derived that account for the possibility that individuals have multiple optima and, hence, account for the possibility that individuals jump between their optimal income levels when the tax schedule is perturbed. The magnitude of these effects is quantified, thereby augmenting the optimal tax formulas from Saez (2001) with additional "jumping effect" terms. The paper provides a partial characterization of when individuals with multiple optimal incomes may exist under the optimal tax schedule. Finally, the paper derives a new methodology to simulate optimal income tax schedules with multidimensional heterogeneity. This method is implemented numerically, showing that individuals with multiple optimal income levels can exist under the optimal tax schedule.