# Find Elementary Matrices that produce RREF

I have a matrix, M, which I am reducing to RREF using RowReduce[M]. I would like to know the elementary matrices which perform the row reducing operations on M during RowReduce[M] to reduce it to RREF, or at least the product of them all.

Update April 29 2023 There is a newer version here Show steps in finding the reduced row echelon form of a symbolic matrix which now works on symbolic matrices also.

edit May 12, 2022. Added frame around pivot element as it moves to make it easier to see.

I've updated displayRREF to now handle all cases (no solution, unique solution, infinite number of solutions).

It also now does RREF only on a matrix on its own if no b vector is given.

But if a b is given as well, then it will also solve the system $$Ax = b$$.

I've kept the original answer below, but that old code can now be replaced by this newer version. One day I might make this a resource function when I have sometime.

Since I wrote this in code cell. Posting the code here might not display as easy to read as in the notebook. So here is a link to the notebook also RREF.nb

First, here are some usage examples.

## Source code

 displayMatrix[A_?(MatrixQ[#]&),dashed_?BooleanQ]:= If[dashed,displayMatrix[A,Last@Dimensions[A]],Print[MatrixForm[A]]]
displayMatrix[A_?(MatrixQ[#]&),nCols_Integer]:= Print[A[[All,1;;nCols]]//matWithDiv[nCols,Background->LightOrange]]
displayRREF[Ain_?(MatrixQ[#]&),dashed_:False]:=Module[
{multiplier,j,i,pivotRow,pivotCol,nRows,nCols,p,tmp,startAgain,A=Ain,n,m,pivotsFound={},keepEliminating,nIter},
displayMatrix[A,dashed];
{nRows,nCols} = Dimensions[A];

keepEliminating=True;

n=1; m=1;
nIter = 0;
While[keepEliminating,
nIter++;
If[nIter>100, (*safe guard*)
Return["Internal error. Something went wrong. Or very large system?",Module]
];

If[dashed && m==nCols,
keepEliminating=False
,
Print["Pivot is A(",n,",",m,")"];
Print@makeNiceMatrix[A,{n,m},dashed];
If[A[[n,m]] != 0,
If[A[[n,m]] != 1,
A[[n,All]] = A[[n,All]]/A[[n,m]];
Print["Making the pivot 1 using using row(",n,")= row(",n,")/A(",n,",",m,")"]
(*Print@makeNiceMatrix[A,{n,m},dashed]*)
];
If[n<nRows,
Do[
If[A[[j,m]] != 0,
multiplier = A[[j,m]]/A[[n,m]];
Print["Zeroing out element A(",j,",",m,") using row(",j,")=",multiplier,"*row(",n,")-row(",j,")"];
A[[j,m;;]] = A[[j,m;;]] - multiplier*A[[n,m;;]];
Print@makeNiceMatrix[A,{n,m},dashed];
]
,{j,n+1,nRows}
];
];
pivotsFound = AppendTo[pivotsFound,{n,m}];
If[n==nRows,
keepEliminating=False
,
n++;
m++
]
,
(*pivot is zero*)
Print["Pivot is zero"];
(*Print@makeNiceMatrix[A,{n,m},dashed];*)
If[n==nRows&&m==nCols, keepEliminating=False
,
(*pivot is zero. If we can find non-zero pivot row below, then exchange rows*)
If[n<nRows,
p = FirstPosition[A[[n+1;;,m]],_?(#!=0&)];
If[ p===Missing["NotFound"]|| Length[p]==0,
If[m<nCols,
m++
,
keepEliminating=False
]
,
(*found non zero pivot below. Exchange rows*)
tmp = A[[n,All]];
A[[n,All]] = A[[First[p]+n,All]];
A[[First[p]+n,All]] = tmp;
Print["Exchanging row(",n,") and row(",First[p]+n,")"]
(*Print@makeNiceMatrix[A,{n,m},dashed]*)
]
,
If[m<nCols,
m++
,
keepEliminating=False
]
]
]
]
]
];

(*pivotsFound = DeleteDuplicates[pivotsFound];*)
Print[">>>>>>Starting backward elimination phase. The pivots are ",pivotsFound];
Do[
pivotRow=First@entry;
pivotCol=Last@entry;

If[pivotRow>1,
Do[
If[ A[[i,pivotCol]] != 0,
Print["Zeroing out element A(",i,",",pivotCol,") using row(",i,")=row(",i,")-A(",i,",",pivotCol,")*row(",pivotRow,")"];
A[[i,;;]] = A[[i,;;]] - A[[i,pivotCol]]*A[[pivotRow,;;]];
Print@makeNiceMatrix[A,{pivotRow,pivotCol},dashed]
]
,
{i,pivotRow-1,1,-1}
]
]
,
{entry,pivotsFound}
];

{A,pivotsFound[[All,2]]}
]
makeSolutionSpecialCase[A_?(MatrixQ[#]&),b_?(VectorQ[#]&),pivotCols_List]:=Module[

Print["Pivot columns are ",MatrixForm[pivotCols]];

ClearAll[x,t]; (*did not make them local, to prevent \$ from showing in print*)
{nRows,nCols} = Dimensions[A];
Print["There are ",nLeadingVariables," leading variables and ",nFreeVariables," free variables. These are "];

Array[t, nFreeVariables];
Array[x, nCols];

m=0;k=0;
Do[
If[Not[MemberQ[pivotCols,n]],
m++;
Print[x[n]," is a free variable. Let ",x[n],"=",t[m]];
AppendTo[variables,t[m]];
AppendTo[sol,0]
,
AppendTo[variables,x[n]];
AppendTo[sol,b[[++k]]]
]
,{n,1,nCols}
];

freeVariables=(t[#]&/@Range[nFreeVariables]);
Print["Hence the system after RREF is the following>>>>>>"];
Print[MatrixForm[A.variables],"=",MatrixForm[b]];
Print["There is different solution for different value of the free variables."];
Print["Setting free variable ", freeVariables," to zero gives"];
variables=variables/.((t[#]->0)&/@Range[nFreeVariables]);
Print[MatrixForm[A.variables],"=",MatrixForm[b]];
Print["Therefore the final solution is "];
Print[MatrixForm[x[#]&/@Range[nCols]],"=",MatrixForm[sol]]

]
(*version 3/10/2017. Original version*)
(*version 5/7/2022. Make it handle all cases*)
(*version 5/12/2022. Added frame around pivot as it moves*)

displayRREF[A_?(MatrixQ[#]&),b_?(VectorQ[#]&)]:=Module[
{B,(*augmented*)nRows,nCols,nEquations,nVariables = Length@b,rref,pivotCols,matRank,augmentedRank},

{nRows, nCols} = Dimensions[A];
nEquations     = nRows;

If[nEquations != nVariables, Return["Size of b vector is not the same as number of rows in the matrix",Module]];

matRank   = MatrixRank[A];
B = Join[A,Transpose[{b}],2];

Print["Augmented matrix is "];
displayMatrix[B,True];

augmentedRank = MatrixRank[B];

If[matRank<augmentedRank, (*Case No solution*)
Print["System is not consistent, no solution exist. Try using least squares."]
,
(*we must have rank A == rank[A|b]  -- system is consistent. It can have*)
(*infinite solutions or one unique solution*)
If[matRank == nCols,
If[nCols == nRows, (*square. case A*)
Print["System is consistent. One unique solution"];
{rref,pivotCols} = displayRREF[B,True];

If[MatrixRank[A]!=Length[pivotCols], (*Verify*)
Print["Internal Error detected. Pivot columns not same as Rank. Please report bug"]
,
Print["Solution  vector is ", MatrixForm[rref[[All, nCols + 1]]]]
]
,
Print["Internal Error detected. matRank != nCols but not square. "]
]
,
(*-- rank A  < N. Case C *)
Print["System is consistent but infinite solutions."];
{rref,pivotCols} = displayRREF[B,True];
makeSolutionSpecialCase[rref[[All,1;;nCols]],rref[[All,nCols+1]],pivotCols]
]
]
]
makeNiceMatrix[mat_?MatrixQ,pivot_List,dashed_?BooleanQ]:=Module[{g,nRow,nCol},
{nRow,nCol} = Dimensions[mat];
(*g=Grid[mat,Dividers->{dashPosition->{Red,Dashed}},Background->{None,None,{pivot->Pink}}];*)
If[dashed,
g = Grid[mat,Dividers->{nCol->{Red,Dashed}},Frame->{None,None,{pivot->True}},Background->LightOrange]
,
g = Grid[mat,Frame->{None,None,{pivot->True}}]
];
MatrixForm[{{g}}]
]
(*makes a dash line inside Matrix*)
Format[matWithDiv[n_,opts:OptionsPattern[Grid]][m_?MatrixQ]]:=MatrixForm[{{Grid[m,opts,Dividers->{n->{Red,Dashed}}]}}];


## OLD ANSWER (do not use any more. Above is newer version)

Updated to Solve Ax=b and display all transformations to reduced Echelon form, and simplified the API. This now also prints the matrix inverse as by product, since it now uses an full augmented matrix.

The idea is to use LUDecomposition first to obtain the permutation rows used for pivoting. Once the order of the rows used is known, then forward elimination is used to generate Echelon form, then backward elimination is used to produce the final reduced echelon form and the solution vector.

Full augmented matrix is used so that the RHS of the augmented matrix will contain the matrix inverse at the end. Verified same inverse is produced as Mathematica Inverse. This only works on matrices that have non-zero determinant.

Here are some usage examples

## Example 1

(*solve Ax=b *)
mat = {{2, 7, 3}, {1, 3, 2}, {3, 7, 9}};
b = {11, 2, -12};
displayRREF[mat, b]


## Example 2

mat={{17,42,-36},{13,45,-34},{12,47,-35}};
b={213,226,197};
displayRREF[mat,b]


## Example 3

mat={{5,2,18,4},{0,1,2,5},{4,1,12,6},{2,3,8,9}};
b={1,2,3,4};
displayRREF[mat,b]


## Function

displayRREF[mat_?(MatrixQ[#] &), b_?(VectorQ[#] &)] :=
Module[{i, j, multiplier, pivot, augmented, m = Length@mat, lu, p, c, tmp, inverse},
(*version 3/10/2017*)

(*check if matrix is singular.Per Daniel Lichtblau post seen Wolfram site*)
(*this is better method than using Det*)
If[MatrixQ[mat] &&
MatrixRank[mat] == Length[mat] == Length[mat[[1]]] === False,
Return["Sorry, but matrix is singular!"]
];

If[Length@b != Length@mat,
Return[
"Size of b vector not the same as number of rows in A matrix"]
];

{lu, p, c} = LUDecomposition[mat];
tmp = lu SparseArray[{i_, j_} /; j >= i -> 1, {Length@mat,
Length@mat}];
(*Print["Mathematica says Echelon form is ", MatrixForm@tmp];*)
tmp = mat[[p, All]];
(*Print["Mathematica says inverse Matrix is", MatrixForm@Inverse@
mat];*)
augmented = Join[mat, Transpose[{b}], 2];
augmented = ArrayFlatten[{{augmented, IdentityMatrix[Length@tmp]}}];

Print[">>>>>>Starting forward Gaussian elimination phase using ",
augmented[[All, 1 ;; m + 1]] //
matWithDiv[m + 1, Background -> LightOrange] ,
MatrixForm[augmented[[All, m + 2 ;;]]]];
Do[
Print["pivot now is (", pivot, ",", pivot, ")" ];
Do[
multiplier = augmented[[j, pivot]]/augmented[[pivot, pivot]];
Print["will now zero out element (", j, ",", pivot,
") by subtracting ", multiplier, " times row ", pivot,
" from row ", j];
augmented[[j, pivot ;;]] =
augmented[[j, pivot ;;]] -
multiplier*augmented[[pivot, pivot ;;]];
Print[
augmented[[All, 1 ;; m + 1]] //
matWithDiv[m + 1, Background -> LightOrange] ,
MatrixForm[augmented[[All, m + 2 ;;]]]]
, {j, pivot + 1, m}
]
, {pivot, 1, m}
];

Print[">>>>>>Starting backward elimination phase"];

Do[
Do[
multiplier = augmented[[j, pivot]]/augmented[[pivot, pivot]];
Print["will now zero out element (", j, ",", pivot,
") by subtracting ", multiplier, " times row ", pivot,
" from row ", j];
augmented[[j, pivot ;;]] =
augmented[[j, pivot ;;]] -
multiplier*augmented[[pivot, pivot ;;]];
Print[
augmented[[All, 1 ;; m + 1]] //
matWithDiv[m + 1, Background -> LightOrange] ,
MatrixForm[augmented[[All, m + 2 ;;]]]]
, {j, 1, pivot - 1}
]
, {pivot, 2, m}
];

Print[">>>>>>Starting Final phase, convert reduced echelon to     identity matrix"];
Do[
augmented[[j, ;;]] = augmented[[j, ;;]]/augmented[[j, j]]
, {j, 1, m}
];

Print[augmented[[All, 1 ;; m + 1]] //
matWithDiv[m + 1, Background -> LightOrange] ,
MatrixForm[augmented[[All, m + 2 ;;]]]];
(*flip at inverse Matrix now using p from LUDecomposition,
but using column wise*)
Print["Inverse Matrix is ",
MatrixForm[ augmented[[All, m + 2 ;;]] ]];
Print["Solution  vector is ", MatrixForm[augmented[[All, m + 1]]]]
];

(*makes a dash line inside Matrix*)
Format[matWithDiv[n_, opts : OptionsPattern[Grid]][m_?MatrixQ]] :=
MatrixForm[{{Grid[m, opts, Dividers -> {n -> {Red, Dashed}}]}}];

• Is there any way to get to this to work with a singular matrix as RREF certainly works for the augmented system? mat = {{1, -2, 4, 3, 2}, {-7, 14, -28, -12, -23}, {4, -8, 12, 12, 8}, {-3, 6, -4, -1, -14}}; b = {2, 4, 18, -10}; displayRREF[mat, b]
– Moo
May 7, 2022 at 0:08
• @Moo Yes, it is possible. But just a matter to getting the time to do it as I am busy with other things. I hope to be able to do that sometime in the future. May 7, 2022 at 0:10
• @Moo I've updated it to handle all cases. Any bugs please let me know. May 9, 2022 at 5:13
• Thanks! I tried this: mat = {{1, -2, 4, 3, 2}, {-7, 14, -28, -12, -23}, {4, -8, 12, 12, 8}, {-3, 6, -4, -1, -14}}; b = {2, 4, 18, -10}; displayRREF[mat, b]. it worked perfectly Thank you for updating it! (+1). By the way, i think it is great idea to make this a Resource Function!
– Moo
May 9, 2022 at 10:43
• @Nasser Hi! I want to limit the matrix to RREF by adding a multiple of one row to another row instead of column transformation and multiple of the row itself. If RREF cannot be achieved by this method, "Impossible" will be returned. How should I modify your code? thank you! Nov 3, 2022 at 3:30

The following answer is a contribution to enrich Nasser's code a little more:

NiceMatrixExtended[mat_?MatrixQ, pivot_?VectorQ] /;
Length[pivot] == 2 :=
Module[{g, nRow, nCol, highlightedMat}, {nRow, nCol} = Dimensions[mat];
(*A version of the matrix with the pivot highlighted using Item*)
highlightedMat = mat;
highlightedMat[[pivot[[1]], pivot[[2]]]] =
Item[mat[[pivot[[1]], pivot[[2]]]], Background -> LightBlue];
g = Grid[highlightedMat, Background -> LightYellow,
Dividers -> {{nCol - nRow -> Directive[Red, Dashed, Thickness[1.5]],
nCol - nRow + 1 -> Directive[Red, Dashed, Thickness[1.5]]}}];
MatrixForm[{{g}}]]


Testing NiceMatrixExtended:

system = {5 w + 2 x + 3 y - 6 z == 31, w + x - y - 4 z == 2,
2 w - x + 7 y - 5 z == 16, 3 w + 3 x + 2 y - 5 z == 23};
vars = {x, y, z, w};

(*Extended Augmented System Matrix*)
asys = {{2, 3, -6, 5, 31, 1, 0, 0, 0}, {1, -1, -4, 1, 2, 0, 1, 0, 0},
{-1, 7, -5, 2, 16, 0, 0, 1, 0}, {3, 2, -5, 3, 23, 0, 0, 0, 1}};
NiceMatrixExtended[asys, {1, 1}]


https://i.sstatic.net/hv0gQ.png