November 2019

[uau]

Originally Posted by bsquared

Oleg was kind enough to post his solution and code, which still seem like magic to me.

I was actually considering a similar approach, but never got around to implementing it. The basics of it are actually simple if you consider calculating the determinant by making an orthogonal version of the matrix. (As in: go through rows, and for each row subtract a multiple of the row from all following rows such that all scalar products between rows are 0. Determinant is then the product of row vector lengths.)

By such considerations, it's pretty obvious that ideally you'd want all row vectors to have about the same length (using bigger numbers for already long vectors has diminishing payoff). And dot products between rows to be minimal. I verified that the known best matrix for size 8 from OEIS does fulfill those conditions - near same vector lengths (but not exact - square sum from 202 to 208), and dot products between rows from 155 to 158. It didn't seem like Oleg's solution used much math beyond that.

find_best_swap() benefited greatly from uau's hint because it nearly eliminates the innermost determinant at the cost of a single inverse at the beginning; probably he was doing something similar to the above.

Basically, except I didn't search for the best swap, just did the first possible swap. I did test best swap, but searching for it first mainly seemed to slow things down. When always doing first swap found, avoiding calculating the inverse from scratch too was extra helpful as for a new matrix the first searches find a swap very quickly.


Below is the optimization code implemented in Cython:

cdef void add18(double a[18], double b[18], double c):
    cdef int i
    for i in range(18):
        a[i] += c * b[i]

cdef int inverse(int mat[81], double r[9][18]):
    cdef int i, j, k
    cdef double det = 1
    k = 0
    for i in range(9):
        for j in range(9):
            r[i][j] = mat[k]
            k += 1
    for i in range(9):
        for j in range(9, 18):
            r[i][j] = 0
    for i in range(9):
        r[i][i+9] = 1

    cdef double a
    for i in range(9):
        if abs(r[i][i]) < 1e-6:
            for j in range(i+1, 9):
                if abs(r[j][i]) > 2e-6:
                    add18(r[i], r[j], 1)
                    break
            else:
                return 0
        a = 1/r[i][i]
        det *= r[i][i]
        for j in range(18):
            r[i][j] *= a
        r[i][i] = 1 # hopefully better numerical stability
        for j in range(9):
            if i == j:
                continue
            add18(r[j], r[i], -r[j][i])
    return int(abs(det)+0.5)

cdef void adjust_inverse(double r[9][18], int r1, int c1, int r2, int c2, int v):
    cdef int i
    for i in range(9):
        r[i][c1] += v * r[i][r1+9]
    for i in range(9):
        r[i][c2] -= v * r[i][r2+9]
    cdef double a
    if c1 == c2:
        a = 1/r[c1][c1]
        for i in range(c1, 18):
            r[c1][i] *= a
        for i in range(9):
            if i == c1:
                continue
            add18(r[i], r[c1], -r[i][c1])
        return
    if abs(r[c1][c1]) < 1e-6:
        add18(r[c1], r[c2], 1)
    a = 1/r[c1][c1]
    for i in range(18):
        r[c1][i] *= a
    r[c1][c1] = 1 # more stable?
    for i in range(9):
        if i == c1:
            continue
        add18(r[i], r[c1], -r[i][c1])
    a = 1/r[c2][c2]
    for i in range(18):
         r[c2][i] *= a
    r[c2][c2] = 1
    for i in range(9):
        if i == c2:
            continue
        add18(r[i], r[c2], -r[i][c2])

def optim(arr, pairs_in):
    cdef int i, j, k, v
    cdef int m[81]
    cdef int div[81]
    cdef short pairs[3240]
    for i in range(9):
        for j in range(9):
            div[i*9+j] = i
    for i in range(81):
        m[i] = arr[i]
    for i in range(3240):
        pairs[i] = pairs_in[i]
    cdef double inv[9][18]
    cdef int det, d2
    cdef double r
    det = inverse(m, inv)
    if det == 0:
        return 0
    cdef int time, prev, p1, p2
    cdef int r1, c1, r2, c2

    prev = 3239 # If this is only improvement on first pass, not detected
    i = -1

    while True:
        i += 1
        if i == 3240:
            i = 0
        if i == prev:
            break
        p1 = pairs[i] >> 8
        p2 = pairs[i] & 255
        if m[p1] == m[p2]:
            continue
        r1 = div[p1]
        c1 = p1 - 9*r1
        r2 = div[p2]
        c2 = p2 - 9*r2
        v = m[p2] - m[p1]
        r = (1+inv[c1][r1+9]*v)*(1+inv[c2][r2+9]*-v)+inv[c2][r1+9]*inv[c1][r2+9]*v*v
        d2 = int(det * abs(r) + .5)
        if d2 > det:
            m[p1] += v
            m[p2] -= v
            prev = i
            det = d2
            adjust_inverse(inv, r1, c1, r2, c2, v);
    for i in range(81):
        arr[i] = m[i]
    return det

And pure Python calling it:

#!/usr/bin/python3

import pyximport
pyximport.install(language_level=3)

from det_optim import optim

import random
B = []
for i in range(1, 10):
    B.extend([i]*9)

pairs = [0]*3240
k = 0
for i in range(81):
    for j in range(i+1, 81):
        pairs[k] = (i << 8) + j
        k += 1

pl = []
for _ in range(100):
    random.shuffle(pairs)
    pl.append(tuple(pairs))

def unit(B, bestd):
    for _ in range(100):
        random.shuffle(B)
        x = B[:]
        for pairs in pl:
            B = x[:]
            d = optim(B, pairs)
            if d >= bestd:
                bestd = d
                print(d, B)
                assert d < 935000000
    return bestd

bestd = 0
#for _ in range(1):
while True:
    bestd = unit(B, bestd)

If you have Cython installed, you can save these in the same directory and run the second. The first needs to be named "det_optim.pyx" so that import statement in the second finds the right file.

The optimization routine takes the order to try swaps in as an argument. The reason for that was that the optimization of a matrix is fast enough that using Python's random.shuffle() to generate a new random matrix was a significant portion of time used, and trying swaps in a different order seemed to work about as well as always generating a completely new random matrix.

If you comment out the 'while True:' in the second file and uncomment the line above, you can benchmark it to see how long it takes to optimize 10000 matrices (note that the Cython code will be compiled to C, so the first run will use time for that, but it should be cached for subsequent runs).

[LaurV]

Very Nice! Congrats to all.

I am a bit ashamed to say my solution this time hihi