At first, I found that there's an optimal solution where each square's left and upper edge must touch a red point. Otherwise we can simply shift those squares.
This lead to a straightforward bit-mask DP. Suppose we sort the points along the x axis then the y axis. Then for each mask(bit 1 represents a point already been covered), we find the lowest 0 bit(corresponding to the left most uncovered point thus is the left edge of the square), iterate through all points while using those to the upper right of it as possible upper edges, iterate through all possible squares and calculate the new mask to do the DP. This results in O((2^n)*(n^2)*m), which barely passed under the time limit by some careful pruning(ignore points already covered when iterating through points).
As I looked at others' solutions in the practice room, I found a clever technique to optimize. Instead of explicitly determine the position of the squares, we can first simply calculate all the masks using exactly one square. Then using the techniques to iterate through subsets, we can set dp[mask] to min(dp[i]+dp[maks^i], mask|i==mask). Inductively, if two subsets can both be covered by an arbitrary number of squares, so is the whole set. This results in O(3^n).