Given the full matrix ring over a finite field, $M := M_{n \times n}(\mathbb{F}_q)$ for prime $q$ and integer $n$, what can one say about subsets $S$ of $M$ satisfying the condition that:
$A,B \in S$ implies $A - B$ is invertible unless $A = B$ (n.b. I'm not making any assumptions on the invertibility of elements of $S$ themselves).
Is it possible to construct, or show the existence of, relatively large subsets satisfying this condition? I'm particularly interested in the case where $q$ scales asymptotically and $n$ is a fixed, small constant, and by large subsets I'm attempting to construct $S$ with $\vert S \vert \approx O(q^n)$, but suspect this is infeasible.
As for my own attempts, I thought initially that given the collection of all subspaces of $\mathbb{F}_q^n$ of dimension $n-1$, one ought to be able to choose one element from each subspace so that the difference of any two elements is full rank, and thus invertible. This would lead to a set $S$ of size ${n \choose n-1}_q = 1+ \dots + q^{n-1}$, but I am unable to prove that this is either achievable or impossible.
Asymptotically, I've only been able to come up with obvious naive sets of size $q$, but playing around with small cases suggests that larger sets are possible in some cases, just they do not exhibit an obvious pattern (at least not obvious to me).
The largest possible such a set $S$ has exactly $q^n$ elements:
In small cases it is easy to describe an embedding (there are several). For example, if $n=2$ and $q$ is an odd prime we can find a quadratic non-residue $\epsilon\in\Bbb{F}_q$. Then the collection of matrices of the form $$ M(a,b)=\left(\begin{array}{cc} a& \epsilon b\\ b&a\end{array}\right) $$ forms a subfield $L=\Bbb{F}_{q^2}$. The determinant $\det M(a,b)=a^2-\epsilon b^2$ vanishes only when $a=b=0$. Furthermore, the difference of two such matrices is of the same form. This is a generalization of the well known way of representing complex numbers by $2\times 2$ real matrices: $$a+ib\mapsto\left(\begin{array}{cc}a&-b\\ b&a\end{array}\right),$$ where we use $-1$ as a non-square.
When $n=3$, $p\neq3$ we can similarly use a non-cube $\epsilon$ and matrices of the form $$\left(\begin{array}{ccc} a&\epsilon c&\epsilon b\\ b&a&\epsilon c\\ c&b&a\end{array}\right).$$
For larger $n$ the method of describing an explicit set of matrices becomes a bit more complicated. If $m(x)$ is the minimal polynomial of a generator of the extension field, we can use $\Bbb{F}_q$-linear combinations of the powers $A^i$, $i=0,1,2\ldots,n-1$, of the companion matrix $A$ of $m(x)$.