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}\$.

Python3, 216 bytes

def f(v,V):
 t=[[0],[-3,0],[2,-1,0],[-5,-6,-7,0],[4,-7,6,-1,0],[7,4,-5,-2,3,0],[-6,5,4,-3,-2,1,0]]
 D={}
 for x in range(7):
  for y in range(7):D[u]=D.get(u:=t[x][y]if y<len(t[x])else -t[y][x],0)+v[x]*V[y]
 return D

Try it online!

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

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.

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.

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.