This research develops an original algorithm for rich portfolio optimization (ARPO), considering more realistic constraints than those usually analyzed in the literature. Using a matheuristic framework that combines an iterated local search metaheuristic with quadratic programming, ARPO efficiently deals with complex variants of the mean-variance portfolio optimization problem, including the well-known cardinality and quantity constraints. ARPO proceeds in two steps. First, a feasible initial solution is constructed by allocating portfolio weights according to the individual return rate. Secondly, an iterated local search framework, which makes use of quadratic programming, gradually improves the initial solution throughout an iterative combination of a perturbation stage and a local search stage. According to the experimental results obtained, ARPO is very competitive when compared against existing state-of-the-art approaches, both in terms of the quality of the best solution generated as well as in terms of the computational times required to obtain it. Furthermore, we also show that our algorithm can be used to solve variants of the portfolio optimization problem, in which inputs (individual asset returns, variances and covariances) feature a random component. Notably, the results are similar to the benchmark constrained efficient frontier with deterministic inputs, if variances and covariances of individual asset returns comprise a random component. Finally, a sensitivity analysis has been carried out to test the stability of our algorithm against small variations in the input data.