# How to write a more concise RowReduce function that can deal with symbolic matrix?

I asked a similar question in this post before and understand the reason why RowReduce function can't reduce symbolic matrix.

I wrote an RowReduce function that can handle symbolic matrices:

exchanges[v1_, v2_] :=
Module[{t},
t = InversePermutation[PermutationList[FindPermutation[v1, v2]]];
If[t != {},
Select[MapIndexed[First[#2] -> #1 &,
LinearAlgebraLAPACKPermutationToPivot[t]], Apply[Unequal]], {}]]

matrixWrapOperation[mat_, rule_] :=
Module[{B = mat},
If[rule != {},
Do[B[[{Keys[rule[[n]]], Values[rule[[n]]]}, ;;]] =
B[[{Values[rule[[n]]], Keys[rule[[n]]]}, ;;]];
Print["Swap the elements in row", Keys[rule[[n]]], "and row ",
Values[rule[[n]]], "->", MatrixForm@B], {n, 1, Length[rule]}]; B,
B]]

rref[A_?MatrixQ] :=
Module[{sysa, sysm, sysn, sysi, sysj, sysk, sysl, sysL, sysB},
{sysm, sysn} = Dimensions[A]; sysB = A;
Print[MatrixForm@A];

(*Below is the preprocessing sort\[DownArrow]*)
Print["Start preprocessing"];
sysB = matrixWrapOperation[A,
exchanges[
A[[;; , 1]], (SortBy[A /. {0 -> Infinity},
First] /. {Infinity -> 0})[[All, 1]]]];
Print["End of preprocessing"];
(*The above is the preprocessing sort\[UpArrow]*)

sysi = 1(*sysi represents the row number*);
sysk = 1(*sysk represents the column number*); While[sysi <= sysm,
While[sysk < sysn &&
sysB[[sysi ;;, sysk]] == ConstantArray[0, sysm - sysi + 1], sysk++;
If[sysk == sysn && sysB[[sysi ;;, sysk]] == 0, Return[sysB];
Goto[end]]];

(*For[sysl=sysi,sysl\[LessEqual]sysm,sysl++,If[(sysa=sysB[[sysl,
sysk]])\[NotEqual]0&&sysk<sysn&&sysB[[
sysl,(sysk+1);;]]\[Equal]ConstantArray[0,sysn-sysk],sysB[[sysl,
sysk]]=1;Print["Divide all elements in row",sysl," by ",sysa,"->",
MatrixForm@sysB];]];*)

sysa = DeleteCases[sysB[[sysi ;;, sysk]], 0];
If[Length[sysa] == 0, sysL = 1, sysL = First[sysa]];

sysj = sysi;
While[sysj < sysm && sysB[[sysj, sysk]] =!= sysL, sysj++;];

If[sysi != sysj,
If[sysB[[sysj, sysk]] =!= 0,
sysB[[{sysi, sysj}, ;;]] = sysB[[{sysj, sysi}, ;;]];
Print["Swap the elements in row", sysi, "and row ", sysj, "->",
MatrixForm@sysB], Return[sysB]]
];

For[sysl = sysi, sysl <= sysm, sysl++,
If[(sysa = sysB[[sysl, sysk]]) =!= 0 && sysl > sysi,
sysB[[sysl, ;;]] =
sysB[[sysl, ;;]] -
sysB[[sysl, sysk]]/sysB[[sysi, sysk]] sysB[[sysi, ;;]];

Print["The ", sysl, " row of elements plus ", -(sysa/sysL),
" times the ", sysi, " row of elements",
"-> " MatrixForm@sysB]]];
sysi++(*Outermost row loop count*);]; Label[end];

Return[sysB];]


Examples for testing:

A = RandomInteger[{0, 10}, {4, 4}];
rref[A]
rref[{{1, 9, 7, 2}, {0, 2, -(9/2),
2}, {0, -(11/2), -(11/2), -t},
{s, 10, 9, 7}, {1, 2, 4, 3}}]

rref[{{8, 0, 0, 1, 3}, {8, 0, 0, 0, 3},
{6, 0, 0, 2, 4}, {7, 0, 0, 8, 9},
{10, 0, 0, 1, 5}}]

rref[{{1, -1, -1}, {2, a, 1}, {-1, 1, a}}]

rref[{{1, -1, -1, 2, 2}, {2, a, 1, 1, a}, {-1, 1, a, -a - 1, -2}}]


But the above code is too cumbersome. Is there a more concise way to realize the function of rref function?

• One concise possibility is, well, 'RowReduce'. You should therefore explain why it is not suitable. Aug 8, 2020 at 16:21
• @DanielLichtblau For example: the row simplest form of matrix $\left(\begin{array}{ccc} 1 & a & 2 \\ 0 & 1 & 1 \\ -1 & 1 & 1 \end{array}\right)$ should be $\left(\begin{array}{ccc} 1 & -1 & -1 \\ 0 & 1 & 1 \\ 0 & a-2 & 0 \end{array}\right)$ instead of $\left(\begin{array}{lll} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array}\right)$. RowReduce[{{1, a, 2}, {0, 1, 1}, {-1, 1, 1}}] rref[{{1, a, 2}, {0, 1, 1}, {-1, 1, 1}}] Aug 8, 2020 at 22:24
• Perhaps the point to emphasize is that RowReduce seems to make the generic assumption that any expression that is not literally zero is nonzero. Aug 9, 2020 at 11:29
• @Montevideo Perhaps you should refer to your previous question (mathematica.stackexchange.com/q/227929/43522) to clarify the problem. As it stands, this question doesn't explain the problem very well. Aug 11, 2020 at 12:10
• @SjoerdSmit Thank you very much for your advice. I have updated the question. Aug 11, 2020 at 22:18

Perhaps this, but I'm unsure how safe it is to treat a - 2 as zero:

mat = {{1, a, 2}, {0, 1, 1}, {-1, 1, 1}};
RowReduce[mat,
ZeroTest -> (Quiet[Length@Solve[# == 0, Reals] > 0] &)]


$$\left( \begin{array}{ccc} 1 & 0 & 2-a \\ 0 & 1 & 1 \\ 0 & 0 & 2-a \\ \end{array} \right)$$

There's also

UpperTriangularize@First@LUDecomposition@mat


$$\left( \begin{array}{ccc} 1 & a & 2 \\ 0 & 1 & 1 \\ 0 & 0 & 2-a \\ \end{array} \right)$$

The columns are not reduced. Also, I suspect LUDecomposition would divide by a - 2 if needed.

• Thank you very much, but I hope your method can handle the following situation:mat = {{1, -1, -1, 2, 2}, {2, a, 1, 1, a}, {-1, 1, a, -a - 1, -2}}; RowReduce[mat, ZeroTest -> (Quiet[Length@Solve[# == 0, Reals] > 0] &)]. Reference answer: $\left(\begin{array}{ccccc} 1 & -1 & -1 & 2 & 2 \\ 2 & a & 1 & 1 & a \\ -1 & 1 & a & -a-1 & -2 \end{array}\right) \rightarrow\left(\begin{array}{ccccc} 1 & -1 & -1 & 2 & 2 \\ 0 & a+2 & 3 & -3 & a-4 \\ 0 & 0 & a-1 & 1-a & 0 \end{array}\right)$ Aug 10, 2020 at 1:23
• I hope it can handle it, too, but there are doubts. Nonetheless: Take[UpperTriangularize@First@LUDecomposition@ArrayReshape[mat, {5, 5}], 3] Aug 10, 2020 at 2:06
• Thank you very much for your reply, your code works very well. Aug 10, 2020 at 3:26
• @Montevideo I'm not sure exactly what you want. I'm pretty sure the built-in functions as I've shown them, do something close, but not exactly the same. If your requirements differ too greatly or cannot be relaxed, you may have to write your own custom function. I don't have time to do that. I can help show what's done easily, though. For instance, clearing denominators: # * PolynomialLCM @@@ Denominator[#] &@ Take[UpperTriangularize@First@LUDecomposition@ArrayReshape[{{0, 0, c1}, {1, -1, c3}, {-1, 1, c4}}, {3, 3}], 3] Aug 11, 2020 at 15:49
• @Montevideo The desired (first) result seems wrong if $c3+c4=0$ and $c1\ne0$. My result where they are both in the third column comprises $c3+c4\ne0$ and $c1=0$ and the one just mentioned. It's technically not in row-echelon form when $c1=0$, and it's not reduced if both are nonzero. It sort of in the next-to-last step, still comprising all possibilities and waiting for the values of the $c$'s. At least that's how I understand the problem. The only way to get a completely correct result is by cases (a Piecewise result perhaps). Aug 11, 2020 at 23:03