Convert Multivariate Polynomial to Matrix Representation

Question

Let's say I have a multivariate polynomial:

import sympy

x_1, x_2, x_3, x_4 = sympy.symbols('x_1 x_2 x_3 x_4')

expressions = [
           -x_1*x_1-x_2*x_2+x_1*x_2,
           x_2*x_2-x_1*x_2,
           x_1*x_1-x_1*x_2,
           x_1*x_2,
           x_1*x_1-x_1*x_3,
           x_1*x_3,
           x_3*x_3-x_2*x_3,
           -x_2*x_2-x_4*x_4+x_2*x_4,
           x_2*x_2-x_2*x_3,
           x_2*x_3,
           -x_3*x_3-x_4*x_4+x_3*x_4,
           x_3*x_4,
          ]

model = sympy.Poly(sympy.Add(*expressions))
model

# Poly(x_1**2 - x_2*x_3 + x_2*x_4 + 2*x_3*x_4 - 2*x_4**2, x_1, x_2, x_3, x_4, domain='ZZ')

Notice that the variables are [x_1, x_2, x_3, x_4] and so it is possible to represent the coefficients of the polynomial as a 4x4 square matrix where the coefficient of the squared terms (i.e., x_i*x_i) are the diagonal terms along the matrix and the off-diagonal terms depend on the coefficients of x_i*x_j :

[[1, 0,  0,  0],
 [0, 0, -1,  1],
 [0, 0,  0,  2],
 [0, 0,  0, -2]
]

Starting with the sympy polynomial, is it possible to extract the coefficients and construct the corresponding sympy matrix as shown above for ANY polynomial with variables [x_1, x_2, ..., x_N]?

At the end of the day, I'm really hoping to obtain the final matrix as a numpy array so that it can be used for additional computation outside of sympy.

Answer 1

As long as all monomials are of degree 2 (e.i. a_ij * x_i * x_j), this should work:

import sympy

x_1, x_2, x_3, x_4 = sympy.symbols('x_1 x_2 x_3 x_4')

expressions = [
           -x_1*x_1-x_2*x_2+x_1*x_2,
           x_2*x_2-x_1*x_2,
           x_1*x_1-x_1*x_2,
           x_1*x_2,
           x_1*x_1-x_1*x_3,
           x_1*x_3,
           x_3*x_3-x_2*x_3,
           -x_2*x_2-x_4*x_4+x_2*x_4,
           x_2*x_2-x_2*x_3,
           x_2*x_3,
           -x_3*x_3-x_4*x_4+x_3*x_4,
           x_3*x_4,
          ]

model = sympy.Add(*expressions)

Q = sympy.Rational(1, 2) * sympy.hessian(model, (x_1, x_2, x_3, x_4))
sympy.pprint(Q)

It gives:

⎡1   0     0     0 ⎤
⎢                  ⎥
⎢0   0    -1/2  1/2⎥
⎢                  ⎥
⎢0  -1/2   0     1 ⎥
⎢                  ⎥
⎣0  1/2    1    -2 ⎦

Or, if you need an upper triangular matrix:

import numpy as np

h = np.array(sympy.hessian(model, (x_1, x_2, x_3, x_4)), dtype=float)
Q = np.triu(h) - 0.5 * np.diag(np.diag(h))
print(Q)

which gives:

[[ 1.  0.  0.  0.]
 [ 0.  0. -1.  1.]
 [ 0.  0.  0.  2.]
 [ 0.  0.  0. -2.]]