Tikhonov regularization is a popular method to approximate solutions of linear discrete ill-posed problems when the observed or measured data is contaminated by noise. Multi-parameter Tikhonov regularization may improve the quality of the computed approximate solutions. We propose a new iterative method for large-scale multi-parameter Tikhonov regularization with general regularization operators based on a multidirectional subspace expansion. This expansion may be combined with subspace truncation to avoid excessive growth of the search space. Furthermore, we introduce a simple and effective parameter selection strategy based on the discrepancy principle and related to perturbation results.