mersenneforum.org  

Go Back   mersenneforum.org > Fun Stuff > Puzzles

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

7410 Posts
Default

Quote:
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.

Quote:
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:
Code:
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:
Code:
#!/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).
uau is offline   Reply With Quote
Old 2019-12-04, 19:18   #90
LaurV
Romulan Interpreter
 
LaurV's Avatar
 
Jun 2011
Thailand

19·461 Posts
Default

Very Nice! Congrats to all.

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

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 16:48.

Wed Sep 30 16:48:01 UTC 2020 up 20 days, 13:58, 0 users, load averages: 2.03, 1.86, 1.80

Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2020, 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.