Analytic solution for non-flat filter design
If you are not concerned on the phase and just want to approximate a magnitude response, then your first option should be the frequency sampling method implemented in Matlab/Octave fir2() function.

You would provide the frequency grid and corresponding frequency response magnitude at those frequencies.

As you have also mentioned, least-squares approach is another alternative. Indeed by using suitable weights, you can distribute the error according to your priority cirteria.

Magnitude approximation based on LMS adaptive system identification is also a possible option.