Feature request: Sparse matrix bipartite matching behaving like linear sum assignment

Hi all.

I own a small matching service business and am finishing up the beta product.

I use scipy.optimize.linear_sum_assignment and it does precisely what I need it to do, and very quickly. Optimize matches in a bipartite graph with edge weights we have calculated. However as we scale, I know memory will eventually become an issue with this algorithm. While the number of real values in the graph will grow linearly thanks to design choices, the data in the matrix will still grow exponentially. So I will eventually need a more memory efficient algorithm - I can worry about computational complexity when it fits in memory. Each row will be limited to 1000 values - an arbitrary number to manage memory scaling which we might decrease at some point. When we scale to 1m users, that's 99.9% sparsity and 8TB of RAM worth of 64bit floats.

Currently scipy.sparse.csgraph.min_weight_full_bipartite_matching behaves slightly different than linear_sum_assignment. I give a toy example further down. It throws a value error (intentionally) when a row or column is unmatched. It also selects two edges, even when selecting only one can reach a bigger total sum of matched edges. If I have missed a way to do this currently, what I've written below is completely moot and I'd very much like to know how!

These look like intentional design choices I'd like to opt out of by setting a default/phantom edge weight. Adding a parameter with a default "None" argument to min_weight_full_bipartite_matching as to not cause breaking changes, or a separate implementation would both be attractive options.

As for the implementation, I'm not confident whether LAPJV or augmenting path would be better for this kind of memory optimization or changed behavior. Augmenting path is probably way easier, but LAPJV is probably more optimal. All I can think of for LAPJV, is initializing an array for each dimension and have those values represent the default value for the row and column. Those values are then subtracted from when the row/col it represents is. A phantom edge at x, y could be represented by the addition of the default values for x and y - if I have understood the algorithm properly. Honestly the code at https://github.com/scipy/scipy/blob/main/scipy/sparse/csgraph/_matching.pyx#L525 is a bit difficult for me to comprehend without diving super deep.

Thanks for reading and happy new year!
Erik "Barbarossa" Axelsson




Toy example:
```
import numpy as np
from scipy.sparse import coo_array
from scipy.sparse.csgraph import min_weight_full_bipartite_matching
from scipy.optimize import linear_sum_assignment

cases = [0]
# Different sums due to selecting row 1, col 0 and, row 5 col 4 - sum 0.65
# Which excludes selecting row 5, col 0 and row 1, col 4 (no edge) - sum 0.8
# [[0.   0.52 0.   0.74 0.   0.22]
#  [0.54 0.   0.   0.   0.   0.  ]
#  [0.   0.   0.32 0.47 0.   0.  ]
#  [0.   1.48 0.   0.   0.   1.26]
#  [0.33 0.47 0.   0.   0.   0.  ]
#  [0.8  0.   0.   0.   0.11 0.  ]]

cases.append(1)
# Throws ValueError due to collision between col 0 and col 3
# with the only available values on the same row
# [[0.   0.   0.   0.   0.79 0.06]
#  [0.2  0.   0.67 0.87 0.   0.  ]
#  [0.   0.   0.   0.   0.   0.29]
#  [0.   0.   0.   0.   0.83 0.8 ]
#  [0.   0.41 0.09 0.   0.59 0.  ]
#  [0.   0.67 0.83 0.   0.72 0.  ]]

cases.append(5)
# Throws ValueError due to no edge on row 2. Works fine if the row is removed.
# [[0.   1.37 1.05 0.   0.   0.56]
#  [0.   0.62 0.   0.   0.96 0.  ]
#  [0.   0.   0.   0.   0.   0.  ]
#  [0.69 1.   0.   0.   0.   0.83]
#  [0.   0.84 0.09 0.   0.   0.82]
#  [0.   0.56 0.   0.83 0.   0.  ]]


for i in cases:
    print(f"\ncase: {i}")
    np.random.seed(i)
    
    # Example sparse matrix
    row_indices = np.random.randint(0, 6, 15)
    col_indices = np.random.randint(0, 6, 15)
    data = np.random.random(15).round(2)
    print(f"rows: {row_indices}")
    print(f"cols: {col_indices}")

    # Creating a sparse matrix in coo_array format
    matrix = coo_array((data, (row_indices, col_indices)))

    # Linear sum assignment (maximize)
    print("Dense matrix:")
    dense_matrix = matrix.toarray()
    print(dense_matrix)
    rows, cols = linear_sum_assignment(dense_matrix, maximize=True)

    print("Optimal matching using linear_sum_assignment (max):")
    for row, col in zip(rows, cols):
        print(f"Row {row} matches with Column {col}: {dense_matrix[row, col]:.3}")
    
    lsa_sum = dense_matrix[rows, cols].sum()
    print(f"linear_sum_assignment sum: {lsa_sum}")


    # Max weight full bipartite matching
    rows, cols = min_weight_full_bipartite_matching(matrix, maximize=True)
    print("Optimal matching using min_weight_full_bipartite_matching (max):")
    for row, col in zip(rows, cols):
        print(f"Row {row} matches with Column {col}: {dense_matrix[row, col]:.3}")

    fbm_sum = dense_matrix[rows, cols].sum()
    print(f"min_weight_full_bipartite_matching sum: {fbm_sum}")
```
_______________________________________________
SciPy-Dev mailing list -- [email protected]
To unsubscribe send an email to [email protected]
https://mail.python.org/mailman3/lists/scipy-dev.python.org/
Member address: [email protected]