Jelly, 28 bytes

ṙJ$⁽Ç¢œ?N5ẒƇ¤¦_W{JœṖ"$ZṚẎŒḌḋ

A dyadic Link that accepts \$\mathbf v_1\$ on the left and \$\mathbf v_2\$ on the right and yields \$\mathbf v_1\times\mathbf v_2\$

Try it online!

• See the orthonormal multiplication table.
• See two random vectors and two random scalars satisfying \$a\mathbf{v}_1\times b\mathbf{v}_2 = ab\left(\mathbf{v}_1\times\mathbf{v}_2\right)\$ here.
• See three random vectors satisfying \$\left(\mathbf v_1+\mathbf v_2\right)\times\mathbf v_3 = \left(\mathbf v_1\times\mathbf v_3\right)+\left(\mathbf v_2\times\mathbf v_3\right)\$ here.

How?

Uses this multiplication table of orthonormal basis vectors:

ṙJ$⁽Ç¢œ?N5ẒƇ¤¦_W{JœṖ"$ZṚẎŒḌḋ - Link: V1; V2  e.g. abcdefg; tuvwxyz
ṙJ$                          - rotate {V1} left by its 1-indices
                               -> bcdefga cdefgab defgabc efgabcd fgabcde gabcdef abcdefg
   ⁽Ç¢œ?                     - 4502nd lexicographic permutation of {that}
                               -> abcdefg cdefgab efgabcd fgabcde bcdefga gabcdef defgabc
        N5ẒƇ¤¦               - negate values in the 2nd, 3rd, and 5th lists (primes up to 5)
              _W{            - replace the 1st with zeros (subtract [V1])
                 JœṖ"$       - partition each before its 1-index
                      ZṚẎ    - transpose, reverse and tighten
                         ŒḌ  - form a matrix from these diagonals
                               ->  0 -d -g  b -f  e  c
                                   d  0 -e -a  c -g  f
                                   g  e  0 -f -b  d -a
                                  -b  a  f  0 -g -c  e
                                   f -c  b  g  0 -a -d
                                  -e  g -d  c  a  0 -b
                                  -c -f  a -e  d  b  0
                           ḋ - dot-product {V2} (vectorsies)
                                 -> [  0 - ud - vg + wb - xf + ye + zc,
                                      td +  0 - ve - wa + xc - yg + zf,
                                      tg + ue +  0 - wf - xb + yd - za,
                                     -tb + ua + vf +  0 - xg - yc + ze,
                                      tf - uc + vb + wg +  0 - ya - zd,
                                     -te + ug - vd + wc + xa +  0 - zb,
                                     -tc - uf + va - we + xd + yb +  0,
                                    ]

Python 2, 165 146 bytes

def g(x,y):
 x*=2;y*=2;z=[0]*10;p,q,r=0,1,3
 for a in range(7)*3:p,q,r=q,r,p;z[a+p]+=x[a+q]*y[a+r]-x[a+r]*y[a+q]
 exec'print sum(z[p::7]);p+=1;'*7

Attempt This Online!

Changed my approach completely.

This is now a port of Seven-dimensional cross product.

Prints output to stdout.

JavaScript (ES6), 83 bytes

Expects (u)(v).

u=>v=>[0,1,2,4].reduce((p,i)=>p.map(n=>i&&n+u[i%=7]*v[j%=7]-u[j++]*v[i++],j=i*3),u)

Try it online!

Commented

u =>              // u[] = first vector
v =>              // v[] = second vector
[0, 1, 2, 4]      // lookup list
.reduce((p, i) => // for each value i with p[] as the accumulator:
  p.map(n =>      //   for each value n in p[]:
    i &&          //     force to zero if i = 0
    n +           //     otherwise, take the current value n
    u[i %= 7] *   //     and add u[i] * v[j] - u[j] * v[i]
    v[j %= 7] -   //     where both i and j are reduced modulo 7
    u[j++] *      //     and incremented afterwards
    v[i++],       //
    j = i * 3     //     start with j = i * 3, which gives (modulo 7):
                  //     (i,j) = (1,3), (2,6), (4,5)
  ),              //   end of map()
  u               //   start with p[] = u[] (but the content of p[]
                  //   is cleared on the first iteration)
)                 // end of reduce()

Wolfram Language (Mathematica), 62 bytes

s&@Do[s[[t]]+=#[[t=Mod[i+{0,1,3},7,1]]]#2[[t]],{i,s=0#;7}]&

-1 bytes thanks to @att.

Try it online!

 (\[Cross], unicode U+F4A0) is a special character in Mathematica that represents cross product.

First initialize the result \$\mathbf{s}\$ to a zero vector of length 7. And then:

• \$(s_1, s_2, s_4) \leftarrow (s_1, s_2, s_4) + (a_1, a_2, a_4) \times (b_1, b_2, b_4)\$ (where \$\times\$ is the 3D cross product);
• \$(s_2, s_3, s_5) \leftarrow (s_2, s_3, s_5) + (a_2, a_3, a_5) \times (b_2, b_3, b_5)\$;
• \$\ldots\ldots\$
• \$(s_7, s_1, s_3) \leftarrow (s_7, s_1, s_3) + (a_7, a_1, a_3) \times (b_7, b_1, b_3)\$.

Now \$\mathbf{s}\$ becomes the 7D cross product of \$\mathbf{a}\$ and \$\mathbf{b}\$.

Charcoal, 63 35 bytes

ＩＥ⁷ΣＥＥ⁷⊗﹪⁻μι⁷∧λ×⎇＆λ⊖λ§θ⁻ι뱧θ⁺μλ§ημ

Try it online! Link is to verbose version of code. No test suite because that requires renaming the variables. Explanation: Uses the C(+) octonion multiplication table that I found here (warning: no SSL available).

  ⁷                                 Literal integer `7`
 Ｅ                                  Map over implicit range
      ⁷                             Literal integer `7`
     Ｅ                              Map over implicit range
          μ                         Inner index
         ⁻                          Subtract
           ι                        Outer value
        ﹪                           Modulo
            ⁷                       Literal integer `7`
       ⊗                            Doubled
    Ｅ                               Map over values
              λ                     Inner value
             ∧                      Logical And
                  λ                 Inner value
                 ＆                  Bitwise And
                    λ               Inner value
                   ⊖                Decremented
                ⎇                   If true then
                      θ             First input
                     §              Cyclically indexed by
                        ι           Outer value
                       ⁻            Subtract
                         λ          Inner value
                            θ       Else first input
                           §        Cyclically indexed by
                              μ     Inner index
                             ⁺      Plus
                               λ    Inner value
                          ±         Negated
               ×                    Multiplied by
                                 η  Second input
                                §   Indexed by
                                  μ Inner index
   Σ                                Take the sum
Ｉ                                   Cast to string
                                    Implicitly print

PARI/GP, 74 bytes

f(a,b)=[sum(n=0,2,a[u=(i+2^n)%7+1]*b[v=(i+3<<n)%7+1]-b[u]*a[v])|i<-[0..6]]

Attempt This Online!

A port of Arnauld's JavaScript solution.

PARI/GP, 97 bytes

f(a,b)=Vecrev(sum(t=0,6,matdet([x^t,x^u=t++%7,x^v=(t+2)%7;a[t],a[u++],a[v++];b[t],b[u],b[v]])),7)

Attempt This Online!

Using the algorithm in my Mathematica answer.

Maple,  187  176 bytes

(s,t)->eval(s^+.Matrix([[],[-e3],[e2,-e1],[-e5,-e6,-e7],[e4,-e7,e6,-e1],[e7,e4,-e5,-e2,e3],[-e6,e5,e4,-e3,-e2,e1,0]],shape=antisymmetric).t,{seq(e||i=Vector(7,{i=1}),i=1..7)});

If the two vectors are s and t, and Q is the multiplication table as a matrix, then the result in the form a*e1+b*e2+... is s^T.Q.t. Then substitute e1, e2 etc with the unit vectors [1,0,0,0,0,0,0]^T, [0,1,0,0,0,0,0]^T etc. to get the final result.