calculating the wedge product for 2 4D vectors
Solution 1:
The wedge product is not a matrix, although a matrix representation of the coordinates with respect to a specific basis is possible. The wedge product does have an area representation, but it is signed and oriented area that can be thought of as lying in the plane of the span of the wedge vectors. Your numerical example can also be computed readily. I'll attempt to answer some of these questions below in stages. I will start off more generally, assuming that you might be interested in arbitrary general vector spaces (with a possibly non-Euclidean metric), and then for the numerical example, return to a Euclidean space.
Coordinate expansion of a wedge product
Given a basis for a non-Euclidean N-dimensional vector space $ \left\{ {f_\mu} \right\} $, and it's reciprocal basis $ \left\{ {f^\mu} \right\} $ where $ f_\mu \cdot f^\nu = {\delta_\mu}^\nu $, the decomposition of a vector in this space along each of the basis directions has the form (sums implied) $$x = \left( { x \cdot f^\mu } \right) f_\mu.$$ Observe that $ x \cdot f^\mu $ are the coordinates of the vector $ x $ with respect to the basis $ \left\{ {f_\mu} \right\} $. A visualization of such a decomposition with respect to a pair of vectors is illustrated in fig. 1.1.
Writing $$y = \left( { y \cdot f^\mu } \right) f_\mu,$$ let's expand the wedge product $ x \wedge y $ in terms their coordinates for this specific basis. We have $$\begin{aligned}x \wedge y&=\left( { \left( { x \cdot f^\mu } \right) f_\mu } \right) \wedge\left( { \left( { y \cdot f^\nu } \right) f_\nu } \right) \\ &=\left( { x \cdot f^\mu } \right)\left( { y \cdot f^\nu } \right)\left( {f_\mu \wedge f_\nu} \right) \\ &=\sum_{\mu \ne \nu}\left( { x \cdot f^\mu } \right)\left( { y \cdot f^\nu } \right)\left( {f_\mu \wedge f_\nu} \right) \\ &=\sum_{\mu < \nu}\left( { x \cdot f^\mu } \right)\left( { y \cdot f^\nu } \right)\left( {f_\mu \wedge f_\nu} \right)+\sum_{\nu > \mu}\left( { x \cdot f^\nu } \right)\left( { y \cdot f^\mu } \right)\left( {f_\nu \wedge f_\mu} \right) \\ &=\sum_{\mu < \nu}\left( {\left( { x \cdot f^\mu } \right)\left( { y \cdot f^\nu } \right)-\left( { y \cdot f^\mu } \right)\left( { x \cdot f^\nu } \right)} \right)\left( {f_\mu \wedge f_\nu} \right),\end{aligned}$$ where we have used the fact that $ f_\mu \wedge f_\mu = 0 $, for any $ \mu $, have performed a dummy index swap in the second sum, and finally, used the antisymmetry property $ f_\nu \wedge f_\mu = - f_\mu \wedge f_\nu $ to group terms. We see that the coordinates of each of the basis bivectors has a determinant structure $$ x \wedge y =\sum_{\mu < \nu}\begin{vmatrix} x \cdot f^\mu & x \cdot f^\nu \\ y \cdot f^\mu & y \cdot f^\nu\end{vmatrix}f_\mu \wedge f_\nu.$$
Area interpretation of the wedge product
We may reduce this sum to its simplest form by introducing an orthonormal "standard" basis $ \left\{ {e_1, e_2} \right\} $ for the subspace containing the span of $ x, y $ such that $ x $ is colinear with $ e_1 $. That is $$\begin{aligned} x &= \left( { x \cdot e^1 } \right) e_1 \\ y &= \left( { y \cdot e^1 } \right) e_1 + \left( { y \cdot e^2 } \right) e_2.\end{aligned}$$ Given that we are allowing the underlying vector space to be non-Euclidean, by "orthonormal" here, we mean $ e_\mu \cdot e_\mu = \pm 1$, for $ \mu \in \left\{ {1,2} \right\} $, and $ e_1 \cdot e_2 = 0 $. In terms of this basis, the wedge of the two vectors expands to $$ x \wedge y =\begin{vmatrix} x \cdot e^1 & 0 \\ y \cdot e^1 & y \cdot e^2\end{vmatrix}e_1 \wedge e_2.$$ We see that such a basis choice "diagonalizes" the determinant coefficient of the one remaining wedge product basis bivector, leaving just $$ x \wedge y =\left( { x \cdot e^1 } \right)\left( { y \cdot e^2 } \right)e_1 \wedge e_2.$$ If we form a parallelogram by placing $ y $ at the head of $ x $ and then complete a loop back along $ -x $ and $ -y $, as illustrated in fig. 1.2, we see that $ x \cdot e^1 $ is the length of the base of the parallelogram, and $ y \cdot e^2 $ is the length of the projection of $ y $ along the perpendicular to $ x $.
In the example above the wedge product's numerical coefficient is positive, but if we had computed $ y \wedge x $, as illustrated in fig. 1.3, the base vector is now $ -x \cdot e^1 e_1 $, so $$ y \wedge x =\left( { -x \cdot e^1 } \right)\left( { y \cdot e^2 } \right)e_1 \wedge e_2= - (x \wedge y),$$ as required algebraically, so we say that $ x \wedge y $ is a signed area with respect to the orientation of the basic $ \left\{ { e_1, e_2 } \right\} $. In particular, if the wedge is positive, it has the same orientation as the square formed by applying this parallelogram construction recipe to the vectors $ e_1, e_2 $ as illustrated in fig. 1.4.
Observe that the wedge product's signed area interpretation is also basis dependent. For example, had we selected $ f_1 = -e_1, f_2 = e_2 $ as our basis vectors, then the numerical coefficient $ \left( { x \cdot f^1 } \right) \left( { y \cdot f^2 } \right) $ would flip sign accordingly, and we would say that $ x \wedge y $ has a negative sign with respect to the bivector $ f_1 \wedge f_2 $.
The numerical example from the original question
With $\mathbf{a} = \mathbf{i} + 3 \mathbf{j} + 4 \mathbf{k} + 5 \mathbf{l}$ and $\mathbf{b} = 6 \mathbf{i} + 7 \mathbf{j} + 8 \mathbf{k} + 9 \mathbf{l}$, we have $$\begin{aligned}B &= \mathbf{a} \wedge \mathbf{b} \\ &= \left( {2 \mathbf{i} + 3 \mathbf{j} + 4 \mathbf{k} + 5 \mathbf{l}} \right) \wedge \left( { 6 \mathbf{i} + 7 \mathbf{j} + 8 \mathbf{k} + 9 \mathbf{l} } \right) \\ &=14 \mathbf{i} \wedge \mathbf{j} + 16 \mathbf{i} \wedge \mathbf{k} + 18 \mathbf{i} \wedge \mathbf{l} \\ &\quad+18 \mathbf{j} \wedge \mathbf{i} + 24 \mathbf{j} \wedge \mathbf{k} + 27 \mathbf{j} \wedge \mathbf{l} \\ &\quad+24 \mathbf{k} \wedge \mathbf{i} + 28 \mathbf{k} \wedge \mathbf{j} + 36 \mathbf{k} \wedge \mathbf{l} \\ &\quad+30 \mathbf{l} \wedge \mathbf{i} + 35 \mathbf{l} \wedge \mathbf{j} + 40 \mathbf{l} \wedge \mathbf{k} \\ &=(14 -18) \mathbf{i} \wedge \mathbf{j}+ (16-24) \mathbf{i} \wedge \mathbf{k}+ (18-30) \mathbf{i} \wedge \mathbf{l}+ (24-28) \mathbf{j} \wedge \mathbf{k}+ (27-35) \mathbf{j} \wedge \mathbf{l}+ (36-40) \mathbf{k} \wedge \mathbf{l} \\ &=4 \mathbf{j} \wedge \mathbf{i} + 8 \mathbf{k} \wedge \mathbf{i} + 12 \mathbf{l} \wedge \mathbf{i} + 4 \mathbf{k} \wedge \mathbf{j} + 8 \mathbf{l} \wedge \mathbf{j} + 4 \mathbf{l} \wedge \mathbf{k} \\ &= \left( { 8 \sqrt{5} } \right)\frac{\mathbf{j} \wedge \mathbf{i} + 2 \mathbf{k} \wedge \mathbf{i} + 3 \mathbf{l} \wedge \mathbf{i} + \mathbf{k} \wedge \mathbf{j} + 2 \mathbf{l} \wedge \mathbf{j} + \mathbf{l} \wedge \mathbf{k}}{2 \sqrt{5} }.\end{aligned}$$
In the last step, the wedge product has been factored into a scalar part $\sqrt{320} = 8 \sqrt{5}$ and a unit bivector part (assuming a Euclidean metric.) The scalar factor is the area, and the unit bivector (i.e. squares to $\pm 1$) provides the orientation.
To confirm that the bivector on the right is a unit bivector, one can compute the (geometric algebra) square of the numerator: $$\left( { \mathbf{j} \wedge \mathbf{i} + 2 \mathbf{k} \wedge \mathbf{i} + 3 \mathbf{l} \wedge \mathbf{i} + \mathbf{k} \wedge \mathbf{j} + 2 \mathbf{l} \wedge \mathbf{j} + \mathbf{l} \wedge \mathbf{k} } \right)^2= -1 -4 -9 -1 -4 -1 = -20 = -\left( { 2 \sqrt{5} } \right)^2.$$
Now, does an area of $ \sqrt{320} $ make sense? We can compute that area using standard vector algebra methods to confirm that this is the case. First use Gram-Schmidt to compute the portion of $ \mathbf{b} $ that is perpendicular to $ \mathbf{a} $, which is: $$\mathbf{b}_\perp = \mathbf{b} - \left( { \mathbf{a} \cdot \mathbf{b} } \right) \frac{\mathbf{a}}{\mathbf{a} \cdot \mathbf{a}} = (52/27) \mathbf{i} + (8/9) \mathbf{j} -(4/27) \mathbf{k} -(32/27) \mathbf{l}.$$ Using this, we can compute the (squared) area (with $\mathbf{x}^2 \equiv \mathbf{x} \cdot \mathbf{x}$), to find $$\mathbf{a}^2 \mathbf{b}_\perp^2 = 320,$$ as expected.
You also saw that the wedge product magnitude should have the form $$\left\lVert {\mathbf{a}} \right\rVert\left\lVert {\mathbf{b}} \right\rVert \sin\theta.$$ That is easy to show, first calculating $$\cos\theta = \frac{ \mathbf{a} \cdot \mathbf{b} }{ \left\lVert {\mathbf{a}} \right\rVert \left\lVert {\mathbf{b}} \right\rVert },$$ after which we find $$\left\lVert {\mathbf{a}} \right\rVert \left\lVert {\mathbf{b}} \right\rVert \sin\theta = 8 \sqrt{5},$$ as seen above.
You stated that $(18 \mathbf{i}) \wedge \mathbf{j}$ "would effectively give" $18 \mathbf{k}$. That is not correct. All that you can say is: $$ \left(18 \mathbf{i}\right) \wedge \mathbf{j} = 18 \left( \mathbf{i} \wedge \mathbf{j} \right). $$ However, in 3D, it is true that $\mathbf{i} \wedge \mathbf{j}$ is dual to $\mathbf{k}$, but the wedge product is a bivector, not a vector (the fact that the wedge of two vectors has the same numerical coefficients as the cross product of those vectors helps lead to that confusion). It is a different beast. This duality may be represented in various ways. In geometric algebra, for example, this duality can be observed by multiplying by the $\mathbb{R}^3$ pseudoscalar $I = \mathbf{i} \mathbf{j} \mathbf{k}$, for example:
$$\begin{aligned}I \mathbf{k}&=\mathbf{i} \mathbf{j} \mathbf{k} \mathbf{k} \\ &=\mathbf{i} \mathbf{j} \\ &=\mathbf{i} \wedge \mathbf{j},\end{aligned}$$
and $$ \mathbf{x} \wedge \mathbf{y} = I \left( { \mathbf{x} \times \mathbf{y} } \right).$$
(this is one way to extract the $\sin\theta$ representation of the wedge product$)
There are other representations of the duality operation that do not require geometric algebra, such as the hodge dual from differential forms, but the end effect is the same, providing a way of mapping a k-form to a N-k form. In general, such duals have to be treated as different algebraic objects.