Find permutation of a vector that minimizes std of matrix vector product

Question

I'm trying to solve the next problem:

Given symmetric matrix A (12x12) that shows grid of a competition. A vector x (12) of rankings of the teams.

Their product gives a vector that represents total ranking of all teams with whom a team from A plays.

For example: you have 3 teams. Rankings x [1, 2, 3]. Matrix A:

      0 2 1
      2 0 4
      1 4 0

Matrix A is fixed. We need to find permutation of x so that STD(Ax) is minimal.

My previous attempt was to try to check all permutations. But it works quite long since 12!.

import itertools
import numpy as np

A = np.matrix('0,3,1,2,2,2,2,2,2,1,2,2;3,0,3,1,2,3,1,1,2,2,1,2;1,3,0,2,2,2,2,2,2,1,2,2;2,1,2,0,2,1,2,2,2,2,3,2;2,2,2,2,0,2,1,2,2,3,2,1;2,3,2,1,2,0,3,2,1,2,1,2;2,1,2,2,1,3,0,2,1,1,3,3;2,1,2,2,2,2,2,0,3,2,2,1;2,2,2,2,2,1,1,3,0,3,1,2;1,2,1,2,3,2,1,2,3,0,2,2;2,1,2,3,2,1,3,2,1,2,0,2;2,2,2,2,1,2,3,1,2,2,2,0')
min = 1000000
for x in itertools.permutations([2433,2057,1935,1927,1870,1841,1818,1770,1680,1497,1435,1289]):
    x = np.matrix(x).T
    b = A.dot(x)
    cur = np.std(b)
    if cur < min:
        min = cur
        res = x

I know there is scipy minimize but I dont know weather it can cope with permutations of x instead of continuous optimization.

The question is how to solve this task as fast and precise as possible.

Thanks.

score 1 · Accepted Answer · answered May 15 '21 at 16:51

You can write a much faster brute-force implementation.

Firstly, you can use a matrix multiplication rather than many dot product by working on chunks of permutations. Matrix multiplication kernel are heavily optimized and thus run much faster than many dot product.

Moreover, you can partially pre-compute the permutations to speed up the computation even more by splitting permutations in two parts. The idea is to first build an index which contains all the permutation consisting in picking 5 element among the 12. Then, the idea is to find all the permutation of an array of 7 item (the indices and not the values themselves). Finally, all the permutations can be build from the two index.

Note that further optimizations are possible when the two above optimizations are applied together: one can compute the matrix multiplication more efficiently if one part of the permutation is constant.

The resulting algorithm is complex but much more efficient then the original one. Here is the code:

def computeOptim(A):
    mini = 1000000

    permValues = np.array([2433, 2057, 1935, 1927, 1870, 1841, 1818, 1770, 1680, 1497, 1435, 1289])

    # Precompute partial permutations: high and low part of all the permutations.
    loPerms = np.array(list(itertools.permutations(range(7))))
    hiPerms = np.array(list(itertools.permutations(range(12), 5)))

    # Iterate over chunks (of 7!=5040 permutations)
    for hiPerm in hiPerms:
        # Find the remaining index to include in the low-part permutations
        loPermIndices = np.array(list(set(range(12))-set(hiPerm)))

        # Find all the possible low-part permutations for the current 
        # high-part permutation by reindexing the values.
        curLoPerms = loPermIndices[loPerms]

        # Compute the chunks of possible x values
        loPermValues = permValues[curLoPerms]
        hiPermValues = permValues[hiPerm]

        # A matrix multiplcation is used to compute many dot product in a row.
        # Compute effciently  B = A @ X  with  X the matrix containing all the permutations
        hiB = A[:,:len(hiPermValues)] @ hiPermValues[None,:].T
        loB = A[:,len(hiPermValues):] @ loPermValues.T
        B = hiB + loB

        multiCur = np.std(B, axis=0)
        minPos = np.argmin(multiCur)

        if multiCur[0,minPos] < mini:
            mini = multiCur[0,minPos]
            res = np.concatenate((hiPermValues, loPermValues[minPos]))

A = np.matrix('0,3,1,2,2,2,2,2,2,1,2,2;3,0,3,1,2,3,1,1,2,2,1,2;1,3,0,2,2,2,2,2,2,1,2,2;2,1,2,0,2,1,2,2,2,2,3,2;2,2,2,2,0,2,1,2,2,3,2,1;2,3,2,1,2,0,3,2,1,2,1,2;2,1,2,2,1,3,0,2,1,1,3,3;2,1,2,2,2,2,2,0,3,2,2,1;2,2,2,2,2,1,1,3,0,3,1,2;1,2,1,2,3,2,1,2,3,0,2,2;2,1,2,3,2,1,3,2,1,2,0,2;2,2,2,2,1,2,3,1,2,2,2,0')
computeOptim(A)

It succeed to find the optimal solution in 50 seconds on my machine while the original code would take roughly 5h30. As a result this code is about 400 faster.

The optimal solution found is:

mini = 291.80729942892106
res = [2433 1841 1289 1818 2057 1927 1770 1870 1497 1680 1935 1435]

Thank you! Could you provide a link where the trick with the permutations is explained? — Sergei Shumilin, May 16 '21 at 07:09
Well I am not sure there is a link to explain specifically that. I created this for the purpose of solving this question. However, the mathematical idea behind this is simple: picking 12 balls from a bag is strictly equivalent to pick 5 balls and then the 7 others. From the programming point of view, I guess it is possible to extrapolate this algorithm from the standard [recursive one](https://algotree.org/algorithms/recursive/permutations_recursion/) (due to [symbolic expression computing](https://stackoverflow.com/questions/16395704) that are indices here). — Jérôme Richard, May 16 '21 at 12:40
For the 7/5, I tested different values. I choose 6/6 first and found that 7/5 is faster and takes less memory. The idea is that the matrix multiplication should work on big enough matrices so it can be much faster than a dot product. However, the matrices should not be too big because then they do not fit in the CPU cache slowing down the computation. Moreover, the size of `loPerms` and `hiPerms` also matters. Pre-computing all the permutation would take too much memory (and slow too). Ideally, they should be equally big to minimize the memory footprint. 7/5 and 8/4 are good for that. — Jérôme Richard, May 16 '21 at 12:50

Find permutation of a vector that minimizes std of matrix vector product

1 Answers1