I want to evaluate a sum like $\sum\limits_{j > i}^n q_i^k\, q_j^k$, where $n$ is the length of the vector $q^k$. The vector $q^k$ is, for example, $[q_1, q_2, q_3, q_4]$ and $i$ and $j$ correspond to the index of the vector elements.
This was my attempt; I got errors :/
vec = Table[q[i], {i, 4}]
(* {q[1], q[2], q[3], q[4]} *)
With[{l = Length[vec]},
Sum[vec[[i]] vec[[j]], {i, 1, l}, {j, i + 1, l}]
]
(* q[1] q[2] + q[1] q[3] + q[2] q[3] + q[1] q[4] + q[2] q[4] + q[3] q[4] *)
Plus @@ Times @@@ Subsets[vec, {2}]
(* q[1] q[2] + q[1] q[3] + q[2] q[3] + q[1] q[4] + q[2] q[4] + q[3] q[4] *)
Update
For binary vectors (see OP's comment), a fast approach would be:
vec = RandomInteger[{0, 1}, 10^6];
Binomial[Total[vec], 2] // RepeatedTiming
(* {0.0010, 124964252556} *)
(* comparison of the result with Jim Baldwin's answer *)
%[[2]] === (Total[vec^k]^2 - Total[vec^(2 k)])/2
(* True *)
SymmetricPolynomial is the built-in way to get the sums you want:
vec = Table[q[i], {i, 4}];
SymmetricPolynomial[2, vec^k]
(* q[1]^k q[2]^k+q[1]^k q[3]^k+q[2]^k q[3]^k+q[1]^k q[4]^k+q[2]^k q[4]^k+q[3]^k q[4]^k *)
(But @Xavier 's way will help you understand Mathematica better.)
Update
To follow up on @ciao 's comment, if the values are numeric, then there are definitely better ways. For your example, the following works:
n = 4000;
k = 2;
vec = Abs[RandomVariate[NormalDistribution[1, 0.2], n]];
Timing[(Total[vec^k]^2 - Total[vec^(2 k)])/2]
(* {0.`,8.583920493183017`*^6} *)
Timing[SymmetricPolynomial[2, vec^k]]
(* {2.6364169`,8.583920493183037`*^6} *)
Note the difference in timing. The point is that one can determine such sums as you have using just using the sums of the powers of the vector elements.
