Maxima: reduced row echelon form and pseudo inverse.

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);

addrow(m,i,j,k) := block([n,p],
  require_matrix(m, "first", "addrow"),
  require_integer(i, "second", "addrow"),
  require_integer(j, "third", "addrow"),
  require_rational(k, "fourth", "addrow"),
  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],
  require_matrix(m, "first", "addrow"),
  require_integer(i, "second", "addrow"),
  require_rational(k, "fourth", "addrow"),
  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.