Linear Ising penalty functions can be used instead of quadratic ones, as shown here https://doi.org/10.1103/PhysRevResearch.6.043241 and here https://doi.org/10.1088/1367-2630/ad7e4a
Their implementation in Python is available here https://collections.durham.ac.uk/files/r2fq977t82m
We should implement it as an option in toQUBO.