Yes. First of all, your norm() call implies Frobenius: but any kind of least-squares is totally uncalled for in this problem since you just want integer comparison.
To put it lightly, the description of what this is actually supposed to do is troubled. For instance, rather than "labelling" I would call it "substitution". I hope the extended back-and-forth has prepared you for the kind of information reviewers need to know in future questions.
A better algorithm is definitely possible. First, frame the problem as a bipartite graph, where:
- the left subset nodes are unique values from the original
arr2
- the right subset nodes are all possible destination values from 0 through n-1
- the graph starts as being complete and unbalanced; the left cardinality will typically be smaller than the right cardinality
- the edge weights are non-negative integers representing the number of mismatches between the source and destination
The problem becomes: how do you prune this graph such that
- for every left-node, there remains exactly one edge (the left subset becomes 1-regular);
- for every right-node, there remains at most one edge;
- the graph becomes disconnected and acyclic; and
- the edges are chosen such that the total remaining weight is minimised.
This is a known problem - a variant of bipartite maximum matching called Bipartite Min-Cost Perfect Matching, also called the assignment problem. Crucially, the insight from interpreting this as a bipartite matching problem is that it can be solved using mixed-integer linear programming (MIP).
MIP problems are well-understood and have excellent software that already exists to solve them. I recommend glpk and its Python SWIG bindings in swiglpk.
In this strategy,
c, the cost vector, will be non-negative integers capturing your edge weightsz, the objective, is the total distancexs, the structural variables/columns, will be a flattened 2d matrix of boolean edge choicesxr, the auxiliary variables/rows, will be the node degree (sum of edge counts per node)A, the constraint matrix, will be constructed to sum up the node degrees over both axes
xs, when interpreted as a 2D matrix, must be constrained to have exactly one nonzero per (source) row, and between 0 and 1 nonzeros per (dest) column. In glpk these are expressed with fixed and double bounds, respectively.
It takes a little bit of code to express all of this, but it's quite worth it: when I run a problem of n=500 on my machine it completes in about two seconds.
Suggested
import random
from collections import defaultdict
from timeit import timeit
from typing import Optional, Iterator
import swiglpk as lp
class OptError(Exception):
pass
def make_map() -> tuple[tuple[int, list[int]], ...]:
"""
From arr1 and arr2, produce an ordered map in tuple-of-pairs format
where the key is the unique source arr2 value, and the values are all
(potentially non-unique) arr1 values to which a match will be attempted.
"""
targets = defaultdict(list)
for source, dest in zip(arr2, arr1):
targets[source].append(dest)
return tuple(targets.items())
def fill_sources() -> None:
"""
For each of our unique arr2 source values,
- Fill some the GLPK problem's "auxiliary rows" of count n_source. These
rows correspond to constraints on the arr2 source values.
- Name those rows and apply fixed bounds: each source must get exactly one
match.
For each combination of source and candidate destination,
- Reinterpret the source/destination matrix as a flat, horizontal vector
- Fill all of the GLPK problem's "structural columns" of count
n_source*n_dest.
- Name those columns and set the structural variable type as binary-value.
- Set the corresponding objective coefficient to the mismatch cost of
choosing this pair.
- Populate the constraint matrix to enforce that each source gets one dest,
and every dest gets up to one source
"""
row_indices = lp.intArray(3)
row_values = lp.doubleArray(3)
row_values[1] = row_values[2] = 1
# For each source in the map
for y, (source, dests) in enumerate(target_map):
lp.glp_set_row_name(problem, y+1, f'{source}->')
# First n_source auxiliary rows are fixed-bounded at 1
lp.glp_set_row_bnds(
P=problem, i=y+1, type=lp.GLP_FX, lb=1, ub=1,
)
# For each possible destination (not just those in the map)
for candidate_dest in range(n_dest):
x = y*n_dest + candidate_dest + 1
# This structural variable is a Binary Variable
lp.glp_set_col_kind(problem, x, lp.GLP_BV)
lp.glp_set_col_name(problem, x, f"{source}->{candidate_dest}")
# Objective coefficient is the sum of mismatched elements
cost = sum(
candidate_dest != actual_dest
for actual_dest in dests
)
lp.glp_set_obj_coef(problem, x, cost)
# Set the constraint matrix elements (sparse)
row_indices[1] = 1 + y # sum the sources
row_indices[2] = 1 + n_source + candidate_dest # sum the destinations
lp.glp_set_mat_col(
P=problem, j=x, len=2, ind=row_indices, val=row_values,
)
def fill_dests() -> None:
"""
For each of our candidate destination values,
- Fill the rest of the GLPK problem's "auxiliary rows" of count n_dest.
These rows correspond to constraints on the arr1 dest values.
- Name those rows and apply double bounds: each dest must get up to one
match.
"""
# For each possible destination
for candidate_dest in range(n_dest):
y = 1 + n_source + candidate_dest
lp.glp_set_row_name(problem, y, f'->{candidate_dest}')
# Last n_dest auxiliary rows are double-bounded between 0 and 1
lp.glp_set_row_bnds(
P=problem, i=y, type=lp.GLP_DB, lb=0, ub=1,
)
def create() -> None:
"""
Set some problem parameters, including objective characteristics; set the
structural and auxiliary size; call into subroutines to fill all problem
constraints.
"""
lp.glp_set_prob_name(problem, 'proximal_substitution')
lp.glp_set_obj_name(problem, 'distance')
lp.glp_set_obj_dir(problem, lp.GLP_MIN)
lp.glp_add_cols(problem, n_source * n_dest)
lp.glp_add_rows(problem, n_source + n_dest)
fill_sources()
fill_dests()
def check_code(code: int) -> None:
"""
If an optimisation function returned a failure code, attempt to convert
that into a code name and raise an OptError
"""
if code != 0:
codes = {
getattr(lp, k): k
for k in dir(lp)
if k.startswith('GLP_E')
}
raise OptError(f'glpk returned {codes[code]}')
def check_status(expected: int, actual: int) -> None:
"""
If an optimisation status is not as expected, attempt to convert the status
into a status name and raise an OptError
"""
if expected != actual:
statuses = {
getattr(lp, f'GLP_{stat}'): stat
for stat in ('OPT', 'FEAS', 'INFEAS', 'NOFEAS', 'UNBND', 'UNDEF')
}
raise OptError(
f'Optimisation failed: status {statuses.get(actual)}'
)
def solve(log_level: int = lp.GLP_MSG_ON) -> None:
"""
Perform a two-pass linear programming solver, with a simplex basis and then
a mixed-integer programming (MIP) discretisation.
"""
# First pass to calculate the basis segment: do not use the default MIP
# pre-solver because we want to use the faster primal method instead of dual
parm = lp.glp_smcp()
lp.glp_init_smcp(parm)
parm.msg_lev = log_level
check_code(lp.glp_simplex(problem, parm))
check_status(lp.GLP_FEAS, lp.glp_get_prim_stat(problem))
check_status(lp.GLP_OPT, lp.glp_get_status(problem))
# Second (very fast) MIP pass
parm = lp.glp_iocp()
lp.glp_init_iocp(parm)
parm.msg_lev = log_level
check_code(lp.glp_intopt(problem, parm))
check_status(lp.GLP_OPT, lp.glp_mip_status(problem))
def construct_substitution() -> Iterator[tuple[int, int]]:
for x in range(n_source * n_dest):
matched = lp.glp_mip_col_val(problem, x+1)
if matched:
i_source, dest = divmod(x, n_dest)
source, _ = target_map[i_source]
yield source, dest
def output(substitution: dict[int, int], log_file: Optional[str] = None) -> None:
if log_file:
lp.glp_print_mip(problem, log_file)
print('arr2, arr1, arr2new:')
print(' '.join(f'{d:3}' for d in arr2))
print(' '.join(f'{d:3}' for d in arr1))
print(' '.join(
f'{substitution[d]:3}'
for d in arr2
))
def main():
create()
t = timeit(solve, number=1)
print(f'\nSolution in {t:.2}s')
substitution = dict(construct_substitution())
output(substitution)
if __name__ == '__main__':
N = 500
random.seed(0)
arr1 = [random.randrange(N) for _ in range(N)]
arr2 = [random.randrange(N) for _ in range(N)]
target_map = make_map()
n_dest = len(arr1)
n_source = len(target_map)
problem = lp.glp_create_prob()
main()
Output
Shown with GLP_MSG_ALL:
GLPK Integer Optimizer 5.0
319 rows, 23800 columns, 47600 non-zeros
23800 integer variables, all of which are binary
Preprocessing...
319 rows, 23800 columns, 47600 non-zeros
23800 integer variables, all of which are binary
Scaling...
A: min|aij| = 1.000e+00 max|aij| = 1.000e+00 ratio = 1.000e+00
Problem data seem to be well scaled
Constructing initial basis...
Size of triangular part is 319
Solving LP relaxation...
GLPK Simplex Optimizer 5.0
319 rows, 23800 columns, 47600 non-zeros
0: obj = 2.000000000e+02 inf = 1.180e+02 (1)
Perturbing LP to avoid stalling [299]...
352: obj = 1.990000000e+02 inf = 0.000e+00 (0)
Removing LP perturbation [796]...
* 796: obj = 1.020000000e+02 inf = 0.000e+00 (0) 2
OPTIMAL LP SOLUTION FOUND
Integer optimization begins...
Long-step dual simplex will be used
+ 796: mip = not found yet >= -inf (1; 0)
+ 796: >>>>> 1.020000000e+02 >= 1.020000000e+02 0.0% (1; 0)
+ 796: mip = 1.020000000e+02 >= tree is empty 0.0% (0; 1)
INTEGER OPTIMAL SOLUTION FOUND
Writing MIP solution to 'problem-log.txt'...
arr2, arr1, arr2new:
89 111 46 15 128 119 10 152 25 179 100 51 66 91 187 120 145 43 178 172 52 196 14 173 40 41 87 135 64 30 152 113 170 44 3 120 174 104 145 130 79 166 91 99 168 64 39 143 176 3 117 189 20 85 189 11 139 71 34 61 195 123 90 156 73 172 91 151 162 158 33 183 79 99 191 106 166 20 0 152 49 178 85 40 61 57 163 114 96 181 172 145 106 8 102 179 145 107 197 169 181 11 42 114 16 66 179 40 114 135 124 143 154 193 0 9 126 83 79 119 12 106 48 140 162 21 185 33 3 102 173 106 80 0 54 3 183 193 0 172 135 156 25 48 30 155 166 50 77 71 176 46 25 121 101 160 20 5 70 115 29 65 34 167 133 166 165 88 29 39 71 4 10 10 52 174 66 142 80 93 145 10 191 179 155 167 126 182 164 117 163 111 95 137 45 53 96 150 74 2
98 194 107 10 66 130 124 103 77 122 91 149 55 129 35 72 35 193 24 158 64 136 180 154 37 79 25 186 18 175 84 120 143 25 90 111 80 156 163 52 141 122 113 133 66 15 140 3 23 184 102 181 171 160 0 156 126 85 62 186 83 180 16 48 145 56 61 36 139 114 23 20 81 130 125 27 77 141 74 180 31 140 85 138 52 154 140 150 73 113 23 152 98 81 147 61 74 47 48 47 8 156 168 66 121 17 22 173 193 33 38 9 20 179 138 174 100 180 134 70 133 60 55 173 150 107 148 70 115 126 169 164 179 91 21 83 156 29 124 150 161 85 48 62 4 187 69 29 180 56 95 43 85 109 15 25 37 178 56 11 146 162 136 154 174 18 6 31 162 48 155 147 30 100 23 94 29 9 155 5 49 47 183 31 122 53 186 15 173 5 139 108 158 25 66 17 56 18 165 77
98 108 43 10 4 70 30 84 48 61 91 149 12 129 35 111 49 193 24 150 64 136 22 169 138 79 26 33 18 175 84 120 143 32 90 111 94 105 49 52 134 69 129 130 40 18 140 3 95 90 102 181 37 160 181 156 126 56 62 186 83 54 16 85 145 150 129 36 139 114 23 20 134 130 125 60 69 37 74 84 46 24 160 138 186 154 68 66 73 113 150 49 60 81 8 61 49 72 75 47 113 156 168 66 121 12 61 138 66 33 38 3 80 179 74 82 100 180 134 70 133 60 55 87 139 107 148 23 90 8 169 60 155 74 21 90 20 179 74 150 33 85 48 55 175 122 69 29 27 56 95 43 48 109 96 25 37 178 99 11 146 162 62 53 174 69 6 31 146 140 56 147 30 30 64 94 12 9 155 5 49 30 125 61 122 53 100 15 173 102 68 108 158 112 199 17 73 116 165 77
Abridged Log
Problem: proximal_substitution
Rows: 319
Columns: 23800 (23800 integer, 23800 binary)
Non-zeros: 47600
Status: INTEGER OPTIMAL
Objective: distance = 102 (MINimum)
No. Row name Activity Lower bound Upper bound
------ ------------ ------------- ------------- -------------
1 89-> 1 1 =
2 111-> 1 1 =
3 46-> 1 1 =
4 15-> 1 1 =
5 128-> 1 1 =
...
23796 2->195 * 0 0 1
23797 2->196 * 0 0 1
23798 2->197 * 0 0 1
23799 2->198 * 0 0 1
23800 2->199 * 0 0 1
Integer feasibility conditions:
KKT.PE: max.abs.err = 0.00e+00 on row 0
max.rel.err = 0.00e+00 on row 0
High quality
KKT.PB: max.abs.err = 0.00e+00 on row 0
max.rel.err = 0.00e+00 on row 0
High quality
End of output
arr1guaranteed to always have unique elements? \$\endgroup\$