14
\$\begingroup\$

A sample Facebook Meme

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.

\$\endgroup\$
  • 3
    \$\begingroup\$ 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! \$\endgroup\$ – QuadmasterXLII Oct 2 at 2:36
  • 1
    \$\begingroup\$ 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"? \$\endgroup\$ – mschauer Oct 3 at 16:53
  • 2
    \$\begingroup\$ 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. \$\endgroup\$ – Sisyphus Oct 4 at 6:02
  • 1
    \$\begingroup\$ 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. \$\endgroup\$ – mschauer Oct 4 at 6:17
  • 2
    \$\begingroup\$ @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) \$\endgroup\$ – Sisyphus Oct 5 at 0:25
7
\$\begingroup\$

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)
| improve this answer | |
\$\endgroup\$
3
\$\begingroup\$

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

| improve this answer | |
\$\endgroup\$
2
\$\begingroup\$

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

addmat([], [], []).
addmat([HA|TA], [HB|TB], [HC|TC]):-
    addrow(HA, HB, HC),
    addmat(TA, TB, TC).

addrow([], [], []).
addrow([HA|TA], [HB|TB], [HC|TC]):-
    HA + HB #= HC,
    addrow(TA, TB, TC).

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),
    addmat(AX, B, C),
    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
```
| improve this answer | |
\$\endgroup\$
1
\$\begingroup\$

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.

| improve this answer | |
\$\endgroup\$
  • 1
    \$\begingroup\$ 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. \$\endgroup\$ – Anders Kaseorg Oct 4 at 1:32
  • \$\begingroup\$ The next largest solutions the program produces are 9x9, and those are still record sizes so far. \$\endgroup\$ – mschauer Oct 4 at 8:42
  • 2
    \$\begingroup\$ It's a valid and high scoring solution if the title is edited to say N=9 \$\endgroup\$ – QuadmasterXLII Oct 4 at 15:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.