Quadratization problem is, given a system of ODEs with polynomial right-hand side, transform the system to a system with quadratic right-hand side by introducing new variables. Such transformations have been used, for example, as a preprocessing step by model order reduction methods and for transforming chemical reaction networks. We present an algorithm that, given a system of polynomial ODEs, finds a transformation into a quadratic ODE system by introducing new variables which are monomials in the original variables. The algorithm is guaranteed to produce an optimal transformation of this form (that is, the number of new variables is as small as possible), and it is the first algorithm with such a guarantee we are aware of. Its performance compares favorably with the existing software, and it is capable to tackle problems that were out of reach before.