Minimize (firstA_max - firstA_min) + (secondB_max - secondB_min)

The following approach is an attempt to escape the n^2, using an argmin list for the second element of the tuples (lets say the y-part). Where the points are sorted regarding x.

One Observation is that there is an optimum solution where A includes index argmin[0] or argmin[n-1] or both.

1. in get_best_interval_min_max we focus once on including argmin[0] and the next smallest element on y and so one. The we do the same from the max element.

We get two dictionaries {(i,j):(profit, idx)}, telling us how much we gain in y when including points[i:j+1] in A, towards min or max on y. idx is the idx in the argmin array.

2. calculate the objective for each dict assuming max/min or y is not in A.

3. combine the results of both dictionaries, : (i1,j1): (v1, idx1) and (i2,j2): (v2, idx2). result : j2 - i1 + max_y - min_y - v1 - v2.

Constraint: idx1 < idx2. Because the indices in the argmin array can not intersect, otherwise some profit in y might be counted twice.

On average the dictionaries (dmin,dmax) are smaller than n, but in the worst case when x and y correlate [(i,i) for i in range(n)] they are exactly n, and we do not win any time. Anyhow on random instances this approach is much faster. Maybe someone can improve upon this.

import numpy as np
from random import randrange
import time

def get_best_interval_min_max(points):# sorted input according to x dim 
    L = len(points)
    argmin_b = np.argsort([p[1] for p in points])
    b_min,b_max = points[argmin_b[0]][1], points[argmin_b[L-1]][1]
    
    arg = [argmin_b[0],argmin_b[0]]
    res_min = dict()
    for i in range(1,L):
        res_min[tuple(arg)] = points[argmin_b[i]][1] - points[argmin_b[0]][1],i # the profit in b towards min
        if arg[0] > argmin_b[i]: arg[0]=argmin_b[i]
        elif arg[1] < argmin_b[i]: arg[1]=argmin_b[i]
        
    arg = [argmin_b[L-1],argmin_b[L-1]]
    res_max = dict()
    for i in range(L-2,-1,-1):
        res_max[tuple(arg)] = points[argmin_b[L-1]][1]-points[argmin_b[i]][1],i # the profit in b towards max
        if arg[0]>argmin_b[i]: arg[0]=argmin_b[i]
        elif arg[1]<argmin_b[i]: arg[1]=argmin_b[i]
    # return the two dicts, difference along y,     
    return res_min, res_max, b_max-b_min

def argmin_algo(points):
    # return the objective value, sets A and B, and the interval for A in points. 
    points.sort()
    # get the profits for different intervals on the sorted array for max and min
    dmin, dmax, y_diff = get_best_interval_min_max(points)
    key = [None,None]
    res_min = 2e9
    # the best result when only the min/max b value is includes in A
    for d in [dmin,dmax]:
        for k,(v,i) in d.items():
            res = points[k[1]][0]-points[k[0]][0] + y_diff - v
            if res < res_min: 
                key = k
                res_min = res

    # combine the results for max and min. 
    for k1,(v1,i) in dmin.items():
        for k2,(v2,j) in dmax.items():
            if i > j: break # their argmin_b indices can not intersect!
            idx_l, idx_h = min(k1[0], k2[0]), max(k1[1],k2[1]) # get index low and idx hight for combination
            res = points[idx_h][0]-points[idx_l][0] -v1 -v2 + y_diff
            if res < res_min: 
                key = (idx_l, idx_h) # new merged interval
                res_min = res
    return res_min, points[key[0]:key[1]+1], points[:key[0]]+points[key[1]+1:], key

def quadratic_algorithm(points):
    points.sort()
    m, M, res = 1e9, -1e9, 2e9
    idx = (0,0)
    for l in range(len(points)):
        g, G = m, M 
        for r in range(len(points)-1,l-1,-1):
            if r-l+1 < len(points):
                res_n = points[r][0] - points[l][0] + G - g
                if res_n < res:
                    res = res_n
                    idx = (l,r)
            g = min(g, points[r][1])
            G = max(G, points[r][1])
        m = min(m, points[l][1])
        M = max(M, points[l][1])
    return res, points[idx[0]:idx[1]+1], points[:idx[0]]+points[idx[1]+1:], idx

# let's try it and compare running times to the quadratic_algorithm
# get some "random" points
c1=0
c2=0
for i in range(100):
    points = [(randrange(100), randrange(100)) for i in range(1,200)]
    points.sort() # sorted for x dimention
    s = time.time()
    r1 = argmin_algo(points)
    e1 = time.time()
    r2 = quadratic_algorithm(points)
    e2 = time.time()
    c1 += (e1-s)
    c2 += (e2-e1)
    if not r1[0] == r2[0]:
        print(r1,r2)
        raise Exception("Error, results are not equal")
print("time of argmin_algo", c1, "time of quadratic_algorithm",c2)

UPDATE: @Luka proved the algorithm described in this answer is not exact. But I will keep it here because it's a good performance heuristics and opens the way to many probabilistic methods.

I will describe a loglinear algorithm. I couldn't find a counter example. But I also couldn't find a proof :/

Let set A be ordered by first element and set B be ordered by second element. They are initially empty. Take floor(n/2) random points of your set of points and put in set A. Put the remaining points in set B. Define this as a partition.

Let's call a partition stable if you can't take an element of set A, put it in B and decrease the objective function and if you can't take an element of set B, put it in A and decrease the objective function. Otherwise, let's call the partition unstable.

For an unstable partition, the only moves that are interesting are the ones that take the first or the last element of A and move to B or take the first or the last element of B and move to A. So, we can find all interesting moves for a given unstable partition in O(1). If an interesting move decreases the objective function, do it. Go like that until the partition becomes stable. I conjecture that it takes at most O(n) moves for the partition to become stable. I also conjecture that at the moment the partition becomes stable, you will have a solution.