This paper is devoted to the numerical treatment of linear optimality systems (OS) that arise in connection with inverse problems for partial differential equations. If such inverse problems are regularized by Tikhonov regularization, then it follows from standard theory that the associated OS is well-posed, provided that the regularization parameter a is positive and that the involved state equation satisfies suitable assumptions. We explain and analyze how certain mapping properties of the operators appearing in the OS can be employed to define efficient preconditioners for finite element (FE) approximations of such systems. The key feature of the scheme is that the number of iterations needed to solve the preconditioned problem by the minimal residual method is bounded independently of the mesh parameter h, used in the FE discretization, and increases only moderately as alpha -> 0. More specifically, if the stopping criterion for the iteration process is defined in terms of the associated energy norm, then the number of iterations required (in the severely ill-posed case) cannot grow faster than O((ln(alpha))(2)). Our analysis is based on a careful study of the operators involved, which yields the distribution of the eigenvalues of the preconditioned OS. Finally, the theoretical results are illuminated by a number of numerical experiments addressing both a model problem studied by Borzi, Kunisch, and Kwak [SIAM J. Control Optim., 41 (2003), pp. 1477-1497] and an inverse problem arising in connection with electrocardiography [Nielsen, Cai, and Lysaker, Math. Biosci., 210 (2007), pp. 523-553].