Maximize A[i]*B[i] + A[i]*B[j] + A[j]*B[j], i != j, given two unordered lists of positive integers

Question

Can you help me with the algorithm: Given 2 arrays of equal size a[] and b[] with integer numbers greater or equal to 1.

Find non-equal indices i and j (i != j) such that the value -

max(a[i]*b[i] + a[i] * b[j] + a[j] * b[j], a[i]*b[i] + a[j] * b[i] + a[j] * b[j])

will be maximum.

Example:

a = [1, 9, 6, 6] and b = [9, 1, 6, 6].

Maximum will be at i = 2 and j = 3 (zero-based indexing): a[2]*b[2] + a[2]*b[3] + a[3] * b[3] = 6*6+6*6+6*6 = 108

Is there a way to find i and j in less then quadratic time? And the same question for finding objective function value in less then quadratic time?

Thank you!

Answer 1

Here's my attempt at implementing David Eisenstat's idea. I think the i != j restriction made this more complicated but either way, I'd welcome suggestions on improving the code. There's a test at the end against brute force.

The construction of the upper envelope for the lines A[i]*x + A[i]*B[i] relies on Andrew's monotone chain convex hull algorithm applied to the dual points transformed from the lines.

Python code:

# Upper envelope of lines in the plane

from fractions import Fraction
import collections

def get_dual_point(line):
  return (line[0], -line[1])

def get_primal_point(dual_line):
  return (dual_line[0], -dual_line[1])

def get_line_from_two_points(p1, p2):
  if p1[0] == p2[0]:
    return (float('inf'), 0)

  m = Fraction(p1[1] - p2[1], p1[0] - p2[0])
  b = -m * p1[0] + p1[1]

  return (m, b)

def cross_product(o, a, b):
  return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0])


# lower_hull has the structure [(dual_point, original_line, leftmost_x, index)]
def get_segment_idx(lower_hull, x):
  lo = 0
  hi = len(lower_hull) - 1

  # Find the index of the first
  # x coordinate greater than x
  while lo < hi:
    mid = lo + (hi - lo + 1) // 2

    if lower_hull[mid][2] <= x:
      lo = mid
    else:
      hi = mid - 1

  return lo


# Assumes we add points in order of increasing x-coordinates
def add_right_to_lower_hull(lower_hull, point):
  while len(lower_hull) > 0 and lower_hull[-1][0][0] == point[0][0] and lower_hull[-1][0][1] > point[0][1]:
    lower_hull.pop()
  while len(lower_hull) > 1 and cross_product(lower_hull[-2][0], lower_hull[-1][0], point[0]) <= 0:
    lower_hull.pop()
  if not lower_hull or lower_hull[-1][0][0] != point[0][0]:
    lower_hull.append(point)
  # Each segment of the lower hull
  # in the dual plane is a line intersection
  # in the primal plane.
  if len(lower_hull) == 1:
    lower_hull[0][2] = -float('inf')
  else:
    line = get_line_from_two_points(lower_hull[-1][0], lower_hull[-2][0])
    lower_hull[-1][2] = get_primal_point(line)[0]

  return lower_hull


# Assumes we add points in order of decreasing x-coordinates
def add_left_to_lower_hull(lower_hull, point):
  while len(lower_hull) > 0 and lower_hull[0][0][0] == point[0][0] and lower_hull[0][0][1] > point[0][1]:
    lower_hull.popleft()
  while len(lower_hull) > 1 and cross_product(lower_hull[1][0], lower_hull[0][0], point[0]) >= 0:
    lower_hull.popleft()
  if not lower_hull or lower_hull[0][0][0] != point[0][0]:
    lower_hull.appendleft(point)
  # Each segment of the lower hull
  # in the dual plane is a line intersection
  # in the primal plane.
  if len(lower_hull) == 1:
    lower_hull[0][2] = -float('inf')
  else:
    line = get_line_from_two_points(lower_hull[0][0], lower_hull[1][0])
    lower_hull[1][2] = get_primal_point(line)[0]

  return lower_hull


# Maximise A[i] * B[i] + A[i] * B[j] + A[j] * B[j]
def f(A, B):
  debug = False

  if debug:
    print("A: %s" % A)
    print("B: %s" % B)

  best = -float('inf')
  best_idxs = ()

  indexed_lines = [((A[i], A[i] * B[i]), i) for i in range(len(A))]

  # Convert to points in the dual plane
  # [dual_point, original_line, leftmost x coordinate added later, original index]
  dual_points = [[get_dual_point(line), line, None, i] for line, i in indexed_lines]

  # Sort points by x coordinate ascending
  sorted_points = sorted(dual_points, key=lambda x: x[0][0])

  if debug:
    print("sorted points")
    print(sorted_points)

  # Build lower hull, left to right
  lower_hull = []

  add_right_to_lower_hull(lower_hull, sorted_points[0])

  for i in range (1, len(sorted_points)):
    # Query the point before inserting it
    # because of the stipulation that i != j
    idx = sorted_points[i][3]
    segment_idx = get_segment_idx(lower_hull, B[idx])
    m, b = lower_hull[segment_idx][1]
    j = lower_hull[segment_idx][3]
    candidate = m * B[idx] + b + A[idx] * B[idx]

    if debug:
      print("segment: %s, idx: %s, B[idx]: %s" % (segment_idx, idx, B[idx]))

    if candidate > best:
      best = candidate
      best_idxs = (idx, j)

    add_right_to_lower_hull(lower_hull, sorted_points[i])
  
  if debug:
    print("lower hull")
    print(lower_hull)

  # Build lower hull, right to left
  lower_hull = collections.deque()

  lower_hull.append(sorted_points[len(sorted_points) - 1])

  for i in range (len(sorted_points) - 2, -1, -1):
    # Query the point before inserting it
    # because of the stipulation that i != j
    idx = sorted_points[i][3]
    segment_idx = get_segment_idx(lower_hull, B[idx])
    m, b = lower_hull[segment_idx][1]
    j = lower_hull[segment_idx][3]
    candidate = m * B[idx] + b + A[idx] * B[idx]

    if debug:
      print("segment: %s, idx: %s, B[idx]: %s" % (segment_idx, idx, B[idx]))

    if candidate > best:
      best = candidate
      best_idxs = (idx, j)

    add_left_to_lower_hull(lower_hull, sorted_points[i])

  if debug:
    print("lower hull")
    print(lower_hull)

  return best, best_idxs


#A = [1, 9, 6, 6]
#B = [9, 1, 6, 6]

#print("")
#print(f(A, B))


# Test

import random

def brute_force(A, B):
  best = -float('inf')
  best_idxs = ()

  for i in range(len(A)):
    for j in range(len(B)):
      if i != j:
        candidate = A[i] * B[i] + A[i] * B[j] + A[j] * B[j]
        if candidate > best:
          best = candidate
          best_idxs = (i, j)

  return best, best_idxs


num_tests = 500
n = 20
m = 1000

for _ in range(num_tests):
  A = [random.randint(1, m) for i in range(n)]
  B = [random.randint(1, m) for i in range(n)]

  _f = f(A, B)
  _brute = brute_force(A, B)

  if _f[0] != _brute[0]:
    print("Mismatch")
    print(A)
    print(B)
    print(_f, _brute)

print("Done testing.")

Answer 2

Yes, there’s an O(n log n)-time algorithm.

If we look at the function f(x) = max {a[i] b[i] + a[i] x | i}, it traces out a lower (convex) hull. It basically suffices to evaluate f(b[j]) + a[j] b[j] for every j, with the one sticking point that we need i ≠ j.

To fix that problem, we turn to an incremental hull algorithm with amortized O(log n)-time insertions. The canonical algorithm keeps the hull in sorted order, allowing for O(log n)-time queries. Starting with an empty hull, for each index i, we query i and then insert i. Do this once in forward order and once in reverse order, and take the max.