[proofplan]
We count invertible matrices by counting their columns. An $n \times n$ matrix over $\mathbb{F}_q$ is invertible exactly when its ordered list of columns is an ordered basis of the [vector space](/page/Vector%20Space) $\mathbb{F}_q^n$. The first column may be any nonzero vector, and after $k$ linearly independent columns have been chosen, the next column may be any vector outside their $k$-dimensional span. Multiplying the number of choices at each stage gives the stated product.
[/proofplan]
[step:Identify invertible matrices with ordered bases of $\mathbb{F}_q^n$]
Let $V:=\mathbb{F}_q^n$, viewed as an $n$-dimensional vector space over $\mathbb{F}_q$, and let $M_n(\mathbb{F}_q)$ denote the set of all $n \times n$ matrices with entries in $\mathbb{F}_q$. For each matrix $A \in M_n(\mathbb{F}_q)$ and each $j \in \{1,\dots,n\}$, let $c_j(A) \in V$ denote the $j$th column of $A$. Define the column map
\begin{align*}
\Gamma: M_n(\mathbb{F}_q) \to V^n, \qquad A \mapsto (c_1(A),\dots,c_n(A)).
\end{align*}
This map is bijective, since an $n \times n$ matrix is uniquely determined by its ordered list of columns.
We claim that $A \in GL_n(\mathbb{F}_q)$ if and only if $(c_1(A),\dots,c_n(A))$ is an ordered basis of $V$. Let
\begin{align*}
L_A: V \to V, \qquad x \mapsto Ax
\end{align*}
be the [linear map](/page/Linear%20Map) represented by $A$ in the standard basis $(e_1,\dots,e_n)$ of $V$. For every $j \in \{1,\dots,n\}$, we have $L_A(e_j)=c_j(A)$. If $A$ is invertible, then $L_A$ is a vector space isomorphism, so it sends the basis $(e_1,\dots,e_n)$ to the basis $(c_1(A),\dots,c_n(A))$. Conversely, if $(c_1(A),\dots,c_n(A))$ is a basis of $V$, then $L_A$ sends a basis to a basis, hence is an isomorphism, so $A$ is invertible. Therefore $|GL_n(\mathbb{F}_q)|$ equals the number of ordered bases of $V$.
[/step]
[step:Count the possible choices for each successive column]
We count ordered bases $(v_1,\dots,v_n)$ of $V$. Since $|\mathbb{F}_q|=q$, the vector space $V=\mathbb{F}_q^n$ has $q^n$ elements.
For the first vector $v_1$, the only forbidden choice is $0$, so there are $q^n-1$ choices.
Now suppose $k \in \{1,\dots,n-1\}$ and linearly independent vectors $v_1,\dots,v_k \in V$ have already been chosen. Define their span by
\begin{align*}
W_k:=\operatorname{span}_{\mathbb{F}_q}\{v_1,\dots,v_k\} \le V.
\end{align*}
Because $v_1,\dots,v_k$ are linearly independent, the coordinate map
\begin{align*}
\theta_k: \mathbb{F}_q^k \to W_k, \qquad (a_1,\dots,a_k) \mapsto \sum_{i=1}^{k} a_i v_i
\end{align*}
is bijective. Hence $|W_k|=|\mathbb{F}_q^k|=q^k$. The next vector $v_{k+1}$ must lie outside $W_k$, and every vector outside $W_k$ gives a linearly independent list $(v_1,\dots,v_k,v_{k+1})$. Thus there are $q^n-q^k$ choices for $v_{k+1}$.
[guided]
We are building an ordered basis one vector at a time. The condition at stage $k+1$ is not that the next vector be nonzero by itself, but that it avoid all linear combinations of the previously chosen vectors.
Start with $v_1$. Since $V=\mathbb{F}_q^n$ has $q^n$ elements and the zero vector cannot appear in a linearly independent list, $v_1$ has $q^n-1$ possible choices.
Now assume that $k \in \{1,\dots,n-1\}$ and that $v_1,\dots,v_k$ have already been chosen linearly independently. Define the subspace
\begin{align*}
W_k:=\operatorname{span}_{\mathbb{F}_q}\{v_1,\dots,v_k\}.
\end{align*}
Because the chosen vectors are linearly independent, every element of $W_k$ has a unique expression as a linear combination of them. Equivalently, the map
\begin{align*}
\theta_k: \mathbb{F}_q^k \to W_k, \qquad (a_1,\dots,a_k) \mapsto \sum_{i=1}^{k} a_i v_i
\end{align*}
is bijective. Since $\mathbb{F}_q^k$ has $q^k$ elements, this gives $|W_k|=q^k$.
The vector $v_{k+1}$ must not belong to $W_k$. If $v_{k+1}\in W_k$, then it is a linear combination of $v_1,\dots,v_k$, so the list is linearly dependent. If $v_{k+1}\notin W_k$, then no nontrivial relation among $v_1,\dots,v_k,v_{k+1}$ can exist: solving such a relation for $v_{k+1}$ would place $v_{k+1}$ in $W_k$, unless the coefficient of $v_{k+1}$ were zero, and then the independence of $v_1,\dots,v_k$ forces all remaining coefficients to be zero. Therefore precisely the vectors outside $W_k$ are allowed, and their number is
\begin{align*}
|V\setminus W_k|=|V|-|W_k|=q^n-q^k.
\end{align*}
[/guided]
[/step]
[step:Multiply the independent choices to obtain the order]
By the multiplication principle for successive finite choices, the number of ordered bases of $V$ is
\begin{align*}
(q^n-1)(q^n-q)(q^n-q^2)\cdots(q^n-q^{n-1}).
\end{align*}
Equivalently,
\begin{align*}
\prod_{k=0}^{n-1}(q^n-q^k).
\end{align*}
By the identification of $GL_n(\mathbb{F}_q)$ with the set of ordered bases of $V$ through columns, this is exactly $|GL_n(\mathbb{F}_q)|$. Hence
\begin{align*}
|GL_n(\mathbb{F}_q)|=\prod_{k=0}^{n-1}(q^n-q^k).
\end{align*}
This proves the theorem.
[/step]