matrices.mac

/* vim: ft=maxima
*/

/* does the matrix have an all-zero row? Return the index of the highest
 * one if it exists, and -1 if there's no zero row */

hasZeroRow(a) := block([width, height, i, output],
    height : matrix_size(a)[1],
    width : matrix_size(a)[2],
    output : -1,
    for i thru height do (
        if (apply('matrix, [a[i]]) = zeromatrix(1, width)) 
            then (output : i, return())),
    output
);

/* returns true if there is a zero row with a nonzero row below it,
 * otherwise false */

zeroRowsNotAllAtBottom(a) := block([zr, sm, i],
    zr : hasZeroRow(a),
    if (zr = -1) then return(false),
    sm : apply(submatrix, append(makelist(i, i, 1, zr - 1), [a])),
    return(not is(sm = zeromatrix(matrix_size(sm)[1], matrix_size(sm)[2])))
);


/* return the coordinates of a non-1 leading entry, or [-1, -1]
 * if none exists */

nonOneLeadingEntry(a) := block([height, i, c, output],
    height : matrix_size(a)[1],
    output : [-1, -1],
    for i thru height do (
        c : columnOfLeadingEntry(a, i),
        if (not is(c = -1) and not is(a[i][c] = 1)) then (output : [i, c], return())),
    output
);

/* returns the coordinates of a leading entry such that there is a
 * nonzero entry somewhere else in its column if such a leading entry
 * exists, otherwise returns [-1, -1] */

nonZeroEntryInColumnOfLeadingEntry(a) := block([height,output, i, ii, c],
    height : matrix_size(a)[1],
    output : [-1, -1],
    for i thru height do (
        c : columnOfLeadingEntry(a, i),
        if not is(c = -1) then (
            for ii thru height do (
                if (not is(ii = i) and not is(a[ii][c] = 0)) then (output : [i, c], return())))),
    output
);

/* return [i, j] such that i < j and the leading entry in row j is 
 * not to the right of the leading entry in row i, if such a pair
 * exists, otherwise return [-1, -1] */

badLeadingEntryPositions(a) := block([height, i, c, ii, output, cc],
    height : matrix_size(a)[1],
    output : [-1, -1],
    for i thru (height - 1) do (
        c : columnOfLeadingEntry(a, i),
        if not is(c = -1) then (
            for ii : (i + 1) thru height do (
                cc : columnOfLeadingEntry(a, ii),
                if ((not is(cc = -1)) and (cc <= c)) then
                    (output : [i, ii], return())))),
    output
    );

/* return the first column of a with a nonzero entry in row i if this exists,
 * otherwise -1 */
columnOfLeadingEntry(a, i) := block([q, l, j],
    q : matrix_size(a)[2],
    l : -1,
    for j thru q do (
        if not(a[i, j] = 0) then (l : j, return())
        ),
    l
);

/* compute the row reduced echelon form of a */
rref(a) := block([p, l, i, j],
    p : matrix_size(a)[1],
    a : echelon(a),
    for i : 2 thru p do
        (l : columnOfLeadingEntry(a,i),
        if l > 0 then (
            for j thru i - 1 do
                a : rowop(a, j, i, a[j, l])
            )
        ),
    a
);

 
/* return true if the matrix a is in row reduced echelon form, otherwise
 * false */
isrref(a) := is(a = rref(a));

/* return a string consisting of an HTML<p> containing a <ul> describing why the input matrix is not in RREF */
why_not_rref(mx) := block([mzr, whynot, ble, badLEcol, blep],
    mzr : zeroRowsNotAllAtBottom(mx),

    whynot : "<p><ul>",
    if mzr then whynot : sconcat(whynot, "<li>Your matrix has a row of zeroes which is not at the bottom. Any zero rows of a RREF matrix must be at the bottom. </li>"),

    ble : nonOneLeadingEntry(mx),
    if not (is(ble = [-1, -1])) then whynot : sconcat(whynot, "<li>Your
matrix has a leading entry at position ", ble[1], ", ", ble[2], " which isn't equal to 1. In a RREF matrix all leading entries must be 1. </li>"),

    badLEcol : nonZeroEntryInColumnOfLeadingEntry(mx),
    if not is(badLEcol = [-1, -1]) then whynot : sconcat(whynot, "<li>Your matrix has a leading entry in position ", badLEcol[1], ", ", badLEcol[2], " and there is another nonzero entry in the same column. In a RREF matrix, if a column contains a leading entry, every other entry in that column must be zero. </li>"),

    blep : badLeadingEntryPositions(mx),

    if not is(blep = [-1, -1]) then whynot : sconcat(whynot, "<li>Your matrix has leading entries in row ", blep[1], " and row ", blep[2], ".  The leading entry in row ", blep[2], " is to the left of the leading entry in row ", blep[1], ".  In a RREF matrix, if rows i and j have leading entries and i < j then the leading entry in row j must be to the right of the one in row i.</li>"),

    whynot : sconcat(whynot, "</ul></p>"),
    whynot
);


/* return true if the matrix a is diagonal, otherwise false */
isdiag(a) := block([li, i],
    li : makelist(a[i][i], i, 1, min(matrix_size(a)[1], matrix_size(a)[2])),
    is(a = apply(diag_matrix, li))
);

/* return true if all entries of a are zero */
iszeromx(a) := is(a = zeromatrix(matrix_size(a)[1], matrix_size(a)[2]));

/* return true if v is an eigenvector of a */
isevec(v,a) :=  not(iszeromx(v)) and is(rank(addcol(v, a.v)) = 1); 


/* row ops and elementary matrices */

/* lij(a,l,i,j) adds l times row i to row j of a, modifying a in place */
lij(a, l, i, j) := block(a[j] : l* a[i] + a[j]);
/* corresponding nxn elementary matrix */
e(i, j, l, n) := block([a], a : ident(n),  lij(a, l, i, j), a);

/* sw(a,i,j) swaps rows i and j of a, modifying a in place */
sw(a, i, j) := block(temp : a[i], a[i] : a[j], a[j] : temp);
/* corresponding nxn elementary matrix */
p(i, j, n) := block([a], a : ident(n), sw(a, i, j), a);

/* mu(a,i,l) multiplies row i of a by l, modifying a in place */
mu(a, i, l) := block(a[i] : l * a[i]);
/* corresponding nxn elementary matrix */
delta(i, l, n) := block([a], a : ident(n), mu(a, i, l), a);


/* linear algebra */

/* isLI(v1, v2, ...) is true if the vectors vi are linearly independent,
 * otherwise false */
isLI([vecs]) := block([n,m],
    n : length(vecs),
    m : apply(addcol, vecs),
    is(n = rank(m))
);

/* we represent a subspace by matrix whose columns are a spanning set
 * for the subspace. */

/* dim is an alias for rank */
dim(a) := rank(a);

/* create all of R^n */
fullColumnSpaceOfDimension(n) := ident(n);

/* create the zero subspace of R^n */
zeroSubspace(n) := zeromatrix(n, 1);

/* are the subspaces rep by a and b equal? */
equal_subspaces(a,b) := block([ra, rb, rab],
    if not is(matrix_size(a)[1] = matrix_size(b)[1])
        then return(false),
    ra : rank(a),
    rb : rank(b),
    rab : rank(mat_unblocker(matrix([a,b]))),
    is(ra = rb and ra = rab)
);

/* return a matrix whose columns are a basis for the subspace
 * represented by a, unless a is the zero matrix in which case return a
 * column of zeroes.  */
refine(a) := block([s],
    if (a = zeromatrix(matrix_size(a)[1], matrix_size(a)[2]))
        then return (zeromatrix(matrix_size(a)[1], 1)), 
    s : args(columnspace(a)),
    apply(addcol, s)
);

/* compute transpose(a).a and wrap the result as a matrix if necessary - 
 * maxima converts 1x1 matrices to numbers, we don't want this */
safeAtA(a) := block(
    if (matrix_size(a)[2] = 1)
        then return(matrix([transpose(a).a])),
    transpose(a). a
);


/* orthogonal projection onto column space of a matrix a whose
 * columns are assumed linearly independent */
proj(a) := block([aa],
    aa : refine(a),
    aa . invert(safeAtA(aa)) . transpose(aa)
);


/* is the column space as big as possible? */
spansWholeSpace(a) := is(rank(a) = matrix_size(a)[1]);

/* take a matrix, return matrix representation of its kernel. Don't call
 * this on matrices with zero kernel */
kerMX(a) := apply(addcol, args(nullspace(a)));

/* take the matrix representations a and b of two subspaces A and B,
 * return the matrix rep of A \cap B */
sum_subspacesMX(a,b) := mat_unblocker(matrix([a,b]));

/* compute a matrix representation for the intersection of the subspaces
 * represented by the matrices a and b. The idea is to use the fact that
 * the kernel of
 * [ X ]
 * [ Y ]
 * is the intersection of the kernels of X and Y, so we only have to
 * produce a matrix whose kernel is the column space (image) of a.
 * If A has columns spanning
 *   ker (a^T)= im(a)^\perp
 * then 
 *   ker A^T = im a
 * (certainly im a is contained in ker A^T, then check dimensions). 
 * So we return the kernel of
 * [ A^T ]
 * [ B^T ]
 */
intersect_subspacesMX(a,b) := block([x, y, z],
    if spansWholeSpace(a)
        then return (b)
    else if spansWholeSpace(b)
        then return (a),
    x : kerMX(transpose(a)),
    y : kerMX(transpose(b)),
    z : mat_unblocker(matrix([transpose(x)], [transpose(y)])),
    kerMX(z)
);

/* is this vector v in this subspace a? */
inn(v, a) := is(rank(a) = rank(addcol(v,a)));

/* random permutation matrix */
rpm(n) := block([z, l, p, i], 
    z : zeromatrix(n,n),
    l : makelist(i,i,1,n),
    p : rand_selection(l, n),
    for i : 1 thru n do 
        z[i][p[i]] : 1,
    z
);

/* random lower unitriangular matrix with entries in -maxi...maxi */
random_lut_mx(n, maxi) := block([id, i, j],
    id : ident(n),
    for i : 2 thru n do (
        for j : 1 thru i-1 do (
            id[i][j] : rand_with_step(-maxi, maxi, 1)
            )
        ),
    id
);

/* random element of SL_n(Z) */
random_slnz(n, maxi) := random_lut_mx(n, maxi) . rpm(n) . random_lut_mx(n, maxi);

/* random mx of given size with integer entries from -a..a */
rand_mx(nrows, ncols, a) := block([entries, z, i, j],
    entries : makelist(i, i, -a, a),
    z : zeromatrix(nrows, ncols),
    for i : 1 thru nrows do (
        for j : 1 thru ncols do (
            z[i][j] : rand(entries)
            )
        ),
    z
);

/* random mx of given size with all entries distinct */
rand_mx_distinct(nrows, ncols) := block([count, entries, z, i, j],
    entries : makelist(i, i, -nrows*ncols, nrows*ncols),
    entries : rand_selection(entries, nrows*ncols),
    count : 1,
    z : zeromatrix(nrows, ncols),
    for i : 1 thru nrows do (
        for j : 1 thru ncols do (
            z[i][j] : entries[count],
            count : count + 1
            )
        ),
    z
);
/* where does mx A differ from mx B?  return a list of coordinates where
 * they differ e.g. [[1, 2], [5, 5]] */
where_differ(A, B) := block([i, j, differs],
    differs : [],
    for i : 1 thru matrix_size(A)[1] do (
        for j : 1 thru matrix_size(A)[2] do (
            if not is(A[i][j] = B[i][j]) then (
                differs : cons([i, j], differs)
            )
        )
    ),
    differs
);

/* convert output from where_differ to a printable string */
differences_string(li) := block([output, i],
    output : "",
    for i : 1 thru length(li) - 1 do (
        output : sconcat(output, "(", li[i][1], ", ", li[i][2], "), ")
        ),
    output : sconcat(output,"(", li[length(li)][1], ", ", li[length(li)][2], ")"),
    output 
);

/* does the matrix product AB make sense? */
can_multiply(A, B) := is(matrix_size(A)[2] = matrix_size(B)[1]);

/* produce a string showing how to calculate the correct value of the
 * i,j entry of AB */
work_out_product_entry(A, B, i, j) := block([output, k],
    output : "",
    for k : 1 thru matrix_size(A)[2] - 1 do(
        output : sconcat(output, "(", A[i][k], "\\times ", B[k][j], ") + ")
        ),
    output : sconcat(output, "(", A[i][matrix_size(A)[2]], "\\times ", B[matrix_size(A)[2]][j], ")"),
    output
);

/* return a random linear combination of the first r rows of A */
random_lc(r, A, scale) := block([output, i],
    output : makelist(0, i, 1, matrix_size(A)[2]),
    for i : 1 thru r do (
        output : output + rand_with_step(-scale, scale, 1) * A[i]
        ),
    output
);

/* produce a random m by n matrix with rank r, assumed <= min(m,n). maxi controls
 * size of entries to some extent */
random_with_rank(m, n, r, maxi) := block([larger, smaller, li, z, rw, i, row_order],
    larger : max(m, n),
    smaller : min(m, n),
    li : random_slnz(larger, maxi),
    z : zeromatrix(smaller, larger), /* rows are easier to work with*/
    /* put the first r rows of li into random positions in z, the
     * eventual output matrix. They're LI, so this will guarantee rank at
     * least r
     */
    row_order : rand_selection(makelist(i, i, 1, smaller), smaller),
    for rw : 1 thru r do (
        z[row_order[rw]] : li[rw]
        ),
    /* fill the rest of z up with linear combinations of the rows
     * already added, guaranteeing the rank is exactly r
     */
    for rw : r + 1 thru smaller do (
        z[row_order[rw]] : random_lc(r, li, 3)
        ),
    if (m > smaller) then (z : transpose(z)),
    z
);

random_with_rank2(m, n, r) := block([i, z],
    z : zeromatrix(m, n),
    for i : 1 thru r do (
        z[i][i] : 1
        ),
    random_slnz(m, 3) . z . random_slnz(n, 3)
);

/* multiply a rational matrix mx by the lcm of the denominators of its
 * entries, producing an integer matrix */
clear_denoms(mx) := block([i, denoms],
    denoms : [],
    for i : 1 thru matrix_size(mx)[1] do (
        for j : 1 thru matrix_size(mx)[2] do (
            denoms : cons(denom(mx[i][j]), denoms)
            )
        ),
    lcm(denoms) * mx
);

/* factor out the gcd of the entries of an integer matrix */
divide_out_gcd(mx) := block([i, j, entries_gcd],
    entries_gcd : mx[1][1],
    for i : 1 thru matrix_size(mx)[1] do (
        for j : 1 thru matrix_size(mx)[2] do (
            entries_gcd : gcd(entries_gcd, mx[i][j])
        )
    ),
    (1 / entries_gcd) * mx
);

/* get an integer vector not in the column space of A. Assumes A is a
 * rational matrix that doesn't have full rank */
vector_not_in_image(A) := block([pA, pAc, basis_for_image_complement],
    pA : proj(A),
    pAc : ident(matrix_size(A)[1]) - pA,
    basis_for_image_complement : clear_denoms(refine(pAc)),
    col(basis_for_image_complement, 1)
);

/* get a small nonzero integer vector in the column space of A. maxi controls
 * how small */
small_vector_in_image(A, maxi) := block([pA, output, i],
    pA : refine(proj(A)),
    output : rand_with_prohib(-maxi, maxi, [0]) * col(pA, 1),
    for i : 2 thru matrix_size(pA)[2] do (
        output : output + rand_with_step(-maxi, maxi, 1) * col(pA, i)
        ),
    output
);

image_vector_and_preimage(A, maxi) := block([],
    preim : rand_mx(matrix_size(A)[2], 1, maxi),
    im : A . preim,
    [im, preim]
);