A Hodge dual computation on a $4$-dimensional Riemannian manifold Let $(M,g)$ be a $4$-dimensional smooth Riemannian manifold.  I am trying to understand the following exterior algebra computation:
Let $x^1,x^2,x^3,x^4$ be local coordinates on $M$ such that the Riemannian volume form of $g$ is $\mathrm{d}x^1\wedge\mathrm{d}x^2\wedge\mathrm{d}x^3\wedge\mathrm{d}x^4$.  
Then there exist a local frame $\alpha_1,\alpha_2,\alpha_3,\alpha_4$ of $1$-forms such that the following identities hold:
$$
\begin{aligned}
\ast_g(\mathrm{d}x^1\wedge\mathrm{d}x^2)&= \alpha_3\wedge\alpha_4\\
\ast_g(\mathrm{d}x^1\wedge\mathrm{d}x^3)&= \alpha_4\wedge\alpha_2\\
\ast_g(\mathrm{d}x^1\wedge\mathrm{d}x^4)&= \alpha_2\wedge\alpha_3\\
\end{aligned}\qquad
\begin{aligned}
\ast_g(\mathrm{d}x^3\wedge\mathrm{d}x^4)&= \alpha_1\wedge\alpha_2\\
\ast_g(\mathrm{d}x^4\wedge\mathrm{d}x^2)&= \alpha_1\wedge\alpha_3\\
\ast_g(\mathrm{d}x^2\wedge\mathrm{d}x^3)&= \alpha_1\wedge\alpha_4\\
\end{aligned}
\tag1
$$
(The pattern is that the indices on both sides of every equation are complementary.)
In fact, the solution of these equations (which is unique up to replacing each $\alpha_i$ by $-\alpha_i$) is given by 
$$
\alpha_i = g_{ij}\,\mathrm{d}x^j\qquad\qquad\text{where}\quad 
g = g_{ij}\,\mathrm{d}x^i\mathrm{d}x^j.
$$

Question: Why does the formula $\alpha_i = g_{ij}\,\mathrm{d}x^j$ solve the equations?

Here is what I understood: 
$\mathrm{d}x^1,\mathrm{d}x^2$ are $g$-orthogonal to $\alpha_3,\alpha_4$. Since $\ast_g(\mathrm{d}x^1\wedge\mathrm{d}x^2)$ is decomposable*, we can write 
$$
\ast_g(\mathrm{d}x^1\wedge\mathrm{d}x^2)=\beta_1 \wedge \beta_2$$ 
for some one-forms $\beta_i$, which implies $\mathrm{d}x^1,\mathrm{d}x^2$ are $g$-orthogonal to $\beta_1,\beta_2$. Thus $\text{span}\{\beta_1,\beta_2\}=\text{span}\{\alpha_3,\alpha_4\}$, so $\beta_1 \wedge \beta_2=f\alpha_3 \wedge \alpha_4$ for some function $f$. We now need to prove that $f=1$, which is (up to a sign) equivalent to the statement
$$
\| \ast_g(\mathrm{d}x^1\wedge\mathrm{d}x^2)\|=\|\alpha_3 \wedge \alpha_4\|.
$$
Since the Hodge dual operator is an isometry, this is equivalent to
$$
\| \mathrm{d}x^1\wedge\mathrm{d}x^2\|=\|\alpha_3 \wedge \alpha_4\|.
$$
This is where the assumption $\det(g_{ij})=1$ is supposed to enter. I tried to expand both sides in terms of the $g_{ij}$, but so far I don't see how the result follows.
Perhaps there is another easier way to see this.
For the interested, this computation came up in a question about the existence of "higher-order" harmonic coordinates.
Comment: Of course, everything here is "pointwise", i.e. this is really a result about $4$-dim inner product spaces. I have kept the manifold notation since it might be more familiar. 

*The Hodge star operator preserve decomposability of elements.
 A: You have not defined the Hodge dual precisely, so I am gonna assume that you use the "determinant" convention, meaning that $$ (\omega\wedge\eta)_{i_1...i_kj_1...j_l}=\frac{(k+l)!}{k!l!}\omega_{[i_1...i_k}\eta_{j_1...j_l]}. $$ In this convention $$ (\star\omega)_{i_{k+1}...i_n}=\frac{1}{k!}\epsilon_{i_1...i_n}\omega^{i_1...i_k} $$ where indices are raised and lowered with $g^{ij}$ and $g_{ij}$ respectively and $$ \epsilon_{i_1...i_n}=\sqrt{\det g}\pi_{i_1...i_n} $$ with $\pi$ being the Levi Civita symbol.
If $\omega\in\Omega^k(M)$, the relationship between the tensor and differential form representations is $$ \omega=\omega_{i_1...i_k}dx^{i_1}\otimes...\otimes dx^{i_k}=\frac{1}{k!}\omega_{i_1...i_k}dx^{i_1}\wedge...\wedge dx^{i_k} $$ where the summation is unrestricted in both cases and the array of components is completely antisymmetric.
How does one evaluate $\star dx^{i_1}\wedge...\wedge dx^{i_k}$ then? Well, $$ dx^{i_1}\wedge...\wedge dx^{i_k}=\delta^{i_1}_{j_1}...\delta^{i_k}_{j_k}dx^{j_1}\wedge...\wedge dx^{j_k}=\frac{k!}{k!}\delta^{i_1}_{[j_1}...\delta^{i_k}_{j_k]}dx^{j_1}\wedge...\wedge dx^{j_k}, $$ so the tensor components are $$ k!\delta^{i_1}_{[j_1}...\delta^{i_k}_{j_k]}. $$
The Hodge dual is then $$ (\star dx^{i_1}\wedge...\wedge dx^{i_k})_{j_{k+1}...j_n}=\delta^{i_1}_{[j_1}...\delta^{i_k}_{j_k]}\epsilon^{j_1...j_n}=g^{i_1j_1}...g^{i_k j_k}\epsilon_{j_1...j_n}=\epsilon^{i_1...i_k}\ _{j_{k+1}...j_n}, $$ hence $$ \star dx^{i_1}\wedge...\wedge dx^{i_k}=\frac{1}{k!} \epsilon^{i_1...i_k}\ _{j_{k+1}...j_n} dx^{j_{k+1}}\wedge...\wedge dx^{j_n}=\sum_{<}\epsilon^{i_1...i_k}\ _{j_{k+1}...j_n} dx^{j_{k+1}}\wedge...\wedge dx^{j_n},$$ where the sum with $<$ underneath is restricted summation (to lexicographically ordered indices).
Let us apply this in your case:
$$ \star dx^i\wedge dx^j=\sum_{k<l}\epsilon^{ij}\ _{kl}dx^k\wedge dx^l=\sum_{k<l}\epsilon^{ijkl}\alpha_k\wedge\alpha_l $$ where $\alpha_k=g_{ik}dx^i$.
Now if we fix the indices and use $\sqrt{\det g}=1$, we obtain that $$ \star dx^1\wedge dx^2=\sum_{k<l}\pi^{12kl}\alpha_k\wedge\alpha_l=\pi^{1234}\alpha_3\wedge\alpha_4=\alpha_3\wedge\alpha_4, $$ and so on. In fact, one may obtain a bit more general result without having to fix the indices.
Since if any two indices agree, we'll get a vanishing expression, let us assume that $i$ and $j$ are different. We can then assume that $k,l$ are different from $i,j$, thus allowing two different values for $k,l$. If $k,l$ can take on two different values, then for $\sum_{k<l}\epsilon^{ijkl}\alpha_k\wedge\alpha_l$ there is only one term, we don't actually sum. If we then further assume that we order indices such that $i,j,k,l$ is an even permutation of $1,2,3,4$, we can also get rid of $\pi$, since it will be identially 1. So we obtain $$ \star dx^i\wedge dx^j=\alpha_k\wedge\alpha_l, $$ with $i,j,k,l$ all different and are positive permutations of $1,2,3,4$.
