share/linearalgebra/linearalgebra.mac

   1 /*
   2    Author Barton Willis
   3    University of Nebraska at Kearney
   4    Copyright (C) 2005, Barton Willis
   5
   6    Brief Description: Maxima code for finding the nullspace and column space
   7    of a matrix, along with code that triangularizes matrices using the method
   8    described in ``Eigenvalues by row operations,'' by Barton Willis.
   9
  10    This is free software  you can redistribute it and/or
  11    modify it under the terms of the GNU General Public License,
  12    http://www.gnu.org/copyleft/gpl.html.
  13
  14    This software has NO WARRANTY, not even the implied warranty of
  15    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
  16 */
  17
  18 put('linearalgebra,1,'version);
  19
  20 /* The translator generates code that doesn't work when
  21    translate_fast_arrays is true. So set translate_fast_arrays
  22    to false.
  23 */
  24
  25 eval_when([translate, compile], translate_fast_arrays : false);
  26
  27 eval_when(translate,
  28           declare_translated(ptriangularize_with_proviso,locate_matrix_entry,hankel,toeplitz,
  29                              nullspace,require_integer,columnspace,rowswap,rowop,require_symbol,
  30                              mat_fullunblocker,mat_trace,mat_unblocker, column_reduce,good_pivot,
  31                              hipow_gzero, mat_norm))$
  32
  33 eval_when([batch,load,translate,compile],
  34   put('infolevel, 'silent, 'linearalgebra),
  35   load("load-linearalgebra-lisp-files"));
  36
  37 require_integer(i, pos, fn) :=
  38   if integerp(i) then true else
  39      error("The", pos, "argument to" ,fn, "must be an integer");
  40
  41 require_symbol(x, pos, fn) :=
  42   if symbolp(x) or subvarp(x) then true else
  43      error("The", pos, "argument to", fn, "must be a symbol");
  44
  45 request_rational_matrix(m, pos, fn) :=
  46   if every('identity, map(lambda([s], every('ratnump,s)), args(m))) then true else
  47     print("Some entries in the matrix are not rational numbers. The result might be wrong.");
  48
  49 dotproduct(a,b) := block([scalarmatrixp : true],
  50   require_matrix(a,"first", "dotproduct"),
  51   require_matrix(b,"second", "dotproduct"),
  52   if matrix_size(a) # matrix_size(b) or second(matrix_size(a)) # 1 then
  53      error("Arguments to dotproduct must be vectors with the same size"),
  54   ctranspose(a) . b);
  55
  56 nullspace(m) := block([nr, nc, acc : set(), proviso : true, pv, prederror : false],
  57
  58   require_unblockedmatrix(m, "first", "nullspace"),
  59   /* nc and nr are the sizes of the transpose of m */
  60
  61    [nc, nr] : matrix_size(m),
  62    m : triangularize(addcol(transpose(m), ident(nr))),
  63    for row : 1 thru nr do (
  64      pv : locate_matrix_entry(m, row, 1, row, nc, lambda([s], compare(s,0) # "="), 'bool),
  65
  66      if not(listp(pv)) then (
  67        acc : adjoin(transpose(genmatrix(lambda([ii,j], m[row,j+nc]),1,nr)),acc))
  68      else (
  69        proviso : proviso and notequal(m[row, second(pv)],0))),
  70
  71    if proviso # true then (
  72      print("Proviso: ",proviso),
  73      put(concat(outchar, linenum), proviso, 'proviso)),
  74
  75    subst('span, 'set, subset(acc, lambda([s], some(lambda([x], compare(x,0) # "="), s)))));
  76
  77
  78 nullity(m) := length(nullspace(m));
  79
  80 orthogonal_complement([v]) :=  block([sz],
  81   require_unblockedmatrix(first(v)," ","orthogonal_complement"),
  82   sz : matrix_size(first(v)),
  83   if second(sz) # 1 then error("Each argument must be a column vector"),
  84   for vi in v do (
  85     require_matrix(vi," ","orthogonal_complement"),
  86     if matrix_size(vi) # sz then error("Each vector must have the same dimension")),
  87   nullspace(transpose(rreduce('addcol, v))));
  88
  89 locate_matrix_entry(m, r1, c1, r2, c2, fn, rel) := block([im, cm, mf, e,
  90      nr, nc, ok : true, frel],
  91   require_unblockedmatrix(m, "first", "locate_matrix_entry"),
  92   require_integer(r1, "second", "locate_matrix_entry"),
  93   require_integer(r2, "second", "locate_matrix_entry"),
  94   require_integer(c1, "third", "locate_matrix_entry"),
  95   require_integer(r2, "fourth", "locate_matrix_entry"),
  96   nr : length(m),
  97   nc : length(first(m)),
  98
  99   if (c1 > c2) or (r1 > r2) or (r1 < 1) or (c1 < 1) or (r2 < 1) or (c2 < 1) or
 100      (c1 > nc) or (c2 > nc) or (r1 > nr) or (r2 > nr) then error("Bogus submatrix"),
 101
 102   mf : if rel = 'min then 'inf else if rel = 'max then 'minf else 0,
 103   frel : if rel = 'max then ">" else if  rel = 'min then "<" else if rel = 'bool
 104     then lambda([x,y],x) else error("Last argument must be 'max' or 'min'"),
 105   im : 0,
 106   cm : 0,
 107   for i : r1 thru r2 while ok do (
 108     for j : c1 thru c2 while ok do (
 109        e : apply(fn, [m[i,j]]),
 110        if is(apply(frel, [e, mf])) then (
 111          im : i,
 112          cm : j,
 113          if rel = 'bool then ok : false,
 114          mf : e))),
 115   if im > 0 and cm > 0 then [im, cm] else false);
 116
 117 columnspace(a) := block([m, nc, nr, cs : set(), pos, x, proviso : true, algebraic : true, c_min],
 118   require_unblockedmatrix(a, "first","columnspace"),
 119   [nr, nc] : matrix_size(a),
 120   if nc = 0 or nr = 0 then (
 121     error("The argument to 'columnspace' must be a matrix with one or more rows and columns")),
 122
 123   m : triangularize(a),
 124   c_min : 1,
 125   for i : 1 thru nr while c_min <= nc do (
 126    if listp(pos : locate_matrix_entry(m, i, c_min, i, nc, lambda([x], x # 0), 'bool)) then (
 127      c_min : second(pos) + 1,
 128      if constantp(x : part(part(m, first(pos)),second(pos))) = false then (
 129        proviso : proviso and notequal(ratsimp(x), 0)),
 130      cs : adjoin(col(a,second(pos)) ,cs))),
 131
 132   if proviso # true then (
 133     print("Proviso: ",proviso),
 134     put(concat(outchar, linenum), proviso, 'proviso)),
 135   funmake('span, listify(cs)));
 136
 137 /* Maxima's rank function doesn't complain with symbolic matrices.
 138 I don't approve of rank(matrix([a,b],[c,d])) --> 2. This version will
 139 print the proviso warnings.
 140 */
 141
 142 linalg_rank(m) := length(columnspace(m));
 143
 144 rowswap(m,i,j) := block([n, p, r],
 145   local(p),
 146   require_matrix(m, "first", "rowswap"),
 147   require_integer(i, "second", "rowswap"),
 148   require_integer(j, "third", "rowswap"),
 149   n : length(m),
 150   if (i < 1) or (i > n) or (j > n) then error("Array index out of bounds"),
 151   p : copymatrix(m),
 152   r : p[i],
 153   p[i] : p[j],
 154   p[j] : r,
 155   p);
 156
 157 columnswap(m,i,j) := transpose(rowswap(transpose(m),i,j));
 158
 159 /*  row(m,i) <-- row(m,i) - theta * row(m,i) */
 160
 161 rowop(m,i,j,theta) := block([p : copymatrix(m), listarith : true],
 162   local(p),
 163   p[i] : p[i] - theta * p[j],
 164   p);
 165
 166 /*  col(m,i) <-- col(m,i) - theta * col(m,i) */
 167
 168 columnop(m,i,j, theta) := transpose(rowop(transpose(m),i, j, theta));
 169
 170 hipow_gzero(e,x) := block([n : hipow(e,x)], if n > 0 then n else 'inf);
 171 good_pivot(e,x) := freeof(x,e) and e # 0;
 172
 173 ptriangularize(m,v) := block([p : ptriangularize_with_proviso(m,v)],
 174    if not(emptyp(p[2])) then print("Proviso: ",p[2]),
 175    p[1]);
 176
 177 ptriangularize_with_proviso(m,v) := block([nr, nc, proviso : [], mp],
 178   require_unblockedmatrix(m, "first", "ptriangularize"),
 179   require_symbol(v, "second", "ptriangularize"),
 180   nr : length(m),
 181   nc : length(first(m)),
 182   for i : 1 thru min(nr, nc) do (
 183      mp : column_reduce(m,i,v),
 184      m : first(mp),
 185      proviso : append(proviso, second(mp))),
 186   proviso : setify(proviso),
 187   [expand(m), ratsimp(proviso)]);
 188
 189 column_reduce(m,i,x) := block([nc, nr, pos, q, proviso : []],
 190    /* require_matrix(m, "first", "column_reduce"), */
 191    nc : length(first(m)),
 192    nr : length(m),
 193    /* require_integer(i, "second", "column_reduce"),
 194    if (i < 1) or (i > (length(first(m)))) then error("Bogus matrix column"), */
 195
 196    while not(good_pivot(m[i,i],x)) and (i+1 <= nr) and
 197         listp(pos : locate_matrix_entry(m,i+1,i,nr,i, lambda([e],e # 0), 'bool)) do (
 198      pos : locate_matrix_entry(m,i,i,nr,i,lambda([e], good_pivot(e,x)), 'bool),
 199      if listp(pos) then (
 200        m : rowswap(m,i,pos[1]))
 201      else (
 202        m : expand(m),
 203        pos : locate_matrix_entry(m, i, i, nr, i, lambda([e], hipow(e,x)), 'max),
 204        m : rowswap(m,i,pos[1]),
 205        pos : locate_matrix_entry(m, i+1, i, nr, i, lambda([e], hipow_gzero(e,x)), 'min),
 206        if listp(pos) then (
 207          q : divide(m[i,i], m[pos[1],i],x),
 208          m : rowop(m,i,pos[1],first(q))))),
 209
 210      pos : locate_matrix_entry(m, i, i, nr, i, lambda([e], good_pivot(e,x)), 'bool),
 211
 212      if listp(pos) then (
 213         m : rowswap(m,i,first(pos)),
 214         for j : i+1 thru nr do (
 215           if not(constantp(ratnumer(m[i,i]))) then (
 216             proviso : cons(ratnumer(m[i,i]) # 0, proviso)),
 217            if not(constantp(ratdenom(m[i,i]))) then (
 218             proviso : cons(ratdenom(m[i,i]) # 0, proviso)),
 219           m : rowop(m,j,i,m[j,i]/m[i,i])))
 220      else (
 221        pos : locate_matrix_entry(m, i, i, nr, i, lambda([e], e # 0), 'bool),
 222        if listp(pos) then m : rowswap(m,i, first(pos))),
 223     [m, proviso]);
 224
 225 mat_norm(m, p) := block([listarith : true],
 226   require_matrix(m,"first","mat_norm"),
 227   if blockmatrixp(m) then m : mat_fullunblocker(m),
 228   if p = 1 then lmax(xreduce("+", map(lambda([s], map('cabs, s)), if args(m) = [] then [[]] else args(m))))
 229   else if p = 'inf then lmax(map(lambda([s], xreduce("+", map('cabs, s))), args(m)))
 230   else if p = 'frobenius then sqrt(xreduce("+", lreduce('append, args(cabs(m)^2))))
 231   else error("Not able to compute the ",p," norm"));
 232
 233 mat_fullunblocker(m) := block(
 234   require_matrix(m, "first", "mat_unblocker"),
 235   if blockmatrixp(m) then (
 236      mat_fullunblocker(lreduce('addrow, args(map(lambda([x], lreduce('addcol, x)), m)))))
 237    else m);
 238
 239 mat_unblocker(m):=block(
 240   require_matrix(m, "first", "mat_unblocker"),
 241   if blockmatrixp(m) then (
 242      lreduce('addrow, args(map(lambda([x], lreduce('addcol, x)), m))))
 243    else m);
 244
 245 mat_trace(m) := block([n, acc : 0],
 246   if not matrixp(m) then funmake('mat_trace, [m]) else (
 247     n : matrix_size(m),
 248     if first(n) # second(n) then error("The first argument to 'mat_trace' must be a square matrix"),
 249     n : first(n),
 250     for i from 1 thru n do (
 251       acc : acc + if matrixp(m[i,i]) then mat_trace(m[i,i]) else m[i,i]),
 252     acc));
 253
 254 kronecker_product(a,b) := (
 255   require_matrix(a,"first", "kronecker_product"),
 256   require_matrix(b,"second", "kronecker_product"),
 257   mat_unblocker(outermap('matrix_element_mult, a,b)));
 258
 259 diag_matrix([d]) := block([f, n, prederror : false],
 260
 261   if every(lambda([s], matrixp(s) and not(blockmatrixp(s))), d) then (
 262     f : lambda([i,j], if i = j then inpart(d,i)
 263                else zeromatrix(first(matrix_size(inpart(d,i))), second(matrix_size(inpart(d,j))))))
 264
 265   else if some(lambda([s], matrixp(s) or blockmatrixp(s)), d) then
 266      error("All arguments to 'diag_matrix' must either be unblocked matrices or non-matrices")
 267
 268   else f : lambda([i,j], if i = j then inpart(d,i) else 0),
 269
 270   n : length(d),
 271   genmatrix(f, n, n));
 272
 273 hilbert_matrix(n) := block([row, mat : []],
 274   mode_declare(n,integer),
 275   require_posinteger(n,"first","hilbert_matrix"),
 276   row : makelist(1/i,i,1,n),
 277   for i : 1 thru n do (
 278     mat : cons(row, mat),
 279     row : endcons(1/(n + i), rest(row))),
 280   funmake('matrix, reverse(mat)));
 281
 282 hankel([q]) := block([col,row,m,n,partswitch : false],
 283   if length(q) > 2 or length(q) < 1 then error("The function 'hankel' requires one or two arguments"),
 284   col : inpart(q,1),
 285   row : if length(q) = 2 then inpart(q,2) else map(lambda([x],0), col),
 286
 287   require_list(row,"first","hankel"),
 288   require_list(col,"second","hankel"),
 289   m : length(row),
 290   n : length(col),
 291   genmatrix(lambda([i,j],if i+j-1 <= n then inpart(col,i+j-1) else inpart(row,i+j-n)),n,m));
 292
 293 toeplitz([q]) := block([col,row,m,n,partswitch : false],
 294   if length(q) > 2 or length(q) < 1 then error("The function 'toeplitz' requires one or two arguments"),
 295   col : inpart(q,1),
 296   row : if length(q) = 2 then inpart(q,2) else map('conjugate, col),
 297
 298   require_list(row,"first","toeplitz"),
 299   require_list(col,"second","toeplitz"),
 300   m : length(row),
 301   n : length(col),
 302   genmatrix(lambda([i,j], if i -j >= 0 then inpart(col, i-j+1) else inpart(row,j-i+1)),n,m));
 303
 304 polytocompanion(p,x) := block([n],
 305   if not polynomialp(p,[x], lambda([e], freeof(x,e))) then
 306      error("First argument to 'polytocompanion' must be a polynomial"),
 307   p : expand(p),
 308   n : hipow(p,x),
 309   p : multthru(p / coeff(p,x,n)),
 310   genmatrix(lambda([i,j], if i-j=1 then 1 else if j = n then -coeff(p,x,i-1) else 0),n));
 311
 312 moore_penrose_pseudoinverse(m) := block([mm, z : gensym(), scalarmatrixp : false],
 313   require_unblockedmatrix(m, "first", "moore_penrose_pseudoinverse"),
 314   mm : m . ctranspose(m),
 315   limit(mm : ctranspose(m) . (mm - z * identfor(mm))^^-1,z,0));