Learning algorithms often rely on effective hyperparameter tuning, particularly in penalized optimization problems where constraints aid in extracting interpretable, low-dimensional data representations. This work explores bi-level and multi-level optimization strategies for Nonnegative Matrix Factorization (NMF), a powerful tool for uncovering latent structures while enforcing nonnegativity. We introduce a family of novel algorithms that jointly address factorization and hyperparameter tuning across different divergence measures and domains. Our first contribution, AltBi, formulates Kullback-Leibler NMF as a bi-level optimization problem where penalty hyperparameters are adaptively updated through alternating minimization steps. This algorithm demonstrates enhanced interpretability and reconstruction performance on signal-domain datasets, as shown in [1]. Building on this, AltBi-J integrates a diversity-based penalty into the Frobenius norm framework, yielding solutions with significantly higher sparsity and better components disentanglement, particularly suited for structured data [2]. Extending our approach to the spectral domain, we propose SHINBO, a bi-level optimization algorithm for Itakura-Saito NMF. SHINBO introduces adaptive, row-wise penalty tuning that effectively isolates periodic components from noisy signals, with effective application to fault detection in mechanical systems [3]. Motivated by learning applications involving multiple constraints, we generalize bi-level methodology into a multi-level optimization paradigm capable of handling nested, nonconvex problems. This hierarchical perspective offers a robust theoretical and algorithmic foundation for advancing optimization-driven learning in structured and high-dimensional environments. This is a joint work with Flavia Esposito and Nicoletta Del Buono from the University of Bari Aldo Moro, Andersen Ang from the University of Southamptopn, and Rafal Zdunek from Politechnika Wrocławska.