# (Row?) reduction of a system of equations

For a mathematical problem I am trying to solve, I have (say) $$(N+1)$$ inhomogeneous linear equations in $$M$$ variables of the form

$$\sum_{a=1}^{M}c_{i a}p_{a} = \xi_{i} \quad \text{ where } i=1,\ldots, N$$ $$p_{\alpha} - p_{\beta} = r_0 \qquad \text{ where } \alpha, \beta \in \{1, \ldots, M\} \text{ and } \alpha \neq \beta.$$

(So the first $$N$$ equations have RHS $$\{\xi_i\}_{i=1}^{N}$$ and the $$(N{+}1)^{th}$$ equation has RHS $$r_0$$.)

The integer coefficients $$c_{ia}$$ in each equation add up to zero, i.e.

$$\sum_{a=1}^{M}c_{ia} = 0 \quad \text{ for all } a = 1, \ldots N.$$

For a given range of $$r_0$$, i.e. $$r_{0,min} < r_0 < r_{0,max}$$, I would like to reduce as many equations as possible to the form

$$p_{\alpha_i} - p_{\alpha_j} = f_{ij}(r_0, \xi_1, \ldots)$$

and the rest that cannot be reduced in this way should take the form

$$p_{\alpha_{m_1}} - 2 p_{\alpha_{m_2}} + p_{\alpha_{m_3}} = g(r_0, \xi_1, \ldots)$$

for some pairs $$(\alpha_i, \alpha_j)$$ and some triples $$(\alpha_{m_1}, \alpha_{m_2}, \alpha_{m_3})$$.

Example

Take the matrix

$$M = \left( \begin{array}{cccccccc} -1 & 1 & 0 & 0 & 0 & 1 & -1 & \xi _1 \\ 1 & -1 & 1 & 0 & 0 & 0 & -1 & \xi _2 \\ 0 & 0 & 1 & -1 & 1 & 0 & -1 & \xi _4 \\ 1 & 0 & 0 & 0 & 1 & -1 & -1 & \xi _6 \\ 0 & 0 & 0 & 0 & -1 & 0 & 1 & r_0 \\ \end{array} \right)$$

which is equivalent to the linear system

$$\left( \begin{array}{c} -p_1+p_2+p_6-p_7=\xi _1 \\ p_1-p_2+p_3-p_7=\xi _2 \\ p_3-p_4+p_5-p_7=\xi _4 \\ p_1+p_5-p_6-p_7=\xi _6 \\ p_7-p_5=r_0\\ \end{array} \right)$$

Suppose I run

RowReduce[M]

I get the following output

$$\left( \begin{array}{cccccccc} 1 & 0 & 0 & 0 & 0 & -1 & 0 & \xi _6+r_0\\ 0 & 1 & 0 & 0 & 0 & 0 & -1 & \xi _1+\xi _6+r_0\\ 0 & 0 & 1 & 0 & 0 & 1 & -2 & \xi _1+\xi _2 \\ 0 & 0 & 0 & 1 & 0 & 1 & -2 & \xi _1+\xi _2-\xi _4-r_0\\ 0 & 0 & 0 & 0 & 1 & 0 & -1 & -r_0\\ \end{array} \right)$$

which is equivalent to the linear system

$$\begin{pmatrix}p_1-p_6=\xi _6+r_0,\\p_2-p_7=\xi _1+\xi _6+ r_0,\\p_3+p_6-2 p_7=\xi _1+\xi _2,\\p_4+p_6-2 p_7=\xi _1+\xi _2-\xi _4-r_0,\\p_5-p_7=-r_0\end{pmatrix}$$

It is evident that a further simplification is in fact possible: the third and fourth equations can be subtracted to yield an expression for $$p_3 - p_4$$, namely

$$p_3 - p_4 = \xi_4 + r_0$$

One way of stating the problem is that RowReduce doesn't yield a matrix which corresponds to the fully reduced system of equations (where by reduction I mean the end goal outlined above, which is to reduce the maximal set of equations to the form $$p_{\alpha_i} - p_{\alpha_j} = f_{ij}(r_0, \xi_1, \ldots)$$).

Is there a simple way in Mathematica to automate this or obtain the maximally reduced set of equations?

EDIT: I am aware that instead of RowReduce, I should use LUDecomposition. However, doing so still does not maximally reduce the equation set as desired.

Qred = {{-1, 1, 0, 0, 0, 1, -1}, {1, -1, 1, 0, 0, 0, -1}, {0, 0, 1, -1, 1, 0, -1}, {1, 0, 0, 0, 1, -1, -1}, {0, 0, 0, 0, -1, 0, 1}}; ins = {Subscript[[Xi], 1], Subscript[[Xi], 2], Subscript[[Xi], 4], Subscript[[Xi], 6], r0};

(Qnew = MapThread[Append, {Qred, ins}]) // MatrixForm

(Qnewprime1 = RowReduce[Qnew]) // MatrixForm

(Qnewprime2 = LUDecomposition[Qnew][[1]]) // MatrixForm

• Please add code in a form suitable for copy-and-paste. – Daniel Lichtblau Sep 16 at 15:57
• I'm not sure your "reduction" metric is well-defined: As stated in the documentation of RowReduce, it will return a matrix with an identity-matrix-block on the left, which is a clearly defined target form. The simplification you're proposing destroys that form. The question now is, how exactly should your reduction work? (I'm for example not sure that the simples form according to your definition is reached after the step you are proposing, but this will be quite messy to check) – Lukas Lang Sep 16 at 16:03
• @LukasLang It amounts to a variant of an integer programming problem, wherein one minimizes the number of terms present (an l_0 type of "norm"). – Daniel Lichtblau Sep 16 at 16:44
• @DanielLichtblau I see, thanks for clarifying! – Lukas Lang Sep 16 at 16:50
• @DanielLichtblau, thanks for your replies. I have added some sample code. Thanks also to Lukas Lang for the question. – leastaction Sep 17 at 3:38

I'll show an approach that might or might not work in general but at least seems to have potential. A restriction is that (possibly after preprocessing) all the coefficients of the main variables must be integers. This will be a bit indirect with a few steps, and I do not have the "big picture" insight to clean it up at this point.

The idea is to "exponentiate", that is, treat the integer coefficients as exponents. This will give rise to binomial equations. One can then use some undocumented code, or other tactics I will not go into, to compute what is known as a "toric Groebner basis". A term order is chosen that will tend to make overall degree small, and this will correspond to number of terms in the result upon reverting the exponentiation step.

I will start with the original equations though it would save a step to use the matrix-and-vector input from edit #2 of the original post.

lpolys = {-p1 + p2 + p6 - p7 - e1, p1 - p2 + p3 - p7 - e2,
p3 - p4 + p5 - p7 - e4,
p1 + p5 - p6 - p7 - e6, p7 - p5 - r0};
vars = {p1, p2, p3, p4, p5, p6, p7};

{rhs, lhs} = Normal[CoefficientArrays[lpolys, vars]]

(* Out[1333]= {{-e1, -e2, -e4, -e6, -r0}, {{-1, 1, 0, 0, 0,
1, -1}, {1, -1, 1, 0, 0, 0, -1}, {0, 0, 1, -1, 1, 0, -1}, {1, 0, 0,
0, 1, -1, -1}, {0, 0, 0, 0, -1, 0, 1}}} *)


Due to the way CoefficientArrays handles the right hand side, I had to do a bit of finagling with minus signs in the exponentiation step.

exponpolys =
Numerator[Together[Map[Apply[Times, vars^#] &, lhs] - Exp[-rhs]]]

(* Out[1339]= {p2 p6 - E^e1 p1 p7, p1 p3 - E^e2 p2 p7,
p3 p5 - E^e4 p4 p7, p1 p5 - E^e6 p6 p7, -E^r0 p5 + p7} *)


Here we compute the toric Groebner basis.

tgb = GroebnerBasisToricGroebnerBasis[exponpolys, vars,
CoefficientDomain -> RationalFunctions,
MonomialOrder -> DegreeLexicographic]

(* Out[1353]= {-E^r0 p5 + p7, -p3 + E^(e4 + r0) p4, -p2 +
E^(e1 + e6 + r0) p7, -p1 + E^(e6 + r0) p6, -E^(e4 + r0) p4 p6 +
E^(e1 + e2) p7^2} *)


Had we computed the usual Groebner basis, we would have factors that would cancel on both sides, so it would still be usable. Just a bit more work in the next steps.

We now separate left and right sides of the binomials in the basis.

exponpolys2 = Apply[List, tgb, {1}] /. -1*a_ :> a

(* Out[1350]= {{E^r0 p5, p7}, {p3, E^(e4 + r0) p4}, {p2, E^(e1 + e6 + r0) p7},
{p1, E^(e6 + r0) p6}, {E^(e4 + r0) p4 p6, E^(e1 + e2) p7^2}} *)


Now we revert the exponentiation step, in effect going from binomials to linear polynomials.

logs = Map[PowerExpand[Log[#]] &, exponpolys2, {2}] /. Log[a_] :> a

(* Out[1351]= {{p5 + r0, p7}, {p3, e4 + p4 + r0}, {p2, e1 + e6 + p7 + r0},
{p1, e6 + p6 + r0}, {e4 + p4 + p6 + r0, e1 + e2 + 2 p7}} *)


Last step is to reform the two sides in each p[air as an equation.

newpolys = Apply[Subtract, logs, {1}]

(* Out[1352]= {p5 - p7 + r0, -e4 + p3 - p4 - r0, -e1 - e6 + p2 - p7 - r0,
-e6 + p1 - p6 - r0, -e1 - e2 + e4 + p4 + p6 - 2 p7 + r0} *)
`

At this point one obtains mostly linear equations with just two variables. This was due to the minimizing of total degree when these were in "exponentiated" binomial form.

I would think there are simplifications that could be made to this process. But this is the basic idea, in elaborated form.

• Thanks for your detailed reply, Daniel. The objective of the <black box> is to take those $(N+1)$ linear equations and express a certain number $K$ of them into the form $p_i - 2 p_j + p_k = f(r_0, \xi_1, \xi_2, \ldots)$ and the rest $(N-K+1)$ into the form $p_\alpha - p_\beta = g(r_0, \xi_1, \xi_2, \ldots)$. If you wish, the $(N+1)^{th}$ input equation is a "pivot". Indeed these equations arise from a toric geometry even in the my application. And the number $K$ is characteristic of the underlying geometry. Why do you think this approach is not expected to work in general? – leastaction Sep 17 at 17:20
• I think your point is that after the last step of transposing the $r_0$ and $e_i$ terms in the output of your code to the right-hand side, I would get the desired form. Is there a reference for this procedure? I'll be happy to discuss it more particularly because these equations originate as moment map conditions of the resolution of a toric singularity. – leastaction Sep 17 at 17:23
• (1) I should restate matters. I do not know that this will fail in general cases. I'm just not certain it will always give the desired result, if, say, there are other equations in the ideal of total degree 2. But maybe that is ruled out for some reason. – Daniel Lichtblau Sep 17 at 18:31
• (2) As for reference in the context of this problem, I do not know. I kinda rolled my own on this one. A reference for using this basic approach in the context of integer linear programming would be chapter 2 section 8 of the book "An Introduction to Gro"bner Bases" by Adams and Loustaunau. They in turn base it on a 1991 article by Conti and Traverso. There is also this Mathematica notebook by my colleague Devendra Kapadia in our cyber-library. – Daniel Lichtblau Sep 17 at 18:35
• – Daniel Lichtblau Sep 17 at 18:38