 ## Executive summary:

Print NxN matrices $$\A\$$ and $$\B\$$ where $$\ AB = 10A + B \$$. Largest N wins!

## Details:

Entrants in this challenge should choose an N. Submit a program that takes an input integer 'seed', and produces two full rank NxN matrices $$\A\$$ and $$\B\$$, where each element of $$\A\$$ and $$\B\$$ is an integer from 0 to 9, and $$\AB = 10A + B\$$. In addition, matrices A and B must have at least half of their elements nonzero. The output should always be the same for a given seed. Runtime of your program is not limited, however as evidence that it does finish you must post a sample seed, the output of your program given that seed, and the time taken to compute.

Input and output format is not restricted.

## Scoring:

First, largest N wins. Second, ties are broken by earliest submission.

## No hard coding:

In order to prevent hard coding, there must be at least 10 seeds that produce distinct outputs, where outputs are not distinct if they are equal after exchanging rows, exchanging columns, and or transposing the matrix. You don't need to worry about proving this for good faith solutions.

## Math Background / No Rules Below This Line

First of all, for a given candidate A, B is fixed: $$AB = 10A + B\\ A = 10AB^{-1} + I\\ I = 10B^{-1} + A^{-1}\\ B = 10(I - A^{-1})^{-1}\\ B = 10(I + (A - I)^{-1})$$

One approach is to take the eigendecomposition of A.

From before, $$A = Q \Lambda_A Q^{-1}\\ B = 10(I - A^{-1})^{-1}\\ B = 10(Q I Q^{-1} - Q \Lambda_A^{-1} Q^{-1})^{-1}\\ B = Q(\frac{10}{(I - \Lambda_A^{-1})})Q^{-1}\\$$

This shows that A, B share eigenvectors Q, and their eigenvalues are assosciated by $$\Lambda_A \Lambda_B = 10 \Lambda_A + \Lambda_B\\$$

Interestingly, this proves that AB = BA, ie our matrices commute!

We can calculate the determinant of B as $$\10^N \frac{1}{\Pi_i (1 - \frac{1}{\lambda_{A,i}})}\$$ This is very suggestive, but doesn't immediately yield a way to pick a matrix A such that B is small positive integers.

Because A and B share eigenvectors, B can be written as a linear combination of $$\(I, A, A^2 ... ~ A^{N-1})\$$. For example,

$$A = \left(\begin{array}{rrr}1 & 2 & 4 \\1 & 1 & 3 \\1 & 1 & 1\end{array}\right), B = \left(\begin{array}{rrr}7 & 4 & 6 \\3 & 6 & 4 \\1 & 2 & 8\end{array}\right)$$

$$\Lambda_A = \left(\begin{array}{rrr}-\sqrt{6} + 2 & 0 & 0 \\0 & \sqrt{6} + 2 & 0 \\0 & 0 & -1\end{array}\right), \Lambda_B = \left(\begin{array}{rrr}-\sqrt{6} + 2 & 0 & 0 \\0 & \sqrt{6} + 2 & 0 \\0 & 0 & -1\end{array}\right)$$

We observe that $$\\Lambda_B = \Lambda_A^2 - 2 \Lambda_A + 2I\$$, so $$\B = A^2 - 2A + 2\$$. In this case B is forced to be an integer matrix because it is possible to write B as an integer polynomial of A. If we can find large matrices where $$\10(I + (\Lambda_A - I)^{-1})\$$ is an integer polynomial of $$\\Lambda_A\$$, then that might provide an efficient solution.

Finally, if we write $$\A\$$ as $$\E + I\$$, then we can see some simplifications. $$\B = 10 * (E^{-1} + I)\$$, so B is an integer matrix as long as the determinant of E divides 10. With this E representation, our problem can become to find an integer matrix E such that $$E + I \geq 0\\ E^{-1} + I \geq 0\\ |E| = \pm 10\\$$

(Occasionally this will fail if an entry of $$\ E^{-1} + I > 9\$$ )

This represents an improvement, because once we find a [0-9] matrix E with determinant $$\\pm 10 \$$ (such as by random search), then if there are only a small number of negative entries in $$\E^{-1}\$$ we can run determinant-preserving transformations on E such as permuting rows to move those negative entries onto the diagonal, where they are made positive by the $$\+I\$$

In addition, there is a name and wikipedia page for matrices like $$\E^{-1}\$$ where the off diagonal entries are all positive: A Metzler Matrix. Further literature search in that direction might bring up an efficient way to generate positive integer matrices with a determinant $$\\pm 10\$$ whose inverses are Metzler Matrices.

• I changed the scoring criteria to be largest N followed by earliest submission, and added a specific loophole closure to answer Sisyphus's point. Thanks! – QuadmasterXLII Oct 2 '20 at 2:36
• Can the program produce matrices of variable size, as long as I provide a seed to a record with large N and it is "good faith"? – mschauer Oct 3 '20 at 16:53
• I was hopeful that an SMT solver could make light work of this problem, but disappointing results so far. I can only solve up to N ~ 7. – Sisyphus Oct 4 '20 at 6:02
• I don't see the fundamental difference of searching through the very finite space of all 2x2 matrices with random search or searching through the space of Toeplitz matrices with constant off diagonal through random search, but it's up to you... I like my solution because it is out of the box. I did not do any computational work which is not reflected by my program, I manually found a good seed and hardcoding the seed was allowed. – mschauer Oct 4 '20 at 6:17
• @user Here is my code, there's no special ideas in it. On my PC it finds 6x6 matricies in ~2 minutes and 7x7 matricies in a couple of hours. Note it will not always generate a valid matrix, you need to discard the invalid ones (sometimes too many zeroes) – Sisyphus Oct 5 '20 at 0:25

# Python, N = 7

Using variable neighborhood search over elementary row and column operations. Not every seed works quickly (or at all?), but seed 9 finds a solution for N = 7 in a few seconds.

$time python3 search.py 7 9 [[0 3 6 6 1 6 6] [0 1 5 5 0 5 5] [1 3 1 2 1 8 1] [0 0 1 0 0 1 0] [1 7 9 9 0 9 9] [1 3 1 1 1 0 1] [0 0 1 1 0 1 1]] [[5 5 1 1 0 4 0] [1 7 0 0 1 0 0] [0 0 9 4 0 1 5] [0 0 0 5 0 0 5] [2 2 1 1 7 4 0] [0 0 1 1 0 9 0] [0 2 0 0 0 0 0]] real 0m4.909s user 0m4.883s sys 0m0.026s  ### Code import itertools import sys import numpy as np def search(depth, C, D, best): A = C + I B = D + 10 * I x = np.r_[A, B].flatten() score = ( np.maximum(0, np.maximum(x - 9, 1 - x)).sum() - min((A == 0).sum(), -(-n * n // 2)) - min((B == 0).sum(), -(-n * n // 2)) + 2 * n - np.linalg.matrix_rank(A) - np.linalg.matrix_rank(B) ) if best is None or score < best: return C, D, score if depth == 0: return None tried = {} for i in reversed(range(n * (n - 1) * 5)): j = rng.integers(i + 1) mutation, tried[j] = tried.get(j, j), tried.get(i, i) p, mutation = divmod(mutation, (n - 1) * 5) q, mutation = divmod(mutation, 5) q = (p + q + 1) % n C1 = C.copy() D1 = D.copy() if mutation == 0: C1[p, :] += C1[q, :] D1[:, q] -= D1[:, p] elif mutation == 1: C1[p, :] -= C1[q, :] D1[:, q] += D1[:, p] elif mutation == 2: C1[:, p] += C1[:, q] D1[q, :] -= D1[p, :] elif mutation == 3: C1[:, p] -= C1[:, q] D1[q, :] += D1[p, :] elif p < q: C1[[p, q], :] = C1[[q, p], :] D1[:, [p, q]] = D1[:, [q, p]] else: C1[:, [p, q]] = C1[:, [q, p]] D1[[p, q], :] = D1[[q, p], :] ret = search(depth - 1, C1, D1, best) if ret is not None: return ret return None n, seed = map(int, sys.argv[1:]) rng = np.random.default_rng(seed) I = np.eye(n, dtype=int) score = None k = rng.integers(0, 8, n) C = np.diag(np.array([-10, -5, -2, -1, 1, 2, 5, 10])[k]) D = np.diag(np.array([-1, -2, -5, -10, 10, 5, 2, 1])[k]) depth = 1 while score != 0: ret = search(depth, C, D, score) if ret is not None: C, D, score = ret depth = 0 depth += 1 A = C + I B = D + 10 * I assert (A @ B == 10 * A + B).all() print(A) print(B)  N = 2 This is a naive brute force approach in R: randomvalid <- function(N) { while (T) { X <- matrix(sample(10,N*N,T)-1,N) if (sum(X==0) < N*N/2 && qr(X)$rank == N)
return(X)
}
}

tryforever <- function(N) {
i <- 1
while (T) {
i <- i + 1
if (i %% 100000 == 0) print(sprintf("after %d tries...", i))
A <- randomvalid(N)
B <- randomvalid(N)
if(all(A %*% B == 10*A + B)) {
print(A)
print(B)
print(A %*% B)
}
}
}

# run tryforever(2) to find some solutions for N=2


Some solutions for N=2 are:

A  = [[ 3, 2],[ 6, 2]]
B  = [[ 9, 2],[ 6, 8]]
AB = [[39,22],[66,28]]

A  = [[ 6, 5],[ 3, 2]]
B  = [[ 9, 5],[ 3, 5]]
AB = [[69,55],[33,25]]

A  = [[ 6, 5],[ 9, 8]]
B  = [[ 3, 5],[ 9, 5]]
AB = [[63,55],[99,85]]

A  = [[ 3, 4],[ 8, 7]]
B  = [[ 7, 2],[ 4, 9]]
AB = [[37,42],[84,79]]


Valid solutions seem to be very rare, even for N=2 they pop up only once every few minutes. For N=3 I was not able to find any solution within hours by this approach.

So, let's see whether a mathematical or heuristical approach is able to beat my N=2 solution :)

# N = 4, Swi-Prolog

After all the math, my best result so far is a relatively naive solution in prolog.

:- use_module(library(clpfd)).

% N is the dot product of lists V1 and V2.
dot(V1, V2, N) :- maplist(product,V1,V2,P), my_sumlist(P,N).
product(N1,N2,N3) :- N3 #= N1*N2.

my_sumlist([], 0).
my_sumlist([H|T], N) :- my_sumlist(T, X), N #= X + H.

my_scamul(_, [], []).

my_scamul(K, [R|T], [R2|T2]):-
my_rowmul(K, R, R2),
my_scamul(K, T, T2).

my_rowmul(_, [], []).
my_rowmul(K, [V|T], [V2|T2]):-
K * V #= V2,
my_rowmul(K, T, T2).

HA + HB #= HC,

nn(NN) :-
NN #= 4.

isMat([], 0).
isMat([H|T], N):-
nn(NN),
isRow(H, NN),
isMat(T, Ns),
Ns #= N - 1.

isRow([], 0).
isRow([H|T], N):-
H in 1..9,
isRow(T, Ns),
Ns #= N - 1.

det([[E_0, E_1, E_2, E_3],
[E_4, E_5, E_6, E_7],
[E_8, E_9, E_10, E_11],
[E_12, E_13, E_14, E_15]], 4, D):-
D #= E_1*E_11*E_14*E_4 - E_1*E_10*E_15*E_4 - E_11*E_13*E_2*E_4 + E_10*E_13*E_3*E_4 - E_0*E_11*E_14*E_5 + E_0*E_10*E_15*E_5 + E_11*E_12*E_2*E_5 - E_10*E_12*E_3*E_5 - E_1*E_11*E_12*E_6 + E_0*E_11*E_13*E_6 + E_1*E_10*E_12*E_7 - E_0*E_10*E_13*E_7 - E_15*E_2*E_5*E_8 + E_14*E_3*E_5*E_8 + E_1*E_15*E_6*E_8 - E_13*E_3*E_6*E_8 - E_1*E_14*E_7*E_8 + E_13*E_2*E_7*E_8 + E_15*E_2*E_4*E_9 - E_14*E_3*E_4*E_9 - E_0*E_15*E_6*E_9 + E_12*E_3*E_6*E_9 + E_0*E_14*E_7*E_9 - E_12*E_2*E_7*E_9.

% Matrix multiplication with matrices represented
% as lists of lists. M3 is the product of M1 and M2
mmult(M1, M2, M3) :- transpose(M2,MT), maplist(mm_helper(MT), M1, M3).
mm_helper(M2, I1, M3) :- maplist(dot(I1), M2, M3).

problem(A, B, C, ABf):-
nn(NN),
isMat(A, NN),
isMat(B, NN),
my_scamul(10, A, AX),
mmult(A, B, C),
append(A, B, AB),
flatten(AB, ABf).

main(Seed):-
problem(A, B, C, ABf),
labeling([ff, random_value(Seed)], ABf),
det(A, 4, Determinant),
write_ln(A),
write_ln(B),
write_ln(C),
write_ln(Determinant),
Determinant #\= 0.


Sample Output with seed 20:

?- time(main(20))
|    .
[[8,8,3,6],[4,2,1,2],[5,3,1,2],[4,2,1,2]]
[[9,1,1,2],[1,5,1,2],[1,1,3,6],[1,5,1,2]]
[[89,81,31,62],[41,25,11,22],[51,31,13,26],[41,25,11,22]]
0
[[8,8,3,6],[4,2,1,2],[5,3,1,2],[4,3,1,1]]
[[9,2,1,1],[1,6,1,1],[1,4,3,3],[1,1,1,6]]
[[89,82,31,61],[41,26,11,21],[51,34,13,23],[41,31,11,16]]
-6
% 582,520,793 inferences, 33.623 CPU in 33.623 seconds (100% CPU, 17324911 Lips)
true
$$$$


# Julia, 9

I identify a class of full rank matrices for which the criterion A*B = 10A + B can be checked in constant time, nxn-matrices of the form

   a b b b . . . b
b a b b       .
b b a b       .
b b b a       .
.       .     .
.         .   .
.           a b
b . . . . . b a


Then I do a random search for those matrices in this class fulfilling the criteria. As these are full, I can make block-matrices out of two random elements of the class get a half-filled larger solution.

using Random, LinearAlgebra
check(a1, a2, b1, b2, n, base) = a2 != 0 && a1*b1*n + a2*b1 + b2*a1 == base*a1 + b1 && a1*b1*n + a2*b1 + b2*a1 + a2*b2 == base*(a1 +a2) + b1 + b2

function generate(base=10)
while true
a1 = rand(1:base-1)
b1 = rand(1:base-1)
a2 = rand(-a1+1:base-1-a1)
b2 = rand(-b1+1:base-1-b1)
n = rand(1:10)
if check(a1, a2, b1, b2, n, base)
return a1*ones(Int, n, n) + a2*I, b1*ones(Int, n, n) + b2*I
end
end
end
nosp(A) = sum(A .== 0) <= size(A,1)*size(A,2)/2
fake(A, B) = rank(A) > 0 && rank(B) > 0 && nosp(A) && nosp(B) && A*B == 10*A + B

Random.seed!(333333333333); A1, B1 = generate()
A2, B2 = generate()
A = [A1 0I; 0I A2]
B =[B1 0I; 0I B2]
fake(A2, B2)

 1  2  2  2  2  2  0  0  0  0  0  0
2  1  2  2  2  2  0  0  0  0  0  0
2  2  1  2  2  2  0  0  0  0  0  0
2  2  2  1  2  2  0  0  0  0  0  0
2  2  2  2  1  2  0  0  0  0  0  0
2  2  2  2  2  1  0  0  0  0  0  0
0  0  0  0  0  0  1  2  2  2  2  2
0  0  0  0  0  0  2  1  2  2  2  2
0  0  0  0  0  0  2  2  1  2  2  2
0  0  0  0  0  0  2  2  2  1  2  2
0  0  0  0  0  0  2  2  2  2  1  2
0  0  0  0  0  0  2  2  2  2  2  1

6  1  1  1  1  1  0  0  0  0  0  0
1  6  1  1  1  1  0  0  0  0  0  0
1  1  6  1  1  1  0  0  0  0  0  0
1  1  1  6  1  1  0  0  0  0  0  0
1  1  1  1  6  1  0  0  0  0  0  0
1  1  1  1  1  6  0  0  0  0  0  0
0  0  0  0  0  0  6  1  1  1  1  1
0  0  0  0  0  0  1  6  1  1  1  1
0  0  0  0  0  0  1  1  6  1  1  1
0  0  0  0  0  0  1  1  1  6  1  1
0  0  0  0  0  0  1  1  1  1  6  1
0  0  0  0  0  0  1  1  1  1  1  6
`

There are not too many valid pairs of them but there are ten distributed over the dimensions for base 10. It also produces only a single 12x12 solution therefore it was agreed to count the method at the level of the second largest solutions it produces, ie 9.

• It’s worth noting that $\left(\left[\begin{smallmatrix}1&2&2&2&2&2\\2&1&2&2&2&2\\2&2&1&2&2&2\\2&2&2&1&2&2\\2&2&2&2&1&2\\2&2&2&2&2&1\end{smallmatrix}\right], \left[\begin{smallmatrix}6&1&1&1&1&1\\1&6&1&1&1&1\\1&1&6&1&1&1\\1&1&1&6&1&1\\1&1&1&1&6&1\\1&1&1&1&1&6\end{smallmatrix}\right]\right)$ is the only such pair of 6×6 matrices. All others are 3×3 or smaller. – Anders Kaseorg Oct 4 '20 at 1:32
• The next largest solutions the program produces are 9x9, and those are still record sizes so far. – mschauer Oct 4 '20 at 8:42
• It's a valid and high scoring solution if the title is edited to say N=9 – QuadmasterXLII Oct 4 '20 at 15:08