Introduction to linearalgebra
=============================

linearalgebra is a collection of functions for linear algebra.

Example:

     (%i1) M : matrix ([1, 2], [1, 2]);
                                 [ 1  2 ]
     (%o1)                       [      ]
                                 [ 1  2 ]
     (%i2) nullspace (M);
                                    [  1  ]
                                    [     ]
     (%o2)                     span([   1 ])
                                    [ - - ]
                                    [   2 ]
     (%i3) columnspace (M);
                                     [ 1 ]
     (%o3)                      span([   ])
                                     [ 1 ]
     (%i4) ptriangularize (M - z*ident(2), z);
                              [ 1   2 - z   ]
     (%o4)                    [             ]
                              [           2 ]
                              [ 0  3 z - z  ]
     (%i5) M : matrix ([1, 2, 3], [4, 5, 6], [7, 8, 9]) - z*ident(3);
                          [ 1 - z    2      3   ]
                          [                     ]
     (%o5)                [   4    5 - z    6   ]
                          [                     ]
                          [   7      8    9 - z ]
     (%i6) MM : ptriangularize (M, z);
                   [ 4  5 - z            6            ]
                   [                                  ]
                   [                2                 ]
                   [     66        z    102 z   132   ]
                   [ 0   --      - -- + ----- + ---   ]
     (%o6)         [     49        7     49     49    ]
                   [                                  ]
                   [               3        2         ]
                   [           49 z    245 z    147 z ]
                   [ 0    0    ----- - ------ - ----- ]
                   [            264      88      44   ]
     (%i7) algebraic : true;
     (%o7)                         true
     (%i8) tellrat (MM [3, 3]);
                              3       2
     (%o8)                  [z  - 15 z  - 18 z]
     (%i9) MM : ratsimp (MM);
                    [ 4  5 - z           6           ]
                    [                                ]
                    [                2               ]
     (%o9)          [     66      7 z  - 102 z - 132 ]
                    [ 0   --    - ------------------ ]
                    [     49              49         ]
                    [                                ]
                    [ 0    0             0           ]
     (%i10) nullspace (MM);
                             [        1         ]
                             [                  ]
                             [   2              ]
                             [  z  - 14 z - 16  ]
                             [  --------------  ]
     (%o10)             span([        8         ])
                             [                  ]
                             [    2             ]
                             [   z  - 18 z - 12 ]
                             [ - -------------- ]
                             [         12       ]
     (%i11) M : matrix ([1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12],
                        [13, 14, 15, 16]);
                            [ 1   2   3   4  ]
                            [                ]
                            [ 5   6   7   8  ]
     (%o11)                 [                ]
                            [ 9   10  11  12 ]
                            [                ]
                            [ 13  14  15  16 ]
     (%i12) columnspace (M);
                                [ 1  ]  [ 2  ]
                                [    ]  [    ]
                                [ 5  ]  [ 6  ]
     (%o12)                span([    ], [    ])
                                [ 9  ]  [ 10 ]
                                [    ]  [    ]
                                [ 13 ]  [ 14 ]
     (%i13) apply ('orthogonal_complement, args (nullspace (transpose (M))));
                                [ 0 ]  [  1  ]
                                [   ]  [     ]
                                [ 1 ]  [  0  ]
     (%o13)                span([   ], [     ])
                                [ 2 ]  [ - 1 ]
                                [   ]  [     ]
                                [ 3 ]  [ - 2 ]

Functions and Variables for linearalgebra
=========================================

 -- Function: addmatrices (F, M_1, ..., M_N)

     Using the function F as the addition function, return the sum of the matrices M_1, ..., M_N. The function F must accept any number of arguments (a Maxima nary function).

     Examples:

          (%i1) m1 : matrix([1,2],[3,4])$
          (%i2) m2 : matrix([7,8],[9,10])$
          (%i3) addmatrices('max,m1,m2);
          (%o3) matrix([7,8],[9,10])
          (%i4) addmatrices('max,m1,m2,5*m1);
          (%o4) matrix([7,10],[15,20])

 -- Function: blockmatrixp (M)

     Return true if and only if M is a matrix and every entry of M is a matrix.

 -- Function: columnop (M, I, J, THETA)

     If M is a matrix, return the matrix that results from doing the column operation C_i <- C_i - THETA * C_j. If M doesn't have a row I or J, signal an error.

 -- Function: columnswap (M, I, J)

     If M is a matrix, swap columns I and J.  If M doesn't have a column I or J, signal an error.

 -- Function: columnspace (M)

     If M is a matrix, return span (v_1, ..., v_n), where the set {v_1, ..., v_n} is a basis for the column space of M.  The span of the empty set is {0}. Thus, when the column space has only one member, return span ().

 -- Function: copy (E)

     Return a copy of the Maxima expression E. Although E can be any Maxima expression, the copy function is the most useful when E is either a list or a matrix; consider:

          (%i1) m : [1,[2,3]]$
          (%i2) mm : m$
          (%i3) mm[2][1] : x$
          (%i4) m;
          (%o4)                      [1,[x,3]]
          (%i5) mm;
          (%o5)                      [1,[x,3]]

     Let's try the same experiment, but this time let MM be a copy of M

          (%i6) m : [1,[2,3]]$
          (%i7) mm : copy(m)$
          (%i8) mm[2][1] : x$
          (%i9) m;
          (%o9)                     [1,[2,3]]
          (%i10) mm;
          (%o10)                    [1,[x,3]]

     This time, the assignment to MM does not change the value of M.

 -- Function: cholesky (M)
 -- Function: cholesky (M, FIELD)

     Return the Cholesky factorization of the matrix selfadjoint (or hermitian) matrix M. The second argument defaults to 'generalring.' For a description of the possible values for FIELD, see lu_factor.

 -- Function: ctranspose (M)

     Return the complex conjugate transpose of the matrix M. The function ctranspose uses matrix_element_transpose to transpose each matrix element.

 -- Function: diag_matrix (D_1, D_2,...,D_N)

     Return a diagonal matrix with diagonal entries D_1, D_2,...,D_N. When the diagonal entries are matrices, the zero entries of the returned matrix are zero matrices of the appropriate size; for example:

          (%i1) diag_matrix(diag_matrix(1,2),diag_matrix(3,4));

                                      [ [ 1  0 ]  [ 0  0 ] ]
                                      [ [      ]  [      ] ]
                                      [ [ 0  2 ]  [ 0  0 ] ]
          (%o1)                       [                    ]
                                      [ [ 0  0 ]  [ 3  0 ] ]
                                      [ [      ]  [      ] ]
                                      [ [ 0  0 ]  [ 0  4 ] ]
          (%i2) diag_matrix(p,q);

                                             [ p  0 ]
          (%o2)                              [      ]
                                             [ 0  q ]

 -- Function: dotproduct (U, V)

     Return the dotproduct of vectors U and V.  This is the same as conjugate (transpose (U)) . V.  The arguments U and V must be column vectors.

 -- Function: eigens_by_jacobi (A)
 -- Function: eigens_by_jacobi (A, FIELD_TYPE)

     Computes the eigenvalues and eigenvectors of A by the method of Jacobi rotations. A must be a symmetric matrix (but it need not be positive definite nor positive semidefinite). FIELD_TYPE indicates the computational field, either floatfield or bigfloatfield. If FIELD_TYPE is not specified, it defaults to floatfield.

     The elements of A must be numbers or expressions which evaluate to numbers via float or bfloat (depending on FIELD_TYPE).

     Examples:

          (%i1) S: matrix([1/sqrt(2), 1/sqrt(2)],[-1/sqrt(2), 1/sqrt(2)]);
                               [     1         1    ]
                               [  -------   ------- ]
                               [  sqrt(2)   sqrt(2) ]
          (%o1)                [                    ]
                               [      1        1    ]
                               [ - -------  ------- ]
                               [   sqrt(2)  sqrt(2) ]
          (%i2) L : matrix ([sqrt(3), 0], [0, sqrt(5)]);
                                [ sqrt(3)     0    ]
          (%o2)                 [                  ]
                                [    0     sqrt(5) ]
          (%i3) M : S . L . transpose (S);
                      [ sqrt(5)   sqrt(3)  sqrt(5)   sqrt(3) ]
                      [ ------- + -------  ------- - ------- ]
                      [    2         2        2         2    ]
          (%o3)       [                                      ]
                      [ sqrt(5)   sqrt(3)  sqrt(5)   sqrt(3) ]
                      [ ------- - -------  ------- + ------- ]
                      [    2         2        2         2    ]
          (%i4) eigens_by_jacobi (M);
          The largest percent change was 0.1454972243679
          The largest percent change was 0.0
          number of sweeps: 2
          number of rotations: 1
          (%o4) [[1.732050807568877, 2.23606797749979],
                                  [  0.70710678118655   0.70710678118655 ]
                                  [                                      ]]
                                  [ - 0.70710678118655  0.70710678118655 ]
          (%i5) float ([[sqrt(3), sqrt(5)], S]);
          (%o5) [[1.732050807568877, 2.23606797749979],
                                  [  0.70710678118655   0.70710678118655 ]
                                  [                                      ]]
                                  [ - 0.70710678118655  0.70710678118655 ]
          (%i6) eigens_by_jacobi (M, bigfloatfield);
          The largest percent change was 1.454972243679028b-1
          The largest percent change was 0.0b0
          number of sweeps: 2
          number of rotations: 1
          (%o6) [[1.732050807568877b0, 2.23606797749979b0],
                          [  7.071067811865475b-1   7.071067811865475b-1 ]
                          [                                              ]]
                          [ - 7.071067811865475b-1  7.071067811865475b-1 ]

 -- Function: get_lu_factors (X)

     When X = lu_factor (A), then get_lu_factors returns a list of the form [P, L, U], where P is a permutation matrix, L is lower triangular with ones on the diagonal, and U is upper triangular, and A = P L U.

 -- Function: hankel (COL)
 -- Function: hankel (COL, ROW)

     Return a Hankel matrix H. The first column of H is COL; except for the first entry, the last row of H is ROW. The default for ROW is the zero vector with the same length as COL.

 -- Function: hessian (F,VARS)

     Return the hessian matrix of F with respect to the variables in the list VARS.  The I,J entry of the hessian matrix is DIFF(F VARS[I],1,VARS[J],1).

 -- Function: hilbert_matrix (N)

     Return the N by N Hilbert matrix. When N isn't a positive integer, signal an error.

 -- Function: identfor (M)
 -- Function: identfor (M, FLD)

     Return an identity matrix that has the same shape as the matrix M.  The diagonal entries of the identity matrix are the multiplicative identity of the field FLD; the default for FLD is GENERALRING.

     The first argument M should be a square matrix or a non-matrix. When M is a matrix, each entry of M can be a square matrix -- thus M can be a blocked Maxima matrix. The matrix can be blocked to any (finite) depth.

     See also zerofor

 -- Function: invert_by_lu (M, (RNG GENERALRING))

     Invert a matrix M by using the LU factorization.  The LU factorization is done using the ring RNG.

 -- Function: kronecker_product (A, B)

     Return the Kronecker product of the matrices A and B.

 -- Function: listp (E, P)
 -- Function: listp (E)

     Given an optional argument P, return true if E is a Maxima list and P evaluates to true for every list element. When listp is not given the optional argument, return true if E is a Maxima list.  In all other cases, return false.

 -- Function: locate_matrix_entry (M, R_1, C_1, R_2, C_2, F, REL)

     The first argument must be a matrix; the arguments R_1 through C_2 determine a sub-matrix of M that consists of rows R_1 through R_2 and columns C_1 through C_2.

     Find a entry in the sub-matrix M that satisfies some property. Three cases:

     (1) REL = 'bool and F a predicate:

     Scan the sub-matrix from left to right then top to bottom, and return the index of the first entry that satisfies the predicate F. If no matrix entry satisfies F, return false.

     (2) REL = 'max and F real-valued:

     Scan the sub-matrix looking for an entry that maximizes F. Return the index of a maximizing entry.

     (3) REL = 'min and F real-valued:

     Scan the sub-matrix looking for an entry that minimizes F. Return the index of a minimizing entry.

 -- Function: lu_backsub (M, B)

     When M = lu_factor (A, FIELD), then lu_backsub (M, B) solves the linear system A X = B.

 -- Function: lu_factor (M, FIELD)

     Return a list of the form [LU, PERM, FLD], or [LU, PERM, FLD, LOWER-CND UPPER-CND], where

     (1) The matrix LU contains the factorization of M in a packed form. Packed form means three things: First, the rows of LU are permuted according to the list PERM.  If, for example, PERM is the list [3,2,1], the actual first row of the LU factorization is the third row of the matrix LU. Second, the lower triangular factor of m is the lower triangular part of LU with the diagonal entries replaced by all ones. Third, the upper triangular factor of M is the upper triangular part of LU.

     (2) When the field is either floatfield or complexfield, the numbers LOWER-CND and UPPER-CND are lower and upper bounds for the infinity norm condition number of M.  For all fields, the condition number might not be estimated; for such fields, lu_factor returns a two item list. Both the lower and upper bounds can differ from their true values by arbitrarily large factors. (See also mat_cond.)

     The argument M must be a square matrix.

     The optional argument FLD must be a symbol that determines a ring or field. The pre-defined fields and rings are:

     (a) generalring -- the ring of Maxima expressions, (b) floatfield --  the field of floating point numbers of the type double, (c) complexfield --  the field of complex floating point numbers of the type double, (d) crering  -- the ring of Maxima CRE expressions, (e) rationalfield -- the field of rational numbers, (f) runningerror -- track the all floating point rounding errors, (g) noncommutingring -- the ring of Maxima expressions where multiplication is the non-commutative dot operator.

     When the field is floatfield, complexfield, or runningerror, the algorithm uses partial pivoting; for all other fields, rows are switched only when needed to avoid a zero pivot.

     Floating point addition arithmetic isn't associative, so the meaning of 'field' differs from the mathematical definition.

     A member of the field runningerror is a two member Maxima list of the form [x,n],where X is a floating point number and n is an integer. The relative difference between the 'true' value of x and x is approximately bounded by the machine epsilon times n. The running error bound drops some terms that of the order the square of the machine epsilon.

     There is no user-interface for defining a new field. A user that is familiar with Common Lisp should be able to define a new field.  To do this, a user must define functions for the arithmetic operations and functions for converting from the field representation to Maxima and back. Additionally, for ordered fields (where partial pivoting will be used), a user must define functions for the magnitude and for comparing field members.  After that all that remains is to define a Common Lisp structure mring.  The file mring has many examples.

     To compute the factorization, the first task is to convert each matrix entry to a member of the indicated field. When conversion isn't possible, the factorization halts with an error message. Members of the field needn't be Maxima expressions.  Members of the complexfield, for example, are Common Lisp complex numbers. Thus after computing the factorization, the matrix entries must be converted to Maxima expressions.

     See also  get_lu_factors.

     Examples:

          (%i1) w[i,j] := random (1.0) + %i * random (1.0);
          (%o1)          w     := random(1.) + %i random(1.)
                          i, j
          (%i2) showtime : true$
          Evaluation took 0.00 seconds (0.00 elapsed)
          (%i3) M : genmatrix (w, 100, 100)$
          Evaluation took 7.40 seconds (8.23 elapsed)
          (%i4) lu_factor (M, complexfield)$
          Evaluation took 28.71 seconds (35.00 elapsed)
          (%i5) lu_factor (M, generalring)$
          Evaluation took 109.24 seconds (152.10 elapsed)
          (%i6) showtime : false$

          (%i7) M : matrix ([1 - z, 3], [3, 8 - z]);
                                  [ 1 - z    3   ]
          (%o7)                   [              ]
                                  [   3    8 - z ]
          (%i8) lu_factor (M, generalring);
                    [ 1 - z         3        ]
                    [                        ]
          (%o8)    [[   3            9       ], [1, 2], generalring]
                    [ -----  - z - ----- + 8 ]
                    [ 1 - z        1 - z     ]
          (%i9) get_lu_factors (%);
                            [   1    0 ]  [ 1 - z         3        ]
                  [ 1  0 ]  [          ]  [                        ]
          (%o9)  [[      ], [   3      ], [                9       ]]
                  [ 0  1 ]  [ -----  1 ]  [   0    - z - ----- + 8 ]
                            [ 1 - z    ]  [              1 - z     ]
          (%i10) %[1] . %[2] . %[3];
                                  [ 1 - z    3   ]
          (%o10)                  [              ]
                                  [   3    8 - z ]

 -- Function: mat_cond (M, 1)
 -- Function: mat_cond (M, inf)

     Return the P-norm matrix condition number of the matrix M. The allowed values for P are 1 and INF.  This function uses the LU factorization to invert the matrix M. Thus the running time for mat_cond is proportional to the cube of the matrix size; lu_factor determines lower and upper bounds for the infinity norm condition number in time proportional to the square of the matrix size.

 -- Function: mat_norm (M, 1)
 -- Function: mat_norm (M, inf)
 -- Function: mat_norm (M, frobenius)

     Return the matrix P-norm of the matrix M.  The allowed values for P are 1, inf, and frobenius (the Frobenius matrix norm). The matrix M should be an unblocked matrix.

 -- Function: matrixp (E, P)
 -- Function: matrixp (E)

     Given an optional argument P, return true if E is a matrix and P evaluates to true for every matrix element. When matrixp is not given an optional argument, return true if e is a matrix.  In all other cases, return false.

     See also blockmatrixp

 -- Function: matrix_size (M)

     Return a two member list that gives the number of rows and columns, respectively of the matrix M.

 -- Function: mat_fullunblocker (M)

     If M is a block matrix, unblock the matrix to all levels. If M is a matrix, return M; otherwise, signal an error.

 -- Function: mat_trace (M)

     Return the trace of the matrix M. If M isn't a matrix, return a noun form. When M is a block matrix, mat_trace(M) returns the same value as does mat_trace(mat_unblocker(m)).

 -- Function: mat_unblocker (M)

     If M is a block matrix, unblock M one level. If M is a matrix, mat_unblocker (M) returns M; otherwise, signal an error.

     Thus if each entry of M is matrix, mat_unblocker (M) returns an unblocked matrix, but if each entry of M is a block matrix, mat_unblocker (M) returns a block matrix with one less level of blocking.

     If you use block matrices, most likely you'll want to set matrix_element_mult to "." and matrix_element_transpose to 'transpose. See also mat_fullunblocker.

     Example:

          (%i1) A : matrix ([1, 2], [3, 4]);
                                      [ 1  2 ]
          (%o1)                       [      ]
                                      [ 3  4 ]
          (%i2) B : matrix ([7, 8], [9, 10]);
                                      [ 7  8  ]
          (%o2)                       [       ]
                                      [ 9  10 ]
          (%i3) matrix ([A, B]);
                               [ [ 1  2 ]  [ 7  8  ] ]
          (%o3)                [ [      ]  [       ] ]
                               [ [ 3  4 ]  [ 9  10 ] ]
          (%i4) mat_unblocker (%);
                                   [ 1  2  7  8  ]
          (%o4)                    [             ]
                                   [ 3  4  9  10 ]

 -- Function: nonnegintegerp (N)

     Return true if and only if N >= 0 and N is an integer.

 -- Function: nullspace (M)

     If M is a matrix, return span (v_1, ..., v_n), where the set {v_1, ..., v_n} is a basis for the nullspace of M.  The span of the empty set is {0}. Thus, when the nullspace has only one member, return span ().

 -- Function: nullity (M)

     If M is a matrix, return the dimension of the nullspace of M.

 -- Function: orthogonal_complement (V_1, ..., V_N)

     Return span (u_1, ..., u_m), where the set {u_1, ..., u_m} is a basis for the orthogonal complement of the set (v_1, ..., v_n).

     Each vector V_1 through V_N must be a column vector.

 -- Function: polynomialp (P, L, COEFFP, EXPONP)
 -- Function: polynomialp (P, L, COEFFP)
 -- Function: polynomialp (P, L)

     Return true if P is a polynomial in the variables in the list L, The predicate COEFFP must evaluate to true for each coefficient, and the predicate EXPONP must evaluate to true for all exponents of the variables in L. If you want to use a non-default value for EXPONP, you must supply COEFFP with a value even if you want to use the default for COEFFP.

     polynomialp (P, L, COEFFP) is equivalent to polynomialp (P, L, COEFFP, 'nonnegintegerp).

     polynomialp (P, L) is equivalent to polynomialp (P, L, 'constantp, 'nonnegintegerp).

     The polynomial needn't be expanded:

          (%i1) polynomialp ((x + 1)*(x + 2), [x]);
          (%o1)                         true
          (%i2) polynomialp ((x + 1)*(x + 2)^a, [x]);
          (%o2)                         false

     An example using non-default values for coeffp and exponp:

          (%i1) polynomialp ((x + 1)*(x + 2)^(3/2), [x], numberp, numberp);
          (%o1)                         true
          (%i2) polynomialp ((x^(1/2) + 1)*(x + 2)^(3/2), [x], numberp,
                                                                  numberp);
          (%o2)                         true

     Polynomials with two variables:

          (%i1) polynomialp (x^2 + 5*x*y + y^2, [x]);
          (%o1)                         false
          (%i2) polynomialp (x^2 + 5*x*y + y^2, [x, y]);
          (%o2)                         true

 -- Function: polytocompanion (P, X)

     If P is a polynomial in X, return the companion matrix of P. For a monic polynomial P of degree N, we have P = (-1)^N charpoly (polytocompanion (P, X)).

     When P isn't a polynomial in X, signal an error.

 -- Function: ptriangularize (M, V)

     If M is a matrix with each entry a polynomial in V, return a matrix M2 such that

     (1) M2 is upper triangular,

     (2) M2 = E_N ... E_1 M, where E_1 through E_N are elementary matrices whose entries are polynomials in V,

     (3) |det (M)| = |det (M2)|,

     Note: This function doesn't check that every entry is a polynomial in V.

 -- Function: rowop (M, I, J, THETA)

     If M is a matrix, return the matrix that results from doing the row operation R_i <- R_i - theta * R_j. If M doesn't have a row I or J, signal an error.

 -- Function: rank (M)

     Return the rank of that matrix M. The rank is the dimension of the column space. Example:

          (%i1) rank(matrix([1,2],[2,4]));
          (%o1)                                  1
          (%i2) rank(matrix([1,b],[c,d]));
          Proviso:  {d - b c # 0}
          (%o2)                                  2

 -- Function: rowswap (M, I, J)

     If M is a matrix, swap rows I and J. If M doesn't have a row I or J, signal an error.

 -- Function: toeplitz (COL)
 -- Function: toeplitz (COL, ROW)

     Return a Toeplitz matrix T. The first first column of T is COL; except for the first entry, the first row of T is ROW. The default for ROW is complex conjugate of COL. Example:

          (%i1)  toeplitz([1,2,3],[x,y,z]);

                                            [ 1  y  z ]
                                            [         ]
          (%o1)                             [ 2  1  y ]
                                            [         ]
                                            [ 3  2  1 ]
          (%i2)  toeplitz([1,1+%i]);

                                        [   1     1 - %I ]
          (%o2)                         [                ]
                                        [ %I + 1    1    ]

 -- Function: vandermonde_matrix ([X_1, ..., X_N])

     Return a N by N matrix whose I-th row is [1, X_I, X_I^2, ... X_I^(N-1)].

 -- Function: zerofor (M)
 -- Function: zerofor (M, FLD)

     Return a zero  matrix that has the same shape as the matrix M.  Every entry of the zero matrix is the additive identity of the field FLD; the default for FLD is GENERALRING.

     The first argument M should be a square matrix or a non-matrix. When M is a matrix, each entry of M can be a square matrix -- thus M can be a blocked Maxima matrix. The matrix can be blocked to any (finite) depth.

     See also identfor

 -- Function: zeromatrixp (M)

     If M is not a block matrix, return true if is (equal (E, 0)) is true for each element E of the matrix M.  If M is a block matrix, return true if zeromatrixp evaluates to true for each element of E.