I need to organize a sum of terms in a pseudo-boolean function into a matrix!

Question

I have a long 2-local summation of terms that I want to organize into a matrix such that the coefficients of the 1-local terms (a term with just a scalar coefficient and a single variable) are on the diagonal and the 2 local coefficients (term with scalar and two variables) are in the upper triangular portion of the matrix. This can be in any configuration, which should make the problem much simpler. Here is an portion of the kind of functions I am working with:

-q[1] r[1] - q[2] r[1] - q[3] r[1] - q[4] r[1] - q[1] r[2] - 
 q[2] r[2] - q[3] r[2] - q[4] r[2] - q[1] r[3] - q[2] r[3] - 
 q[3] r[3] - q[4] r[3] - q[1] r[4] - q[2] r[4] - q[3] r[4] - 
 q[4] r[4] - q[1] r[5] - q[2] r[5] - q[3] r[5] - q[4] r[5] - 
 q[1] r[6] - q[2] r[6] - q[3] r[6] - q[4] r[6] - q[1] r[7] - 
 q[2] r[7] - q[3] r[7] - q[4] r[7] - r[1] s[1] - r[2] s[1] - 
 r[3] s[1] - r[4] s[1] - r[5] s[1] - r[6] s[1] - r[7] s[1] - 
 r[8] s[1] - r[1] s[2] - r[2] s[2] - r[3] s[2] - r[4] s[2] - 
 r[5] s[2] - r[6] s[2] - r[7] s[2] - r[8] s[2] - r[1] s[3] - 
 r[2] s[3] - r[3] s[3] - r[4] s[3] - r[5] s[3] - r[6] s[3] - 
 r[7] s[3] - r[8] s[3] - r[1] s[4] - r[2] s[4] - r[3] s[4] - 
 r[4] s[4] - r[5] s[4] - r[6] s[4] - r[7] s[4] - r[8] s[4] - 
 r[1] s[5] - r[2] s[5] - r[3] s[5] - r[4] s[5] - r[5] s[5] - 
 r[6] s[5] - r[7] s[5] - r[8] s[5] - r[1] s[6] - r[2] s[6] - 
 r[3] s[6] - r[4] s[6] - r[5] s[6] - r[6] s[6] - r[7] s[6] - 
 r[8] s[6] - r[1] s[7] - r[2] s[7] - r[3] s[7] - r[4] s[7] - 
 r[5] s[7] - r[6] s[7] - r[7] s[7] - r[8] s[7] - r[1] s[8] - 
 r[2] s[8] - r[3] s[8] - r[4] s[8] - r[5] s[8] - r[6] s[8] - 
 r[7] s[8] - r[8] s[8] - r[1] s[9] - r[2] s[9] - r[3] s[9] - 
 r[4] s[9] - r[5] s[9] - r[6] s[9] - r[7] s[9] - r[8] s[9] - 
 r[1] s[10] - r[2] s[10] - r[3] s[10] - r[4] s[10] - r[5] s[10] - 
 r[6] s[10] - r[7] s[10] - r[8] s[10] - r[1] s[11] - r[2] s[11] - 
 r[3] s[11] - r[4] s[11] - r[5] s[11] - r[6] s[11] - r[7] s[11] - 
 r[8] s[11] - r[1] s[12] - r[2] s[12] - r[3] s[12] - r[4] s[12] - 
 r[5] s[12] - r[6] s[12] - r[7] s[12] - r[8] s[12] - r[1] s[13] - 
 r[2] s[13] - r[3] s[13] - r[4] s[13] - r[5] s[13] - r[6] s[13] - 
 r[7] s[13] - r[8] s[13] - r[1] s[14] - r[2] s[14] - r[3] s[14] - 
 r[4] s[14] - r[5] s[14] - r[6] s[14] - r[7] s[14] - r[8] s[14] - 
 r[1] s[15] - r[2] s[15] - r[3] s[15] - r[4] s[15] - r[5] s[15] - 
 r[6] s[15] - r[7] s[15] - r[8] s[15] + 6s[15]

In this example, there is only 1 1-local term, (6s[15]) so there would only be a 6 on the diagonal. I'm new to Mathematica, so I was wondering if there was a nice/elegant solution, besides just parsing the function and throwing the coefficients into 2 arrays (one for 1-local, one for two) and smooshing it into a matrix.

Answer 1

This should do what you want:

polynomial = -q[1] r[1] - q[2] r[1] - q[3] r[1] - q[4] r[1] - q[1] r[2] - q[2] r[2] - q[3] r[2] - q[4] r[2] - q[1] r[3] - q[2] r[3] - q[3] r[3] - q[4] r[3] - q[1] r[4] - q[2] r[4] - q[3] r[4] - q[4] r[4] - q[1] r[5] - q[2] r[5] - q[3] r[5] - q[4] r[5] - q[1] r[6] - q[2] r[6] - q[3] r[6] - q[4] r[6] - q[1] r[7] - q[2] r[7] - q[3] r[7] - q[4] r[7] - r[1] s[1] - r[2] s[1] - r[3] s[1] - r[4] s[1] - r[5] s[1] - r[6] s[1] - r[7] s[1] - r[8] s[1] - r[1] s[2] - r[2] s[2] - r[3] s[2] - r[4] s[2] - r[5] s[2] - r[6] s[2] - r[7] s[2] - r[8] s[2] - r[1] s[3] - r[2] s[3] - r[3] s[3] - r[4] s[3] - r[5] s[3] - r[6] s[3] - r[7] s[3] - r[8] s[3] - r[1] s[4] - r[2] s[4] - r[3] s[4] - r[4] s[4] - r[5] s[4] - r[6] s[4] - r[7] s[4] - r[8] s[4] - r[1] s[5] - r[2] s[5] - r[3] s[5] - r[4] s[5] - r[5] s[5] - r[6] s[5] - r[7] s[5] - r[8] s[5] - r[1] s[6] - r[2] s[6] - r[3] s[6] - r[4] s[6] - r[5] s[6] - r[6] s[6] - r[7] s[6] - r[8] s[6] - r[1] s[7] - r[2] s[7] - r[3] s[7] - r[4] s[7] - r[5] s[7] - r[6] s[7] - r[7] s[7] - r[8] s[7] - r[1] s[8] - r[2] s[8] - r[3] s[8] - r[4] s[8] - r[5] s[8] - r[6] s[8] - r[7] s[8] - r[8] s[8] - r[1] s[9] - r[2] s[9] - r[3] s[9] - r[4] s[9] - r[5] s[9] - r[6] s[9] - r[7] s[9] - r[8] s[9] - r[1] s[10] - r[2] s[10] - r[3] s[10] - r[4] s[10] - r[5] s[10] - r[6] s[10] - r[7] s[10] - r[8] s[10] - r[1] s[11] - r[2] s[11] - r[3] s[11] - r[4] s[11] - r[5] s[11] - r[6] s[11] - r[7] s[11] - r[8] s[11] - r[1] s[12] - r[2] s[12] - r[3] s[12] - r[4] s[12] - r[5] s[12] - r[6] s[12] - r[7] s[12] - r[8] s[12] - r[1] s[13] - r[2] s[13] - r[3] s[13] - r[4] s[13] - r[5] s[13] -  r[6] s[13] - r[7] s[13] - r[8] s[13] - r[1] s[14] - r[2] s[14] - r[3] s[14] - r[4] s[14] - r[5] s[14] - r[6] s[14] - r[7] s[14] - r[8] s[14] - r[1] s[15] - r[2] s[15] - r[3] s[15] - r[4] s[15] - r[5] s[15] - r[6] s[15] - r[7] s[15] - r[8] s[15] + 6 s[15]

variables = Sort@Variables[polynomial]
(* {q[1], q[2], q[3], q[4], r[1], r[2], r[3], r[4], r[5], r[6], r[7], r[8], s[1], s[2], s[3], s[4], s[5], s[6], s[7], s[8], s[9], s[10], s[11], s[12], s[13], s[14], s[15]} *)

arrays = CoefficientArrays[polynomial, variables];

(matrix = DiagonalMatrix@arrays[[2]] + arrays[[3]]) // MatrixForm

We can recover the original polynomial by dotting with the vector of variables, except the diagonal terms will have an extra factor:

variables.matrix.variables - polynomial // Expand
(* -6 s[15] + 6 s[15]^2 *)