Go Back > Fun Stuff > Puzzles

Thread Tools
Old 2019-12-04, 16:17   #89
Jan 2017

23×11 Posts

Originally Posted by bsquared View Post
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)
                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:
            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:
            add18(r[i], r[c1], -r[i][c1])
    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:
        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:
        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:
        p1 = pairs[i] >> 8
        p2 = pairs[i] & 255
        if m[p1] == m[p2]:
        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:

import pyximport

from det_optim import optim

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

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):

def unit(B, bestd):
    for _ in range(100):
        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).
uau is offline   Reply With Quote
Old 2019-12-04, 19:18   #90
Romulan Interpreter
LaurV's Avatar
Jun 2011

2·3·5·313 Posts

Very Nice! Congrats to all.

I am a bit ashamed to say my solution this time hihi
LaurV is offline   Reply With Quote

Thread Tools

Similar Threads
Thread Thread Starter Forum Replies Last Post
November 2018 Batalov Puzzles 5 2018-12-03 13:31
November 2017 Batalov Puzzles 3 2017-12-08 14:55
November 2016 Xyzzy Puzzles 1 2016-12-06 16:41
November 2015 R. Gerbicz Puzzles 3 2015-12-01 17:48
November 2014 Xyzzy Puzzles 1 2014-12-02 17:40

All times are UTC. The time now is 10:38.

Fri Apr 23 10:38:09 UTC 2021 up 15 days, 5:19, 0 users, load averages: 1.30, 1.72, 1.68

Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2021, Jelsoft Enterprises Ltd.

This forum has received and complied with 0 (zero) government requests for information.

Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.2 or any later version published by the Free Software Foundation.
A copy of the license is included in the FAQ.