Why are matrices indexed by row first? I am wondering is there some greater thought behind indexing matrices in form (row, column), or is it just something that someone did, and everybody else followed?
When I'm looking at the matrix, I find it easiest to imagine it as numbers laid out in 2D grid. From my previous education I usually follow (x,y) order, and it seems like in many other areas we usually deal with x before y. When I'm thinking ${A_{3 1}}$, I'd think x = 3, y = 1, and therefore get 3rd column and 1st row. This would be wrong - I need to do otherwise(y's first).
Why so?
 A: While I believe that the true reason is purely historical, here’s a possible explanation for why the common notation might be more natural, and it boils down to the following facts:


*

*We like to think of the composition $g ∘ f$ of maps $f \colon X → Y$, $g \colon Y → Z$ as “$g$ after $f$” not the other way around (viz. to ensure that “$(g∘f)(x) = g(f(x))$” holds, which is intuitive).

*We like to think of matrices as linear maps in a “covariant way”.

*We like to think of vectors as columns (to avoid wasting precious horizontal space (because we write from left to right, not from top to bottom)).

*We like to cancel stuff out in the middle rather than outside, because that sticks better in my (our?) memory.


So let $R$ be a ring (or a field, take $R = ℚ$ in need) and read on.
In linear Algebra, we often use matrices $A ∈ \operatorname{Mat}_{m×n}(R)$ as maps that are given by multiplication
$$f_A \colon R^n → R^m,~x ↦ Ax,$$
for which we actually need to think of vectors in $R^n$ and $R^m$ as columns.
Now, we could also think of them as maps
$$f_A^{op} \colon R^m → R^n,~x ↦ xA,$$
but then we would


*

*lose the nice “covariant” identity $f_A∘f_B = f_{AB}$ (and would end up instead with the “contravariant” identity $f_B^{op}∘f_A^{op} = f_{AB}^{op}$), and

*need to think of $R^m$ and $R^n$ as rows rather than columns or

*redefine the matrix product to multiply columns of the first matrix onto rows of the second (instead of rows of the first onto columns of the second).


While the last two points might not be an issue, the first one is for many people. So we favor the interpretation $f_A$ over the interpretation $f_A^{op}$.
Okay, so now if you write out a product of a matrix by a vector, you really want to preserve horizontal space, so you rather think of the vector as being a column. We also would like to consider this product to be a special case of a general matrix product, considering vectors as special matrices.
This forces us to multiply matrices by the scheme “dot product of rows of the first matrix by columns of the second matrix”. In particular, the number of rows of a product is determined by its first factor, the number of columns by its second factor.
Now, if we were to annotate by $(j,i)$ the entry at the $j$-th column and $i$-th row of a matrix, the definition of the matrix product
$$(c_{ji})_{n×m} = (a_{ki})_{q×m}·(b_{jk})_{n×q}$$
would be
$$c_{ji} = \sum_{k=1}^q a_{ki}·b_{jk},$$
so you wouldn’t be able to memorize it as “the middle index is cancelling out” which, to me, sounds way more natural than “the outer indices are cancelling out”.
Of course, you sacrifice the naturality of having the two-dimensional grid, reminiscent of the cartesian plane, as you suggested. But since we are already writing matrices top-left to bottom-right (because that’s the way we write) as others already pointed out in the comments, you can’t really think of a matrix in the same way as the cartesian plane anyway, so why try to force it? The way it is, it at least preserves the mathematical positive orientation and if you tilt your head by $-π/2$, you again have your cartesian plane!
A: I don't know the historical reason, but this is how I understand it intuitively based on how vectors and matrices are used in modern mathematics and computer graphics.
Vectors are most commonly notated in column major (vertical) form:
$$
\begin{bmatrix}
x \\
y \\
z \\
\end{bmatrix}
$$
When looking up matrix transformations for example, the overwhelming majority you find will be given using column vectors. For example, from wikipedia:
$$
\begin{bmatrix}
  x' \\
  y' \\
\end{bmatrix} = \begin{bmatrix}
  \cos\theta & -\sin\theta \\
  \sin\theta &  \cos\theta \\
\end{bmatrix}\begin{bmatrix}
  x \\
  y \\
\end{bmatrix}
$$
Since vectors are thought of as columns, it makes sense that you can just consider a matrix to be multiple column vectors together:
$$ \require{color}
\begin{bmatrix}
\colorbox{lightgreen}{$x_1$} & \colorbox{lightblue}{$x_2$} & \colorbox{yellow}{$x_3$} \\
\colorbox{lightgreen}{$y_1$} & \colorbox{lightblue}{$y_2$} & \colorbox{yellow}{$y_3$} \\
\colorbox{lightgreen}{$z_1$} & \colorbox{lightblue}{$z_2$} & \colorbox{yellow}{$z_3$} \\
\end{bmatrix}
$$
In this way, the notation of M-by-N or Row-by-Column makes perfect sense. The dimensionality of the vectors (i.e 2D, 3D) is most important, the number of vectors is in the matrix is of secondary importance.
To illustrate:
$$ \require{color}
\begin{bmatrix}
\colorbox{lightgreen}{$x_1$} & \colorbox{lightblue}{$x_2$} & \colorbox{yellow}{$x_3$} \\
\colorbox{lightgreen}{$y_1$} & \colorbox{lightblue}{$y_2$} & \colorbox{yellow}{$y_3$} \\
\end{bmatrix}
$$
These are 2D vectors and we have 3 of them, therefore this is a $2 \times 3$ matrix.
$$ \require{color}
\begin{bmatrix}
\colorbox{lightgreen}{$x_1$} & \colorbox{lightblue}{$x_2$} & \colorbox{yellow}{$x_3$} & \colorbox{orange}{$x_4$} & \colorbox{pink}{$x_5$} \\
\colorbox{lightgreen}{$y_1$} & \colorbox{lightblue}{$y_2$} & \colorbox{yellow}{$y_3$} & \colorbox{orange}{$y_4$} & \colorbox{pink}{$y_5$} \\
\colorbox{lightgreen}{$z_1$} & \colorbox{lightblue}{$z_2$} & \colorbox{yellow}{$z_3$} & \colorbox{orange}{$z_4$} & \colorbox{pink}{$z_5$} \\
\end{bmatrix}
$$
These are 3D vectors and we have 5 of them, therefore this is a $3 \times 5$ matrix.
