[proofplan]
We prove that row-major vectorization is a linear bijection. Linearity follows directly from the entrywise definitions of matrix addition, scalar multiplication, and coordinate-wise operations in $k^{mn}$. For bijectivity, we construct an explicit inverse map by filling an $m\times n$ matrix row by row from a vector in $k^{mn}$, and then verify both composites are identity maps.
[/proofplan]
[step:Verify that row-major vectorization is linear]
Let $A=(A_{ij})$ and $B=(B_{ij})$ be elements of $M_{m\times n}(k)$, and let $\lambda\in k$. Matrix addition and scalar multiplication in $M_{m\times n}(k)$ are entrywise, so for every $1\le i\le m$ and $1\le j\le n$,
\begin{align*}
(A+B)_{ij}=A_{ij}+B_{ij}
\end{align*}
and
\begin{align*}
(\lambda A)_{ij}=\lambda A_{ij}.
\end{align*}
Using the defining formula for $\operatorname{vec}$, the coordinate of $\operatorname{vec}(A+B)$ with index $(i-1)n+j$ is
\begin{align*}
(\operatorname{vec}(A+B))_{(i-1)n+j}=(A+B)_{ij}=A_{ij}+B_{ij}=(\operatorname{vec} A)_{(i-1)n+j}+(\operatorname{vec} B)_{(i-1)n+j}.
\end{align*}
Since this equality holds for every coordinate of $k^{mn}$, we have
\begin{align*}
\operatorname{vec}(A+B)=\operatorname{vec} A+\operatorname{vec} B.
\end{align*}
Similarly, for every $1\le i\le m$ and $1\le j\le n$,
\begin{align*}
(\operatorname{vec}(\lambda A))_{(i-1)n+j}=(\lambda A)_{ij}=\lambda A_{ij}=\lambda(\operatorname{vec} A)_{(i-1)n+j}.
\end{align*}
Thus
\begin{align*}
\operatorname{vec}(\lambda A)=\lambda\operatorname{vec} A.
\end{align*}
Therefore $\operatorname{vec}$ is $k$-linear.
[/step]
[step:Construct the row-filling inverse map]
Define a map $\Phi:k^{mn}\to M_{m\times n}(k)$ as follows. For a vector $x=(x_1,\dots,x_{mn})\in k^{mn}$, let $\Phi(x)$ be the matrix whose $(i,j)$-entry is
\begin{align*}
(\Phi(x))_{ij}=x_{(i-1)n+j}
\end{align*}
for all $1\le i\le m$ and $1\le j\le n$. This defines a unique element of $M_{m\times n}(k)$, because the formula assigns exactly one scalar in $k$ to each matrix position.
[/step]
[step:Show that the row-filling map is a two-sided inverse]
Let $A=(A_{ij})\in M_{m\times n}(k)$. For every $1\le i\le m$ and $1\le j\le n$,
\begin{align*}
(\Phi(\operatorname{vec} A))_{ij}=(\operatorname{vec} A)_{(i-1)n+j}=A_{ij}.
\end{align*}
Hence $\Phi(\operatorname{vec} A)=A$, so $\Phi\circ\operatorname{vec}=\operatorname{id}_{M_{m\times n}(k)}$.
Now let $x=(x_1,\dots,x_{mn})\in k^{mn}$. For each coordinate index $q\in\{1,\dots,mn\}$, there are unique indices $i\in\{1,\dots,m\}$ and $j\in\{1,\dots,n\}$ such that
\begin{align*}
q=(i-1)n+j.
\end{align*}
For these indices,
\begin{align*}
(\operatorname{vec}(\Phi(x)))_q=(\operatorname{vec}(\Phi(x)))_{(i-1)n+j}=(\Phi(x))_{ij}=x_{(i-1)n+j}=x_q.
\end{align*}
Thus $\operatorname{vec}(\Phi(x))=x$, and therefore $\operatorname{vec}\circ\Phi=\operatorname{id}_{k^{mn}}$.
[guided]
The inverse should undo exactly the row-major ordering used by $\operatorname{vec}$. Since $\operatorname{vec}$ sends the entry in row $i$ and column $j$ to coordinate $(i-1)n+j$, the map $\Phi:k^{mn}\to M_{m\times n}(k)$ is defined by reversing this assignment: for $x=(x_1,\dots,x_{mn})\in k^{mn}$, the matrix $\Phi(x)$ has entries
\begin{align*}
(\Phi(x))_{ij}=x_{(i-1)n+j}.
\end{align*}
We first check the composite $\Phi\circ\operatorname{vec}$. Let $A=(A_{ij})\in M_{m\times n}(k)$. The $(i,j)$-entry of $\Phi(\operatorname{vec} A)$ is obtained by taking coordinate $(i-1)n+j$ of $\operatorname{vec} A$. By the definition of row-major vectorization,
\begin{align*}
(\Phi(\operatorname{vec} A))_{ij}=(\operatorname{vec} A)_{(i-1)n+j}=A_{ij}.
\end{align*}
This holds for every matrix position $(i,j)$, so the matrices have the same entries and hence are equal:
\begin{align*}
\Phi(\operatorname{vec} A)=A.
\end{align*}
Therefore $\Phi\circ\operatorname{vec}=\operatorname{id}_{M_{m\times n}(k)}$.
We next check the other composite, $\operatorname{vec}\circ\Phi$. Let $x=(x_1,\dots,x_{mn})\in k^{mn}$. Every coordinate index $q\in\{1,\dots,mn\}$ corresponds to a unique row-column pair: the unique $i\in\{1,\dots,m\}$ and $j\in\{1,\dots,n\}$ satisfying
\begin{align*}
q=(i-1)n+j.
\end{align*}
This is exactly the Euclidean division of $q-1$ by $n$, with quotient $i-1$ and remainder $j-1$. For this pair, the $q$-th coordinate of $\operatorname{vec}(\Phi(x))$ is
\begin{align*}
(\operatorname{vec}(\Phi(x)))_q=(\operatorname{vec}(\Phi(x)))_{(i-1)n+j}=(\Phi(x))_{ij}=x_{(i-1)n+j}=x_q.
\end{align*}
Since every coordinate agrees, $\operatorname{vec}(\Phi(x))=x$. Thus $\operatorname{vec}\circ\Phi=\operatorname{id}_{k^{mn}}$.
[/guided]
[/step]
[step:Conclude that vectorization is an isomorphism]
The map $\operatorname{vec}:M_{m\times n}(k)\to k^{mn}$ is $k$-linear, and the map $\Phi:k^{mn}\to M_{m\times n}(k)$ satisfies
\begin{align*}
\Phi\circ\operatorname{vec}=\operatorname{id}_{M_{m\times n}(k)}
\end{align*}
and
\begin{align*}
\operatorname{vec}\circ\Phi=\operatorname{id}_{k^{mn}}.
\end{align*}
Therefore $\operatorname{vec}$ is a bijective $k$-[linear map](/page/Linear%20Map), hence an isomorphism of $k$-vector spaces.
[/step]
