Extract linear/quadratic/nonlinear terms from nonlinear expressions

[ZedongPeng]

Summary

I want to extract linear/quadratic/nonlinear terms from nonlinear expressions.

Description

I first tried expressions.args and check the polynomial_degree() of each subexpressions. However, if the subexpression is a constant. An error will be reported since int and float have no polynomial_degree().

Then I tried first transforming expressions into repn and then obtaining the linear/quadratic/nonlinear terms. I got stuck in the second step. The nonlinear_expr is stored in the repn, but linear_expr and quadratic_expr are not stored.

I know that linear_coefs, linear_vars, quadratic_coefs and quadratic_vars are provided. linear_expr and quadratic_expr can be obtained by SUMPRODUCT.

Is there an existing function to do this? I think to_expression might be the best candidate. Can we add some args, like linear_only, quadratic_only for this function?

[jsiirola]

Yes. See if generate_standard_repn will do what you want. You will want to do some testing, though. It is designed to gather linear and quadratic terms for the LP writer, and will "throw up its hands" and dump nonlinear things into a separate subexpression. That said, we don't generally use it for extracting quadratic terms from a general nonlinear expression, so you will want to check that it doesn't accidentally lump some quadratic terms into the nonlinear subexpression (It shouldn't, but that is not somethign that we test).

[michaelbynum]

I can confirm that generate_standard_repn will not always extract linear or quadratic terms from nonlinear expressions:

In [1]: import pyomo.environ as pe

In [2]: m = pe.ConcreteModel()

In [3]: m.x = pe.Var()

In [4]: m.y = pe.Var()

In [5]: from pyomo.repn.standard_repn import generate_standard_repn

In [6]: e = (m.x + m.y + pe.log(m.x)) * m.x

In [8]: print(generate_standard_repn(e))

constant:       0
linear vars:    []
linear var ids: []
linear coef:    []
quadratic vars:    []
quadratic var ids: []
quadratic coef:    []
nonlinear expr:
(x + y + log(x))*x
nonlinear vars: ['x', 'y']



In [9]: e = (m.x + m.y + pe.log(m.x) + 2) * m.x

In [10]: print(generate_standard_repn(e))

constant:       0
linear vars:    []
linear var ids: []
linear coef:    []
quadratic vars:    []
quadratic var ids: []
quadratic coef:    []
nonlinear expr:
(x + y + log(x) + 2)*x
nonlinear vars: ['x', 'y']

[jsiirola]

OK. I am not too surprised. I have been working on new versions of the expression compilers (replacements for generate_standard_repn) that walk the expressions from the "bottom up" instead of "top down". They should help make the routines more robust. Unfortunately, at the moment both those walkers (sitting on the nl-writer-cache-expressions and templated-lp-models branches of my fork) only extract the linear portion of the models. I have a reasonable idea how to extend them to gathering quadratic terms (for the LP writers), but that is something I probably won't get to for several months.

[ZedongPeng]

Thanks. @jsiirola @michaelbynum

[ZedongPeng]

I have one question here. If the coefficient of linear terms is native_numeric_types, why would the accumulation become so complicated? like expr -= - c*v.

https://github.com/Pyomo/pyomo/blob/5683dc3d46e6c4a3a675624f23b0c14d6578be1a/pyomo/repn/standard_repn.py#L180-L188

[bernalde]

I would argue that this issue could be kept open as a placeholder for the walker that Jhon is putting together. It would certainly make Pyomo's life easier in the upcoming world where more and more (MI)LP solvers are supporting quadratic constraints.

[jsiirola]

@bernalde, agreed. Note that the LP writer and generate_standard_repn work fine for quadratic expressions. The limitation @michaelbynum pointed out (which I am not surprised at) is that generate_standard_repn will not reliability identify / extract quadratic terms from general nonlinear expressions.

[michaelbynum]

Well said, @jsiirola. Thank you for clarifying that.

[jsiirola]

Once #2997 is merged, we can close this. The new representation walker is in and should resolve the issues that @ZedongPeng and @michaelbynum reported. There is one thing to note: the default implementation of the QuadraticRepnVisitor does not expand nonlinear expressions (this expansion can lead to significantly larger expressions and is unnecessary for the normal use of the QuadraticRepnVisitor (i.e. in the LP file writer, where general nonlinear expressions will eventually just lead to an error):

>>> from pyomo.repn.quadratic import QuadraticRepnVisitor
>>> walker = QuadraticRepnVisitor({}, {}, {})
>>> import pyomo.environ as pe
>>> m = pe.ConcreteModel()
>>> m.x = pe.Var()
>>> m.y = pe.Var()
>>> e1 = (m.x + m.y + pe.log(m.x)) * m.x
>>> e2 = (m.x + m.y + pe.log(m.x) + 2) * m.x
>>> print(walker.walk_expression(e1))
QuadraticRepn(mult=1, const=0, linear={}, quadratic=None, nonlinear=(log(x) + (x + y))*x)
>>> print(walker.walk_expression(e2))
QuadraticRepn(mult=1, const=0, linear={}, quadratic=None, nonlinear=(log(x) + (x + y) + 2)*x)

However, if you enable nonlinear term expansion, then you will get the output you want:

>>> walker.expand_nonlinear_products= True
>>> print(walker.walk_expression(e1))
QuadraticRepn(mult=1, const=0, linear={}, quadratic={(140448631786704, 140448631786704): 1, (140448631786704, 140448631946320): 1}, nonlinear=log(x)*x)
>>> print(walker.walk_expression(e2))
QuadraticRepn(mult=1, const=0, linear={140448631786704: 2}, quadratic={(140448631786704, 140448631786704): 1, (140448631786704, 140448631946320): 1}, nonlinear=log(x)*x)