Maxima: reduced row echelon form and pseudo inverse.

Hi Wolfgang,

Thanks for the fast response and for attaching the program below - much 
appreciated.

I did try this program, but this seems to work only for matrices with 
numerical entries (rationals, etc), but not symbolic variables. Its not clear 
to me whether or how this program can be modified to handle symbolic matrices.

Regards,

C. Frangos.


On Saturday 08 September 2007 08:13, you wrote:
> 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.