Parse large instances with 15K or 30K nodes

[N-Wouda]

The large instances of Arnold, Gendreau and Sörensen (2017) in CVRPLIB currently hang on the edge weight calculations. These instances are very large, ranging between 6-30K nodes. Let's see if we can support such large instances as well.

[leonlan]

I performed a a quick benchmark with the current distance computation + four alternatives:

Script

import numpy as np
import timeit
from tabulate import tabulate

import numba


def current(coords: np.ndarray) -> np.ndarray:
    diff = coords[:, np.newaxis, :] - coords
    square_diff = diff**2
    square_dist = np.sum(square_diff, axis=-1)
    return np.sqrt(square_dist)


def faster_current(pts: np.ndarray) -> np.ndarray:
    sq_sum = np.sum(pts**2, axis=1)
    sq_diff = np.add.outer(sq_sum, sq_sum) - 2 * np.dot(pts, pts.T)
    np.fill_diagonal(sq_diff, 0)  # avoids minor numerical issues
    return np.sqrt(sq_diff)


def numpy_linalg(X: np.ndarray) -> np.ndarray:
    return np.linalg.norm(X[:, np.newaxis] - X, axis=2)


@numba.jit(nopython=True, fastmath=True)
def variant_numba(a: np.ndarray) -> np.ndarray:
    n = a.shape[0]
    out = np.zeros((n, n))

    for i in range(n):
        for j in range(i):
            out[i][j] = euclidean(a[i], a[j])

    return out + out.T


@numba.jit(nopython=True, fastmath=True)
def euclidean(u, v):
    n = len(u)
    dist = 0
    for i in range(n):
        dist += abs(u[i] - v[i]) ** 2
    return dist ** (1 / 2)


def variant_scipy(pts: np.ndarray) -> np.ndarray:
    from scipy.spatial import distance_matrix

    return distance_matrix(pts, pts)


if __name__ == "__main__":
    np.random.seed(1)
    headers = ["# of coords", "Function", "Time (s)"]
    NUMBER = 3

    for num_coords in [5000, 10000, 15000]:
        coords = np.random.randn(num_coords, 2)

        funcs = [
            current,
            faster_current,
            numpy_linalg,
            variant_numba,
            variant_scipy,
        ]

        results = []
        for func in funcs:
            elapsed = timeit.timeit(
                f"{func.__name__}(coords)", globals=globals(), number=NUMBER
            )
            avg = round(elapsed / NUMBER, 2)
            results.append((num_coords, func.__name__, avg))

        print(tabulate(results, headers=headers))

    # current and numpy_linalg are too slow
    for num_coords in [20000, 25000, 30000]:
        coords = np.random.randn(num_coords, 2)

        funcs = [
            faster_current,
            variant_numba,
            variant_scipy,
        ]

        results = []
        for func in funcs:
            elapsed = timeit.timeit(
                f"{func.__name__}(coords)", globals=globals(), number=NUMBER
            )
            avg = round(elapsed / NUMBER, 2)
            results.append((num_coords, func.__name__, avg))

        print(tabulate(results, headers=headers))

Results

Results:

  # of coords  Function          Time (s)
-------------  --------------  ----------
         5000  current               0.4
         5000  faster_current        0.14
         5000  numpy_linalg          0.41
         5000  variant_numba         0.21
         5000  variant_scipy         0.33

  # of coords  Function          Time (s)
-------------  --------------  ----------
        10000  current               1.58
        10000  faster_current        0.44
        10000  numpy_linalg          1.58
        10000  variant_numba         0.36
        10000  variant_scipy         1.19

  # of coords  Function          Time (s)
-------------  --------------  ----------
        15000  current               4.29
        15000  faster_current        0.89
        15000  numpy_linalg          4.94
        15000  variant_numba         0.76
        15000  variant_scipy         2.62

  # of coords  Function          Time (s)
-------------  --------------  ----------
        20000  faster_current        1.54
        20000  variant_numba         1.39
        20000  variant_scipy         4.61

  # of coords  Function          Time (s)
-------------  --------------  ----------
        25000  faster_current        2.37
        25000  variant_numba         2.16
        25000  variant_scipy         7.2

  # of coords  Function          Time (s)
-------------  --------------  ----------
        30000  faster_current       66.12
        30000  variant_numba         6.8
        30000  variant_scipy        10.86

---

For instances up to 25K, the faster_current function is pretty fast. At 30K it becomes much slower than the numba or scipy - could this be due to memory issues? I'm not sure. Anyway, I think most users don't solve 30K instances anyway and the 60s-ish runtime is acceptable.
I also don't want vrplib to rely on heavy dependencies like numba or scipy, so sticking to the numpy version is preferred.

I'll implement the faster_current version soon.

---

Another option is to write some C++ code.

See this for setting up Poetry with a custom build.

[N-Wouda]

The scipy variant seems to scale very well and is faster than what we have now. Looking at the implementation here, the relevant bits are just a few lines of code. We could implement our calculation similarly?

[N-Wouda]

On my machine (Windows with an x86/i7 CPU from last year), the results look something like this:

  # of coords  Function          Time (s)
-------------  --------------  ----------
         5000  current               0.62
         5000  faster_current        0.28
         5000  numpy_linalg          0.63
         5000  variant_scipy         0.63

        10000  current               2.45
        10000  faster_current        1.14
        10000  numpy_linalg          2.48
        10000  variant_scipy         2.4

        15000  current               5.68
        15000  faster_current        2.71
        15000  numpy_linalg          6.16
        15000  variant_scipy         5.94

        20000  current               10.3
        20000  faster_current        5.71
        20000  variant_scipy        12.82

        25000  faster_current       10.2
        25000  variant_scipy        26.54

        30000  faster_current       27.49
        30000  variant_scipy        41.49

Here, faster_current is the clear winner, but scipy doesn't seem to do so well. Perhaps we should just implement faster_current and be done with it.

[leonlan]

Yeah agreed!