(Row?) reduction of a system of equations

Question

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.

EDIT #2: Code added below

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

Answer 1

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 = GroebnerBasis`ToricGroebnerBasis[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.