While your proof is basically correct, I would not consider it my favourite proof of this fact, for the following reasons:
- It uses the multiplicativity of determinants, which is a much less elementary property than invariance under transposition.
- It works only for matrices over a field, while the definition of the determinant, and the invariance under transposition, require no more than a commutative ring. There is a general principle that identities of this kind, if they hold over fields, must in fact be valid over commutative rings as well, but understanding that kind of argument requires an additional level of mathematical maturity; it is in most cases (like this one) much easier to just give a proof that is itself valid in the more general setting.
- It needs to single out the non-invertible case. If you would allow for a triangular matrix (for which the result is obvious) at the end of the product of elementary operations, you would not need to single out this case (and as a bonus you would in most cases need far fewer elementary operations). See this answer to a similar question where the use of row operations was imposed (but the answer does not say what kind of matrix one is reducing to).
How I would prove this depends on the definition of the determinant has been given. If it is defined using the Leibniz formula (which in my opinion is the right definition to give, although a different motivation should first be given), then the proof just amounts to showing that every permutation has the same sign as its inverse, which is quite elementary (and can by the way be done in a similar way as your proof, but using a decomposition into transpositions).
If on the other hand one has defined the determinant as unique $n$-linear alternating function on an $n$ dimensional vector space taking value $1$ on a given ordered basis, (which still requires using Leibniz or a substitute for proving existence) then invariance under transposition is harder to see. Here one needs to pass to the dual vector space, but it seems necessary to use the fact that if $v_1,\ldots,v_n$ are vectors and $\alpha_1,\ldots,\alpha_n$ are linear forms, then the determinant
$$
\begin{vmatrix}\alpha_1(v_1)&\ldots&\alpha_1(v_n)\\
\vdots&\ddots&\vdots\\\alpha_n(v_1)&\ldots&\alpha_n(v_n)\end{vmatrix}
$$
is not only $n$-linear and alternating in $v_1,\ldots,v_n$ for fixed $\alpha_1,\ldots,\alpha_n$ (as is clear from the definition) but also $n$-linear and alternating in $\alpha_1,\ldots,\alpha_n$ (for fixed $v_1,\ldots,v_n$); I can see no easy argument for this, other than using the Leibniz formula. Once this is established, it is easy to show that a linear operator $f$ induces the same scalar factor on $n$-linear forms when applied to $n$-uplets of vectors as it does when applied (on the right) to $n$-uplets of linear forms, in other words its determinant is the same as that of its transpose (in essence, both cases correspond to inserting an $f$ between each $\alpha_i$ and $v_j$ in the above determinant).
This still does not provide a context in which using row operations to prove this fact would be a natural choice. It is not clear to me from the question whether Artin actually suggests this, and if so why. The only reason I can think of using such an argument is when somebody commands thou shall not use the Leibniz formula.
Let $A=(a_{ij})=(\langle x_{i}, x_{j} \rangle)$ be the Gram matrix of the $x_{i}$'s. Since $A$ is symmetric and has entries in $\{0,1\}$, you can think of $A$ as the adjacency matrix of a graph on vertex set $V=\{1,2,\dots,n\}$ and edge set $E=\{\{i,j\}:x_{i}=x_{j}\}$ (i.e., vertices $i$ and $j$ are adjacent iff $x_{i}=x_{j}$). Call this graph $G_{A}.$ It's worth noting that because the diagonal of $A$ is all $1$'s, each vertex has a self-loop, but we can disregard these since the Hafnian never looks at the diagonal entries of $A.$
Now, suppose $k$ is odd and that $x_{i_{1}},x_{i_{2}},\dots,x_{i_{k}}$ are all copies of the standard basis vector $e_{i}.$ The Hafnian of $A$ counts perfect matchings of $G_{A}$, where self-loops are not allowed in a matching. The vertices $i_{1},i_{2},\dots,i_{k}$ are a clique of size $k$ in the graph, and none of the $i_{j}$'s have neighbors outside of this clique. Since $k$ is odd, there is no way to match the vertices in this clique, and thus $G_{A}$ has no perfect matchings.
You can also see this from the definition of the Hafnian, $\text{haf}(A)=\frac{1}{(n/2)!2^{(n/2)}}\sum_{\sigma \in S_{n}}\prod_{i=1}^{n-1}A_{\sigma(i)\sigma(i+1)}$, with $S_{n}$ the set of permutations of $\{1,2,\dots,n\}.$ For any $\sigma \in S_{n}$, the product $\prod_{i=1}^{n-1}A_{\sigma(i)\sigma(i+1)}$ must be $0$ since $\sigma$ must map an $i_{j}$ to some element not in $\{i_{1},i_{2},\dots,i_{k}\}.$
Best Answer
$$\begin{align} \det(AB) &= \det(\color{red}{E_1} \color{orange}{E_2} \color{blue}{E_3} \cdots E_k B) \\ &= \det(\color{red}{E_1}) \det(\color{orange}{E_2} \color{blue}{E_3} \cdots E_k B) \tag{1} \\ &= \det(\color{red}{E_1}) \det(\color{orange}{E_2}) \det(\color{blue}{E_3} \cdots E_k B) \tag{2} \\ & \ \; \vdots \\ &= \det(E_1) \det(E_2) \cdots \det(E_k) \det(B) \tag{$k$} \\ &= \det(A)\det(B) \end{align}$$ where in each of the steps $(1), (2), \dots, (k)$ we use lemma 1, since the $E_i$'s are elementary matrices; and in the last step we use lemma 2 assuming that $A = E_1\cdots E_k$.