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 18put('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 25eval_when([translate, compile], translate_fast_arrays : false); 26 27eval_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 33eval_when([batch,load,translate,compile], 34 put('infolevel, 'silent, 'linearalgebra), 35 load("load-linearalgebra-lisp-files")); 36 37require_integer(i, pos, fn) := 38 if integerp(i) then true else 39 error("The", pos, "argument to" ,fn, "must be an integer"); 40 41require_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 45request_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 49dotproduct(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 56nullspace(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 78nullity(m) := length(nullspace(m)); 79 80orthogonal_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 89locate_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 117columnspace(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. 138I don't approve of rank(matrix([a,b],[c,d])) --> 2. This version will 139print the proviso warnings. 140*/ 141 142linalg_rank(m) := length(columnspace(m)); 143 144rowswap(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 157columnswap(m,i,j) := transpose(rowswap(transpose(m),i,j)); 158 159/* row(m,i) <-- row(m,i) - theta * row(m,i) */ 160 161rowop(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 168columnop(m,i,j, theta) := transpose(rowop(transpose(m),i, j, theta)); 169 170hipow_gzero(e,x) := block([n : hipow(e,x)], if n > 0 then n else 'inf); 171good_pivot(e,x) := freeof(x,e) and e # 0; 172 173ptriangularize(m,v) := block([p : ptriangularize_with_proviso(m,v)], 174 if not(emptyp(p[2])) then print("Proviso: ",p[2]), 175 p[1]); 176 177ptriangularize_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 189column_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 225mat_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 233mat_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 239mat_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 245mat_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 254kronecker_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 259diag_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 273hilbert_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 282hankel([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 293toeplitz([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 304polytocompanion(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 312moore_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)); 316