mersenneforum.org November 2019
 Register FAQ Search Today's Posts Mark Forums Read

2019-12-04, 16:17   #89
uau

Jan 2017

4F16 Posts

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

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:
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
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
return
if abs(r[c1][c1]) < 1e-6:
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
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

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

 2019-12-04, 19:18 #90 LaurV Romulan Interpreter     Jun 2011 Thailand 112·73 Posts Very Nice! Congrats to all. I am a bit ashamed to say my solution this time hihi

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

All times are UTC. The time now is 13:15.

Wed Oct 21 13:15:50 UTC 2020 up 41 days, 10:26, 1 user, load averages: 1.03, 1.34, 1.46