Bad matrix multiplication that gives the right answer

Question

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.

Answer 1

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)

Answer 2

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

Answer 3

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
```

Answer 4

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.