title: Best Position for a Service Centre
slug: best-position-for-a-service-centre
difficulty: hard
leetcode_id: 1515
leetcode_url: https://leetcode.com/problems/best-position-for-a-service-centre/
categories:
  - arrays
  - math
patterns:
  - slug: greedy
    is_optimal: true

function_signature: "def get_min_dist_sum(positions: list[list[int]]) -> float:"

test_cases:
  visible:
    - input: { positions: [[0, 1], [1, 0], [1, 2], [2, 1]] }
      expected: 4.00000
    - input: { positions: [[1, 1], [3, 3]] }
      expected: 2.82843
  hidden:
    - input: { positions: [[0, 0]] }
      expected: 0.00000
    - input: { positions: [[0, 0], [2, 0]] }
      expected: 2.00000
    - input: { positions: [[0, 0], [0, 1], [1, 0], [1, 1]] }
      expected: 2.82843
    - input: { positions: [[1, 1], [1, 1], [1, 1]] }
      expected: 0.00000
    - input: { positions: [[0, 0], [50, 50], [100, 100]] }
      expected: 141.42136
    - input: { positions: [[1, 2], [3, 4], [5, 6], [7, 8]] }
      expected: 11.31371

description: |
  A delivery company wants to build a new service center in a new city. The company knows the positions of all the customers in this city on a 2D-Map and wants to build the new center in a position such that **the sum of the Euclidean distances to all customers is minimum**.

  Given an array `positions` where `positions[i] = [x_i, y_i]` is the position of the i<sup>th</sup> customer on the map, return *the minimum sum of the Euclidean distances* to all customers.

  In other words, you need to choose the position of the service center `[x_centre, y_centre]` such that the following formula is minimised:

  `∑ sqrt((x_centre - x_i)² + (y_centre - y_i)²)`

  Answers within `10^-5` of the actual value will be accepted.

constraints: |
  - `1 <= positions.length <= 50`
  - `positions[i].length == 2`
  - `0 <= x_i, y_i <= 100`

examples:
  - input: "positions = [[0,1],[1,0],[1,2],[2,1]]"
    output: "4.00000"
    explanation: "Choosing [x_centre, y_centre] = [1, 1] makes the distance to each customer = 1, so the sum of all distances is 4, which is the minimum possible."
  - input: "positions = [[1,1],[3,3]]"
    output: "2.82843"
    explanation: "The minimum possible sum of distances = sqrt(2) + sqrt(2) = 2.82843. The optimal centre lies on the line between the two points."

explanation:
  intuition: |
    Imagine you're dropping a pin on a map, and rubber bands connect it to every customer location. The pin naturally settles at the position where the total tension (sum of distances) is minimised. This point is called the **geometric median** or **Fermat point**.

    Unlike the *centroid* (arithmetic mean of coordinates), which minimises the sum of **squared** distances, the geometric median minimises the sum of **actual** distances. There's no closed-form formula for this — we must use iterative optimisation.

    Think of it like this: start somewhere reasonable (the centroid is a good initial guess), then repeatedly nudge the point in the direction that reduces the total distance. Each step, we compute the gradient of our objective function and move against it.

    The key insight is that while the objective function (sum of Euclidean distances) isn't simple to solve analytically, it's **convex** — meaning any local minimum is the global minimum. This guarantees that gradient-based methods will converge to the optimal solution.

  approach: |
    We solve this using **Gradient Descent** (or alternatively, Weiszfeld's algorithm):

    **Step 1: Define the objective function**

    - `f(x, y)` = sum of Euclidean distances from `(x, y)` to all customer positions
    - `f(x, y) = ∑ sqrt((x - x_i)² + (y - y_i)²)`

    &nbsp;

    **Step 2: Initialise the starting point**

    - Compute the centroid of all positions as a reasonable starting point
    - `x = mean(x_i)`, `y = mean(y_i)`

    &nbsp;

    **Step 3: Compute the gradient**

    - The partial derivatives tell us which direction increases the objective:
    - `∂f/∂x = ∑ (x - x_i) / dist_i` where `dist_i = sqrt((x - x_i)² + (y - y_i)²)`
    - `∂f/∂y = ∑ (y - y_i) / dist_i`
    - We move in the *opposite* direction to decrease the total distance

    &nbsp;

    **Step 4: Iteratively update the position**

    - Use a learning rate `α` that decreases over time for stable convergence
    - `x_new = x - α * ∂f/∂x`
    - `y_new = y - α * ∂f/∂y`
    - Repeat until the change is smaller than the required precision (`10^-7`)

    &nbsp;

    **Step 5: Return the minimum sum**

    - Compute and return `f(x_final, y_final)`

    &nbsp;

    The algorithm is guaranteed to converge because the sum of Euclidean distances is a convex function, meaning there's exactly one global minimum.

  common_pitfalls:
    - title: Confusing Geometric Median with Centroid
      description: |
        The centroid (arithmetic mean of x and y coordinates separately) minimises the sum of **squared** distances, not actual distances. For this problem, you cannot simply compute:
        - `x_centre = mean(x_i)`, `y_centre = mean(y_i)`

        While the centroid is a good starting point for gradient descent, it's rarely the optimal answer. For example, with points at `[[0,0], [0,3], [4,0]]`, the centroid is `(4/3, 1)`, but the geometric median is different.
      wrong_approach: "Return the centroid directly"
      correct_approach: "Use centroid as starting point, then optimise with gradient descent"

    - title: Fixed Learning Rate Causing Oscillation
      description: |
        Using a constant learning rate can cause the algorithm to oscillate around the minimum without converging precisely enough to meet the `10^-5` tolerance.

        The learning rate should decrease over iterations. A common strategy is to start with a larger rate (e.g., `1.0`) and multiply by a decay factor (e.g., `0.999`) each iteration, or halve it when no improvement is made.
      wrong_approach: "Use fixed learning rate α = 0.01"
      correct_approach: "Use adaptive or decaying learning rate"

    - title: Division by Zero at Customer Positions
      description: |
        When computing the gradient, if the current position exactly matches a customer position, `dist_i = 0` causes division by zero.

        Handle this by adding a small epsilon (e.g., `10^-10`) to the distance, or skip that customer's contribution when the distance is effectively zero.
      wrong_approach: "Compute gradient without checking for zero distance"
      correct_approach: "Add epsilon to distance or handle zero-distance case"

  key_takeaways:
    - "**Geometric median vs centroid**: The centroid minimises squared distances; the geometric median minimises actual distances — they require different algorithms"
    - "**Gradient descent for continuous optimisation**: When no closed-form solution exists, iterative methods like gradient descent can find optimal solutions"
    - "**Convexity guarantees convergence**: The sum of Euclidean distances is convex, so gradient descent will find the global minimum"
    - "**Adaptive learning rate**: For numerical precision, decrease the step size as you approach the solution"

  time_complexity: "O(n × k) where n is the number of positions and k is the number of iterations. Typically k is in the hundreds to achieve the required precision."
  space_complexity: "O(1). We only store the current position and gradient values, regardless of input size."

solutions:
  - approach_name: Gradient Descent
    is_optimal: true
    code: |
      import math

      def get_min_dist_sum(positions: list[list[int]]) -> float:
          def total_distance(x: float, y: float) -> float:
              """Calculate sum of Euclidean distances from (x, y) to all positions."""
              return sum(
                  math.sqrt((x - px) ** 2 + (y - py) ** 2)
                  for px, py in positions
              )

          def compute_gradient(x: float, y: float) -> tuple[float, float]:
              """Compute partial derivatives of the distance sum."""
              grad_x, grad_y = 0.0, 0.0
              for px, py in positions:
                  dist = math.sqrt((x - px) ** 2 + (y - py) ** 2)
                  if dist > 1e-10:  # Avoid division by zero
                      grad_x += (x - px) / dist
                      grad_y += (y - py) / dist
              return grad_x, grad_y

          # Start at centroid (good initial guess)
          x = sum(p[0] for p in positions) / len(positions)
          y = sum(p[1] for p in positions) / len(positions)

          # Gradient descent with decaying learning rate
          learning_rate = 1.0
          decay = 0.999
          epsilon = 1e-7

          for _ in range(100000):
              grad_x, grad_y = compute_gradient(x, y)

              # Update position (move against gradient)
              new_x = x - learning_rate * grad_x
              new_y = y - learning_rate * grad_y

              # Check for convergence
              if abs(new_x - x) < epsilon and abs(new_y - y) < epsilon:
                  break

              x, y = new_x, new_y
              learning_rate *= decay

          return total_distance(x, y)
    explanation: |
      **Time Complexity:** O(n × k) — Each iteration computes distances to all n points, running for up to k iterations.

      **Space Complexity:** O(1) — Only stores current position and gradient values.

      We start at the centroid and iteratively move in the direction that decreases total distance. The decaying learning rate ensures we take smaller steps as we approach the optimum, achieving the required precision.

  - approach_name: Weiszfeld's Algorithm
    is_optimal: true
    code: |
      import math

      def get_min_dist_sum(positions: list[list[int]]) -> float:
          def total_distance(x: float, y: float) -> float:
              """Calculate sum of Euclidean distances from (x, y) to all positions."""
              return sum(
                  math.sqrt((x - px) ** 2 + (y - py) ** 2)
                  for px, py in positions
              )

          # Start at centroid
          x = sum(p[0] for p in positions) / len(positions)
          y = sum(p[1] for p in positions) / len(positions)

          epsilon = 1e-7

          for _ in range(1000):
              # Compute weights (inverse distances)
              weights = []
              for px, py in positions:
                  dist = math.sqrt((x - px) ** 2 + (y - py) ** 2)
                  # Handle case when we're exactly on a point
                  weights.append(1.0 / max(dist, epsilon))

              # Weighted average gives new position
              total_weight = sum(weights)
              new_x = sum(w * p[0] for w, p in zip(weights, positions)) / total_weight
              new_y = sum(w * p[1] for w, p in zip(weights, positions)) / total_weight

              # Check convergence
              if abs(new_x - x) < epsilon and abs(new_y - y) < epsilon:
                  break

              x, y = new_x, new_y

          return total_distance(x, y)
    explanation: |
      **Time Complexity:** O(n × k) — Similar to gradient descent, each iteration visits all n points.

      **Space Complexity:** O(n) — Stores weights for all positions.

      Weiszfeld's algorithm is a specialised iterative method for the geometric median. Each iteration computes a weighted average where weights are inverse distances. Points closer to the current estimate have higher influence, naturally pulling the estimate toward the optimal position.

  - approach_name: Simulated Annealing
    is_optimal: false
    code: |
      import math
      import random

      def get_min_dist_sum(positions: list[list[int]]) -> float:
          def total_distance(x: float, y: float) -> float:
              return sum(
                  math.sqrt((x - px) ** 2 + (y - py) ** 2)
                  for px, py in positions
              )

          # Start at centroid
          x = sum(p[0] for p in positions) / len(positions)
          y = sum(p[1] for p in positions) / len(positions)
          best_dist = total_distance(x, y)

          # Simulated annealing
          temp = 100.0
          cooling_rate = 0.9999

          for _ in range(100000):
              # Random step proportional to temperature
              new_x = x + (random.random() - 0.5) * temp
              new_y = y + (random.random() - 0.5) * temp

              new_dist = total_distance(new_x, new_y)

              # Accept if better, or probabilistically if worse
              if new_dist < best_dist:
                  x, y = new_x, new_y
                  best_dist = new_dist
              elif random.random() < math.exp((best_dist - new_dist) / temp):
                  x, y = new_x, new_y

              temp *= cooling_rate

          return best_dist
    explanation: |
      **Time Complexity:** O(n × k) — Similar iteration count to gradient descent.

      **Space Complexity:** O(1) — Only stores current and best positions.

      Simulated annealing is a probabilistic approach that can escape local minima. While the geometric median problem is convex (no local minima), this approach still works and demonstrates an alternative optimisation technique. It's less deterministic than gradient descent but more robust for non-convex problems.
