# [Maxima] Maxima: reduced row echelon form and pseudo inverse.

Wolfgang Lindner LindnerW at t-online.de
Sat Sep 8 01:13:00 CDT 2007

```dear C. Frangos,

here is a version of rref by Antoine Chambert-Loir/France.
antoine.chambert-loir at univ-rennes1.fr
http://perso.univ-rennes1.fr/antoine.chambert-loir/

It was mentioned on this list some weeks ago.

>--------------------------------------
> (1) reduced row echelon form.

/*
Description : forme echelonne par lignes d'une matrice rectangulaire
(a  coefficients dans un corps commutatif).
*/

request_rational_matrix(m, pos, fn) :=
if every('identity, map(lambda([s], every('ratnump,s)), args(m))) then true else
print("Some entries in the matrix are not rational numbers. The result might be wrong.");

rowswap(m,i,j) := block([n, p, r],
require_matrix(m, "first", "rowswap"),
require_integer(i, "second", "rowswap"),
require_integer(j, "third", "rowswap"),
n : length(m),
if (i < 1) or (i > n) or (j < 1) or (j > n)
then error("Array index out of bounds"),
p : copymatrix(m),
r : p[i],
p[i] : p[j],
p[j] : r,
p);

n : length(m),
if (i < 1) or (i > n) or (j < 1) or (j > n)
then error("Array index out of bounds"),
p : copymatrix(m),
p [i] : p[i] + k * p[j],
p);

rowmul(m,i,k) := block([n,p],
n : length(m),
if (i < 1) or (i > n) then error("Array index out of bounds"),
p : copymatrix(m),
p [i] : k * p[i],
p);

rref(m):= block([p,nr,nc,i,j,k,pivot,pivot_row,debug],
debug : 0,
request_rational_matrix(m," ","rref"),
nc: length(first(m)),
nr: length(m),
if nc = 0 or nr = 0 then
error ("The argument to 'rref' must be a matrix with one or more rows and columns"),

p:copymatrix(m),
ci : 1, cj : 1,
while (ci<=nr) and (cj<=nc) do
(
if (debug = 1) then (
disp(p),
print("curseur en ligne ",ci," et colonne ",cj)),
pivot_row : 0, pivot : 0,
for k : ci thru nr do (
if ( abs(p[k,cj]) > pivot ) then (
pivot_row : k,
pivot : abs(p[k,cj]))),
if (debug = 1) then
print("colonne ",cj," : pivot trouve ligne ", pivot_row,", valeur : ",pivot),
if (pivot = 0) then (cj : cj +1)
else (
p : rowswap(p,ci,pivot_row),
if (debug = 1) then      print (".. Echange : ",p),
p : rowmul(p,ci,1/p[ci,cj]),
if (debug = 1) then     print (".. Normalisation : ",p),
for k : 1 thru nr do (
if not (k=ci) then (p : addrow (p,k,ci,-p[k,cj]))),
ci : ci+1, cj : cj+1)),
p
);

>--------------------------------------
> (2) the Moore-Penrose pseudo inverse.

As for the MPI (=general inverse Ginv) I don't know of a MAXIMA implementation.
But John Fox gives a function on the R list, see the thread
https://stat.ethz.ch/pipermail/r-help/2007-September/139923.html :

RREF <- function(X, ...) GaussianElimination(X, ...)
# returns the reduced row-echelon form of X

Ginv <- function(A, tol=sqrt(.Machine\$double.eps), verbose=FALSE,
fractions=FALSE){
# return an arbitrary generalized inverse of the matrix A
# A: a matrix
# tol: tolerance for checking for 0 pivot
# verbose: if TRUE, print intermediate steps
# fractions: try to express nonintegers as rational numbers
m <- nrow(A)
n <- ncol(A)
B <- GaussianElimination(A, diag(m), tol=tol, verbose=verbose,
fractions=fractions)
L <- B[,-(1:n)]
AR <- B[,1:n]
C <- GaussianElimination(t(AR), diag(n), tol=tol, verbose=verbose,
fractions=fractions)
R <- t(C[,-(1:m)])
AC <- t(C[,1:m])
ginv <- R %*% t(AC) %*% L
if (fractions) fractions (ginv) else round(ginv, round(abs(log(tol,
10))))
}

Maybe this info is of some help.

HTH Wolfgang

"C. Frangos" <cfrangos at telkomsa.net> schrieb:
> I am interested in obtaining Maxima programs that can compute for matrices
> with symbolic variables the following:
>
> (1) reduced row echelon form.
> (2) the Moore-Penrose pseudo inverse.
>
> Any assistance would be much appreciated.
>
> Regards,
>
> C. Frangos.

```