# Best Position for a Service Centre
**Difficulty:** HARD
[External](https://leetcode.com/problems/best-position-for-a-service-centre)
Canonical: https://scaleengineer.com/dsa/problems/best-position-for-a-service-centre
**Patterns:** [Math](https://scaleengineer.com/dsa/patterns/math), [Geometry](https://scaleengineer.com/dsa/patterns/geometry), [Randomized](https://scaleengineer.com/dsa/patterns/randomized)
**Data structures:** Array
**Companies:** [Citadel](https://scaleengineer.com/companies/citadel)
---
## Problem
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] = [xi, yi]` is the position of the `ith` 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 `[xcentre, ycentre]` such that the following formula is minimized:

![](https://assets.glich.co/dsa/best-position-for-a-service-centre/image0.jpg) 

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

**Example 1:**

![](https://assets.glich.co/dsa/best-position-for-a-service-centre/image1.jpg) 

**Input:** positions = [[0,1],[1,0],[1,2],[2,1]]
**Output:** 4.00000
**Explanation:** As shown, you can see that choosing [xcentre, ycentre] = [1, 1] will make the distance to each customer = 1, the sum of all distances is 4 which is the minimum possible we can achieve.

**Example 2:**

![](https://assets.glich.co/dsa/best-position-for-a-service-centre/image2.jpg) 

**Input:** positions = [[1,1],[3,3]]
**Output:** 2.82843
**Explanation:** The minimum possible sum of distances = sqrt(2) + sqrt(2) = 2.82843

**Constraints:**

* `1 <= positions.length <= 50`
* `positions[i].length == 2`
* `0 <= xi, yi <= 100`

# Approaches
## Hill Climbing with Shrinking Step Size
This approach treats the problem as finding the minimum point on a 2D surface. We start at an arbitrary point (e.g., the center of the coordinate space [50, 50]) and iteratively try to move to a neighboring point with a lower total distance. If we find a better neighbor, we move there. If all neighbors are worse, it means we are at a local minimum for the current step size. We then reduce the step size to explore the vicinity of the current point more finely. This process is repeated until the step size becomes smaller than the required precision.
**Time:** O(N * K) where N is the number of customer positions and K is the total number of steps taken. The number of steps K depends on the initial position, the step reduction factor, and the geometry of the points. It's generally larger than for more sophisticated methods. · **Space:** O(1) extra space, besides the storage for the input.
**Pros:** Relatively simple to understand and implement.; Guaranteed to find the global minimum because the objective function is convex.
**Cons:** Can be slow to converge if the starting point is far from the minimum or if the function's landscape forms a long, narrow valley.; It only explores four directions at each step, which is less efficient than moving along the true gradient.
### Explanation
The core idea is to perform a local search on a grid. We start with a large step size and a candidate center point. We explore the neighbors of this point (up, down, left, right). If any neighbor results in a smaller sum of Euclidean distances, we move our candidate center to that neighbor and repeat the exploration from the new point. If no neighbor offers improvement, we've found the best point for the current step size. We then shrink the step size (e.g., by half) and repeat the process. This allows us to first find the general area of the minimum and then zoom in to find a more precise location. The process terminates when the step size is negligibly small. Since the objective function is convex, this hill-climbing method is guaranteed to find the global minimum.

```java
class Solution {
    public double getMinDistSum(int[][] positions) {
        double current_x = 50.0;
        double current_y = 50.0;
        double step = 50.0;
        
        int[] dx = {0, 0, 1, -1};
        int[] dy = {1, -1, 0, 0};

        double minTotalDist = calculateTotalDistance(current_x, current_y, positions);

        while (step > 1e-7) {
            boolean movedInIteration = true;
            while(movedInIteration) {
                movedInIteration = false;
                for (int i = 0; i < 4; i++) {
                    double next_x = current_x + step * dx[i];
                    double next_y = current_y + step * dy[i];
                    double dist = calculateTotalDistance(next_x, next_y, positions);
                    if (dist < minTotalDist) {
                        minTotalDist = dist;
                        current_x = next_x;
                        current_y = next_y;
                        movedInIteration = true;
                    }
                }
            }
            step /= 2.0;
        }
        return minTotalDist;
    }

    private double calculateTotalDistance(double x, double y, int[][] positions) {
        double totalDist = 0;
        for (int[] pos : positions) {
            totalDist += Math.sqrt(Math.pow(x - pos[0], 2) + Math.pow(y - pos[1], 2));
        }
        return totalDist;
    }
}
```
### Algorithm
*   Initialize the service center position `(x, y)` to a starting point, for instance, the center of the bounding box `(50.0, 50.0)`.
*   Initialize a `step` size, for example, `50.0`.
*   Calculate the initial minimum distance `min_dist` for the starting `(x, y)`.
*   Define the directions to check: `(0, 1)`, `(0, -1)`, `(1, 0)`, `(-1, 0)`.
*   Loop while `step` is greater than a small epsilon (e.g., `1e-7`):
    *   Repeatedly search for a better point at the current step size by checking the four cardinal directions.
    *   If a neighbor `(new_x, new_y)` offers a smaller total distance, move the current point to this neighbor `(x, y) = (new_x, new_y)` and continue searching from there.
    *   If no neighbor offers improvement, break the inner search loop.
    *   Reduce the step size, e.g., `step /= 2`.
*   Return the final `min_dist`.

## Gradient Descent
Gradient descent is a first-order iterative optimization algorithm for finding the minimum of a function. The idea is to repeatedly take steps in the opposite direction of the gradient of the function at the current point, as this is the direction of steepest descent. For this problem, the function to minimize is the sum of Euclidean distances, which is a differentiable (almost everywhere) and convex function.
**Time:** O(N * K) where N is the number of positions and K is the number of iterations. K is typically smaller than for the hill-climbing approach. · **Space:** O(1) extra space.
**Pros:** Converges faster than the simple hill-climbing approach because it moves in the direction of the steepest descent.; A standard and well-understood optimization technique.
**Cons:** The choice of learning rate and decay schedule can be tricky and may require tuning for good performance.; Still an iterative approximation method.
### Explanation
First, we need the gradient of our objective function `f(x, y) = Σ sqrt((x - xᵢ)² + (y - yᵢ)²)`. The partial derivatives are:
`∂f/∂x = Σ (x - xᵢ) / dᵢ`
`∂f/∂y = Σ (y - yᵢ) / dᵢ`
where `dᵢ` is the distance from `(x, y)` to `(xᵢ, yᵢ)`.

The algorithm starts with an initial guess for the center `(x, y)`, often the centroid of all customer positions. Then, it iteratively updates the position by moving it a small amount in the direction opposite to the gradient. The size of this move is determined by a learning rate (or step size). A common strategy is to use a decaying learning rate, starting large to cover more ground and becoming smaller to fine-tune the position. The process continues for a fixed number of iterations or until the position converges. A special case to handle is when the candidate center `(x, y)` coincides with a customer position, making a distance `dᵢ` zero. This can be avoided by adding a small epsilon to the denominator.

```java
class Solution {
    public double getMinDistSum(int[][] positions) {
        double current_x = 0.0;
        double current_y = 0.0;
        int n = positions.length;
        for (int[] pos : positions) {
            current_x += pos[0];
            current_y += pos[1];
        }
        current_x /= n;
        current_y /= n;

        double step = 1.0;
        double epsilon = 1e-7;

        for (int i = 0; i < 10000; i++) {
            double grad_x = 0.0;
            double grad_y = 0.0;
            
            for (int[] pos : positions) {
                double dx = current_x - pos[0];
                double dy = current_y - pos[1];
                double dist = Math.sqrt(dx * dx + dy * dy);
                if (dist > 1e-9) { // Avoid division by zero
                    grad_x += dx / dist;
                    grad_y += dy / dist;
                }
            }

            double prev_x = current_x;
            double prev_y = current_y;
            current_x -= step * grad_x;
            current_y -= step * grad_y;
            
            if (Math.abs(current_x - prev_x) < epsilon && Math.abs(current_y - prev_y) < epsilon) {
                break;
            }
            
            step *= 0.999; // Decay step size
        }

        return calculateTotalDistance(current_x, current_y, positions);
    }

    private double calculateTotalDistance(double x, double y, int[][] positions) {
        double totalDist = 0;
        for (int[] pos : positions) {
            totalDist += Math.sqrt(Math.pow(x - pos[0], 2) + Math.pow(y - pos[1], 2));
        }
        return totalDist;
    }
}
```
### Algorithm
*   Initialize the service center position `(x, y)`. A good starting point is the centroid of all positions: `x = (Σ xᵢ) / n`, `y = (Σ yᵢ) / n`.
*   Initialize a learning rate `alpha` (e.g., `1.0`) and a decay factor (e.g., `0.999`).
*   Loop for a fixed number of iterations (e.g., 10000) or until the change in position is negligible:
    *   Initialize gradient components `grad_x = 0`, `grad_y = 0`.
    *   For each customer position `(xᵢ, yᵢ)`:
        *   Calculate the distance `dᵢ` from `(x, y)` to `(xᵢ, yᵢ)`.
        *   If `dᵢ` is close to zero, skip this point to avoid division by zero.
        *   Update the gradient: `grad_x += (x - xᵢ) / dᵢ`, `grad_y += (y - yᵢ) / dᵢ`.
    *   Update the position: `x -= alpha * grad_x`, `y -= alpha * grad_y`.
    *   Decay the learning rate: `alpha *= decay_factor`.
*   After the loop, calculate and return the total distance from the final `(x, y)` to all customer positions.

## Weiszfeld's Algorithm for Geometric Median
The problem of finding a point that minimizes the sum of Euclidean distances is a classic problem of finding the 'geometric median'. Weiszfeld's algorithm is a specialized iterative method designed specifically for this problem. It's a form of iteratively re-weighted least squares (IRWLS) and generally converges faster than standard gradient descent for this particular problem.
**Time:** O(N * K) where N is the number of positions and K is the number of iterations until convergence. K is typically very small, making this the most efficient approach among the three. · **Space:** O(1) extra space.
**Pros:** Specifically designed for the geometric median problem, making it very efficient.; Converges quickly in practice (linear convergence rate).; The update rule is simple and elegant.
**Cons:** The derivation is more mathematically involved than general-purpose methods.; Requires special handling for the case where an iterate coincides with a data point.
### Explanation
The algorithm is derived by setting the gradient of the sum-of-distances function to zero. This leads to an iterative update formula for the candidate center `(x, y)`:
`x_next = (Σ xᵢ / dᵢ) / (Σ 1 / dᵢ)`
`y_next = (Σ yᵢ / dᵢ) / (Σ 1 / dᵢ)`
where `dᵢ` is the distance from the current candidate center `(x_current, y_current)` to the customer point `(xᵢ, yᵢ)`.

Each term `1/dᵢ` can be seen as a 'weight' for the point `(xᵢ, yᵢ)`. Points closer to the current center are given higher weight. The next candidate center is the weighted average of all customer positions.

We start with an initial guess (e.g., the centroid) and repeatedly apply this update rule. The sequence of points generated by this process converges to the geometric median. The iteration stops when the change between successive points is smaller than a desired precision. A special case arises if an iterate lands on one of the customer points, which would cause division by zero. In this situation, we can check if that point is the true median. A simpler programmatic solution is to add a tiny epsilon to the distance `dᵢ` to prevent division by zero.

```java
class Solution {
    public double getMinDistSum(int[][] positions) {
        double current_x = 50.0;
        double current_y = 50.0;
        
        double prev_x, prev_y;
        double epsilon = 1e-7;

        do {
            prev_x = current_x;
            prev_y = current_y;

            double numerator_x = 0.0;
            double numerator_y = 0.0;
            double denominator = 0.0;

            for (int[] pos : positions) {
                double dist = Math.sqrt(Math.pow(current_x - pos[0], 2) + Math.pow(current_y - pos[1], 2));
                
                if (dist < 1e-9) {
                    continue;
                }
                
                numerator_x += pos[0] / dist;
                numerator_y += pos[1] / dist;
                denominator += 1 / dist;
            }

            if (denominator == 0) break;

            current_x = numerator_x / denominator;
            current_y = numerator_y / denominator;

        } while (Math.sqrt(Math.pow(current_x - prev_x, 2) + Math.pow(current_y - prev_y, 2)) > epsilon);

        return calculateTotalDistance(current_x, current_y, positions);
    }

    private double calculateTotalDistance(double x, double y, int[][] positions) {
        double totalDist = 0;
        for (int[] pos : positions) {
            totalDist += Math.sqrt(Math.pow(x - pos[0], 2) + Math.pow(y - pos[1], 2));
        }
        return totalDist;
    }
}
```
### Algorithm
*   Initialize the service center position `(x, y)` to the centroid of the customer positions.
*   Loop until the position converges:
    *   Store the current position: `prev_x = x`, `prev_y = y`.
    *   Initialize sums for the numerator and denominator of the update formula: `numerator_x = 0`, `numerator_y = 0`, `denominator = 0`.
    *   For each customer position `(xᵢ, yᵢ)`:
        *   Calculate the distance `dᵢ` from `(x, y)` to `(xᵢ, yᵢ)`.
        *   If `dᵢ` is very small, handle this special case (e.g., by stopping, as the point is likely the median).
        *   Calculate the weight `wᵢ = 1 / dᵢ`.
        *   Update the sums: `numerator_x += xᵢ * wᵢ`, `numerator_y += yᵢ * wᵢ`, `denominator += wᵢ`.
    *   Update the position: `x = numerator_x / denominator`, `y = numerator_y / denominator`.
    *   Check for convergence: if the distance between `(x, y)` and `(prev_x, prev_y)` is less than a small epsilon, break the loop.
*   Return the total distance from the final `(x, y)` to all customer positions.

# Solutions
### Java

```java
class Solution {
public
  double getMinDistSum(int[][] positions) {
    int n = positions.length;
    double x = 0, y = 0;
    for (int[] p : positions) {
      x += p[0];
      y += p[1];
    }
    x /= n;
    y /= n;
    double decay = 0.999;
    double eps = 1 e - 6;
    double alpha = 0.5;
    while (true) {
      double gradX = 0, gradY = 0;
      double dist = 0;
      for (int[] p : positions) {
        double a = x - p[0], b = y - p[1];
        double c = Math.sqrt(a * a + b * b);
        gradX += a / (c + 1 e - 8);
        gradY += b / (c + 1 e - 8);
        dist += c;
      }
      double dx = gradX * alpha, dy = gradY * alpha;
      if (Math.abs(dx) <= eps && Math.abs(dy) <= eps) {
        return dist;
      }
      x -= dx;
      y -= dy;
      alpha *= decay;
    }
  }
}

```

### CPP

```cpp
class Solution {
public:
  double getMinDistSum(vector<vector<int>> &positions) {
    int n = positions.size();
    double x = 0, y = 0;
    for (auto &p : positions) {
      x += p[0];
      y += p[1];
    }
    x /= n;
    y /= n;
    double decay = 0.999;
    double eps = 1e-6;
    double alpha = 0.5;
    while (true) {
      double gradX = 0, gradY = 0;
      double dist = 0;
      for (auto &p : positions) {
        double a = x - p[0], b = y - p[1];
        double c = sqrt(a * a + b * b);
        gradX += a / (c + 1e-8);
        gradY += b / (c + 1e-8);
        dist += c;
      }
      double dx = gradX * alpha, dy = gradY * alpha;
      if (abs(dx) <= eps && abs(dy) <= eps) {
        return dist;
      }
      x -= dx;
      y -= dy;
      alpha *= decay;
    }
  }
};

```

### Python

```python
class Solution:
    def getMinDistSum(self, positions: List[List[int]]) -> float: n = len(positions) x = y = 0 for x1, y1 in positions: x += x1 y += y1 x, y = x / n, y / n decay = 0.999 eps = 1e-6 alpha = 0.5 while 1: grad_x = grad_y = 0 dist = 0 for x1, y1 in positions: a = x - x1 b = y - y1 c = sqrt(a * a + b * b) grad_x += a / (c + 1e-8) grad_y += b / (c + 1e-8) dist += c dx = grad_x * alpha dy = grad_y * alpha x -= dx y -= dy alpha *= decay if abs(dx) <= eps and abs(dy) <= eps: return dist

```
