Matrices with a more general indexing than by integers A matrix is usually understood as an rectangular array of objects (for the most part numbers) arranged in rows and columns. If the objects belong to a set $\mathcal S$, then one can write $M(m,n;\mathcal S)$ for the set of $(m \times n)$-matrices $A$ with entries in $\mathcal S$. Such a matrix $A$ has $m$ rows and $n$ columns. The entry of $A$ occurring in row $i$ and column $j$ is often denoted as $a_{ij}$ and one writes $A = (a_{ij})$.
In linear algebra matrices are used to represent linear maps $f : V \to W$ between finite-dimensional vector spaces $V,W$ with respect to bases $\mathfrak V  =\{v_1,\ldots, v_n\}$ of $V$ and $\mathfrak W  =\{w_1,\ldots, w_m\}$ of $W$. We have $f(v_j) = \sum_{i=1}^m a_{ij}w_i$ with unique $a_{ij}$ and therefore $f(\sum_{j=1}^n\lambda_jv_j) =  \sum_{j=1}^n \lambda_jf(v_j) = \sum_{i,j} \lambda_j a_{ij}w_i$. The $(m\times n)$-matrix $(a_{ij})$ is the matrix representation of $f$ with respect to $\mathfrak V, \mathfrak W$.
For such matrix representations we work only with bases indexed by sets of the special form $I_k = \{1,\ldots,k\}$. Here is my question:

Wouldn't it be more flexible to allow arbitrary finite index sets? That is, to consider matrix sets of the form $M(I, J; \mathcal S) = \mathcal S^{I \times J}$ with arbitrary finite sets $I, J$. A matrix $A \in M(I, J; \mathcal S)$ is then an indexed collection $A 
  =(a_{(i,j)}) \in \mathcal S^{I \times J}$. We can still regard it as a rectangular array of objects of $\mathcal S$ arranged in rows and columns, although these do not have integer row numbers and column numbers.
Of course this concept is not a big innovation. But does is occur somewhere in the literature?

This question was motivated by the answers to  Chain rule for differentiation yields conflicting dimensions. The question deals with a matrix valued function $f : \mathbb R \to M(n,n;\mathbb R)$. Clearly $M(n,n;\mathbb R)$ is isomorphic to $\mathbb R^{n^2}$, but it does not have a canonical basis indexed by $\{1\ldots,n^2\}$. Instead it has a natural base consisting of the matrices $(E_{ij}) \in M(n,n;\mathbb R)$ where the $E_{ij}$ have an entry $1$ in row $i$ and column $j$, all other entries being $0$. The index set of this natural basis is $I_{m,n} = \{1,\ldots,m\} \times \{1,\ldots,n\}$. When considering the derivative $Df \mid_x$ at $x \in \mathbb R$, which naturally is a linear map $\mathbb R \to M(n,n;\mathbb R)$ , the mentioned answers are wavering around by saying that we have to "flatten matrices" or to "identify $M(n,n;\mathbb R)$ with $\mathbb R^{n^2}$" to get a matrix in $M(n^2,1;\mathbb R)$. One can do this, but I think it is uncessary and may even cause confusion (the OP of the above question seems to mix up $M(n^2,1;\mathbb R)$ with $M(n,n;\mathbb R)$ when flattening $Df\mid_x)$). In my opinion it is much more transparent to say that the Jacobian of $f$ at $x$ is a matrix in $M(I_{m,n},1;\mathbb R)$.
 A: Briefly, you're correct about the usefulness of arbitrary index sets; arguably, in fact, standard expositions of linear algebra do exactly as you suggest.[1] But to me the major pitfall is not indexing-by-integers versus indexing-by-finite-sets: The issue is that index sets in linear algebra are ordered even though few books note the fact explicitly.
When we speak of bases in linear algebra, we always[2] mean ordered bases. In other notation, the index set of a vector or of a basis of a finite-dimensional vector space is not the unordered set $\{1, \dots, n\}$ (say), but the ordered set $(1, \dots, n)$. When we represent a linear transformation as a matrix, ordering of bases is crucial.
Any talk of "flattening" a matrix to a column implicitly fixes a bijection from $\{1, \dots, m\} \times \{1, \dots, n\}$ to $(1, \dots, mn)$, i.e., an ordering of the double-index set.

Note 1. A bit tangentially, there is a perspective in real linear algebra where vectors are real-valued functions on an index set. This extends for example to integral operators, where
$$
Tf(x) = \int K(x, y) f(y)\, dy
$$
is a direct analogy to the matrix multiplication formula $y^{i} = \sum_{j} a_{j}^{i} x^{j}$.
Note 2. There are instances where ordering is immaterial, such as when talking about spans, but far more instances, such as the act of writing vectors and matrices as arrays, orientation, and determinants, where ordering matters but is left unsaid or half-said.
A: I decided to write an extended comment in form of a community wiki to Andrew D. Hwang's  excellent answer. I only consider the finite-dimensional case although there are analogues for infinite dimension.
The "indexing-problem" arises already on the level of representing vectors of a finite-dimensional vector space $V$ via their coordinates with respect to a basis $\mathfrak V$ of $V$. A basis of $V$ is usually defined as a subset $\mathfrak V \subset V$ which is linearly independent and spans $V$. In most definitions $\mathfrak V$ is written in the form $\{v_1,\ldots,v_n\}$. This implicitly says that we actually consider ordered bases which are finite ordered sequences $(v_1,\ldots,v_n)$ rather than finite (unordered) sets. In other words, we choose a bijecton $I_n = \{1,\ldots,n\} \to \mathfrak V$ so that the integers $i = 1,\ldots,n$ are used to index the elements of $\mathfrak V$. This seems to be quite natural because explicitly writing down the elements of $\mathfrak V$ is usually done in form of a list, i.e. in form of a linear sequential arrangement (either from left to right  or top-down).
In my opinion working with ordered bases is conceptually unneccesary for most purposes of linear algebra. An exception is the concept of orientation. In fact, each element $\mathbf x \in V$ can be written as $\mathbf x = \sum_{v \in \mathfrak V} x_v v$ with unique $x_v \in K$. In other words, when working with a basis $\mathfrak V$ we can avoid the identification $V \approx K^n$ and use the canonical isomorphism $V \approx K^{\mathfrak V}$. That is, the  coordinate representation of $\mathbf x$ with respect to $\mathfrak V$ is written abstractly as $(x_v)_{v \in \mathfrak V}$ and not as $(x_i)_{i \in I_n}$.
The problem here is that if we want to explicitly write down $(x_v)_{v \in \mathfrak V}$ for concrete numerical values of the $x_v$, then we have to agree upon a definite arrangement pattern for the $x_v$ on a sheet of paper which allows to identify the index $v$ from the position of the element $x_v \in K$ in the arrangement pattern. The usual way is of course to do it in form of a linear sequential arrangement either from left to right, giving row vectors,  or top-down, giving column vectors. This is nothing else than the choice of a bijection $I_n \to  \mathfrak V$. However, we could also agree to arrange the $x_v$ at $n$ positions of a circle (which requires to specify a function from $\mathfrak V$ to the circle). In some cases we even have a non-linear arrangement pattern by definition of $V$. This is the case for $V  =M(m,n;K)$. This vector space has a canonical basis in form of the $E_{ij}$ (see the question). The coordinate representation of a matrix $A $ with respect to this basis is naturally written in matrix pattern which gives us again $A$. Of course we can flatten it, but this would be unnatural and arbitrary.
