Bad matrix multiplication that gives the right answer 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 &amp; 2 &amp; 4 \\1 &amp; 1 &amp; 3 \\1 &amp; 1 &amp; 1\end{array}\right), B = \left(\begin{array}{rrr}7 &amp; 4 &amp; 6 \\3 &amp; 6 &amp; 4 \\1 &amp; 2 &amp; 8\end{array}\right)$$
$$\Lambda_A = \left(\begin{array}{rrr}-\sqrt{6} + 2 &amp; 0 &amp; 0 \\0 &amp; \sqrt{6} + 2 &amp; 0 \\0 &amp; 0 &amp; -1\end{array}\right), \Lambda_B = \left(\begin{array}{rrr}-\sqrt{6} + 2 &amp; 0 &amp; 0 \\0 &amp; \sqrt{6} + 2 &amp; 0 \\0 &amp; 0 &amp; -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 &gt; 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.
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 &lt; 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 &lt; 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)