Thanks for the link @stevengj
Long story short: My objective function f(x::Vector) is an error function of the form ||\textbf{J}_T - \textbf{J}_{exp}(x) ||, where \textbf{J}_T is a target matrix and \textbf{J}_{exp}(x) is a simulated matrix. The latter is calculated from the phonon modes (\vec b) and frequencies (\lambda) of certain 2D systems, i.e. J^{(m,n)}_{exp}(x)= h(\vec b(x), \lambda(x)). To obtain \vec b(x), \lambda(x) I need to diagonalize some Hessian matrix \textbf{A}(x) describing the potentials of my system. So unfortunately it is not a more simple eigenvalue problem. I am not familiar with adjoint methods, so I would need to see if they could be applicable to my problem.
(I skipped the physical (fun!) details and whole expression of f(x) of my problem to keep it brief)