Linear algebra:

 

§         Eigenvalues/Eigenvectors, Decompositions, and PseudoInverse

§         Special Matrices

§         Other functions

 

 

w        Eigenvalues/Eigenvectors, Decompositions, and PseudoInverse

§         Cholesky(mat) computes the Cholesky decomposition of a square symmetric positive-definite matrix mat and returns a lower-triangular matrix L such that mat=L×LT
Example: Cholesky(Hilbert(3)) Þ [[1,0,0][1/2,Ö(3)/6,0][1/3,Ö(3)/6,Ö(5)/30]]
Notes: Cholesky returns a lower-triangular matrix instead of an upper-triangular one, because trying to transpose the lower-triangular matrix results in conj() being applied to all symbolic elements, which looks ugly. (PowerExp can be used to remove most but not all of the conj()'s.) Complex elements are not allowed, since a symmetric positive-definite matrix must be real.

§         Eigenval(matrix) returns the eigenvalues of matrix, including repeated ones
Needs: mZeros, Sort
Example: Eigenval([[1,1,-1][0,0,2][0,-1,3]]) Þ {1,1,2}
Note: Eigenval tries to return exact eigenvalues whenever possible, but it calls the built-in eigVl() function for approximate matrices. Non-symmetric matrices can have eigenvalues with nonzero imaginary parts.

§         Eigenvec(matrix) returns the generalized eigenvectors of matrix
Needs: DelElem, Eigenval, list2eqn, MemberQ, NullVecs
Example: Eigenvec([[1,1,-1][0,0,2][0,-1,3]]) Þ [[1,@1,0][0,2,1][0,1,1]]

§         GenEigVl(mat1,mat2) returns the generalized eigenvalues of the matrix pencil mat1-l×mat2
Needs: Coef, Cubic, Degree, NullVecs
Examples: GenEigVl([[3,6][1,-4]],[[-10,6][0,10]]) Þ
                                                               {-8/25-Ö(194)/50×i, -8/25+Ö(194)/50×i},
                  GenEigVl([[1,0,0,0][0,2,0,0][0,0,3,0][0,0,0,4]],
                                         [[1,0,0,0][0,1,0,0][0,0,0,0][0,0,0,0]]) Þ {¥, ¥, 1, 2},
                  GenEigVl(diag(seq(i,i,1,3)), identity(3)) Þ {1, 3, 2},
                  GenEigVl([[0,0][0,0]],[[0,0][0,0]]) Þ {undef, undef},
                  GenEigVl([[0,0][0,1]],[[0,0][0,0]]) Þ
                                         "Error: Singular intermediate matrix P = [[1,0][0,0]]"
Notes: Generalized eigenvalues can be infinite. The generalized eigenvalue problem is called regular if the determinant det(mat1-l×mat2) does not vanish identically for all values of l, otherwise it is called singular. There are n generalized eigenvalues iff the matrix rank of mat2 is n. Application areas for generalized eigenproblems include optimization, statistics, signal processing, control theory, spectral imaging, and solid mechanics.

§         GSchmidt(mat) performs Gram-Schmidt orthonormalization on mat
Example: GSchmidt([[1,2][3,4]]) Þ [[Ö(10)/10,3×Ö(10)/10]
                                                            [3×Ö(10)/10,-Ö(10)/10]]
Note: GSchmidt returns a matrix with orthogonal columns that span the same subspace as the columns of mat, provided the columns of mat are linearly independent. Gram-Schmidt orthonormalization is known to be numerically unstable (see “Matrix Computations” by Golub and van Loan), but should be fine for symbolic manipulation.

§         Jordan(matrix) returns the Jordan decomposition {J,S} of matrix such that J=S-1×matrix×S
Needs: Eigenvec
Example: Jordan([[1,2][3,4]]) Þ {j=[[5/2-Ö(33)/2,0][0,Ö(33)/2+5/2]],
                                                       s=[[-(Ö(33)+3)/6,(Ö(33)-3)/6][1,1]]}
Note: Jordan decomposition is often used for evaluating symbolic functions of matrices, such as matrix exponentials.

§         PInverse(mat,tol) returns the pseudoinverse (also called the generalized or Moore-Penrose inverse) of the matrix mat with tolerance tol
Needs: SVD
Example: PInverse([[1,2][3,4][5,6]],1E-10) Þ [[-1.33,-0.333,0.667]
                                                                             [1.08,0.333,-0.417]]
                                                                             (to Float-3 precision)
Note: mat can be non-square. The pseudoinverse equals the inverse when mat is a non-singular (all singular values > 0) square matrix. The sum of squares of the elements of mat×PInverse(mat)-I is minimized, where I is the identity matrix, and so the pseudoinverse is useful for data fitting.

§         SingVals(mat) returns the sorted singular values of the numeric matrix mat
Needs: IsCmplxN, Sort
Examples: SingVals([[i,7][-2,1-i]]) Þ {2.09236143579,7.18484680575},
                  SingVals([[1,0][2,0]]) Þ {0.,2.2360679775}

§         SVD(mat) computes the singular value decomposition of numeric matrix mat
Needs: GSchmidt, Reverse
Example: SVD([[1,2][3,4][5,6]]) Þ {u=[[-.883,-.23][-.241,-.525][.402,-.82]],
                                sigma=[[.514,0][0,9.53][0,0]], v=[[.785,-.62][-.62,-.785]]}
                                                                                       (with Float-3 precision)
Note: mat must be a numeric matrix (no undefined variables). It can be non-square. The sigma matrix contains the diagonalized singular values of mat.

w        Special Matrices

§         Circulnt(col) returns the circulant matrix defined by the column col
Example: Circulnt({1,2,3,4}) Þ [[1,4,3,2][2,1,4,3][3,2,1,4][4,3,2,1]]
Note: A circulant matrix is a special case of the Toeplitz matrix.

§         Frank(n) returns the Frank matrix of order n
Example: Frank(4) Þ [[4,3,2,1][3,3,2,1][0,2,2,1][0,0,1,1]]

§         Hadamard(n) returns a Hadamard matrix of order n
Needs: KProduct
Example: Hadamard(4) Þ [[1,1,1,1][1,-1,1,-1][1,1,-1,-1][1,-1,-1,1]]
Notes: Hadamard matrices exist only for n = 2 or n a multiple of 4. For an n´n Hadamard matrix H, HT×H = n×I, where I is the identity matrix, and det(H) = nn/2. Hadamard matrices are used in various combinatorics problems such as the Reed-Muller error-correcting codes.

§         Hankel(list) returns the corresponding Hankel matrix
Example: Hankel({1,1,1,6,21}) Þ [[1,1,1][1,1,6][1,6,21]]
Note: Hankel matrices are symmetric.

§         Hilbert(n) returns the Hilbert matrix of order n
Example: Hilbert(3) Þ [[1,1/2,1/3][1/2,1/3,1/4][1/3,1/4,1/5]]
Note: The inverse of a Hilbert matrix has integer elements. A Hilbert matrix is a special case of the Hankel matrix.

§         IHilbert(n) returns the inverse of the order-n Hilbert matrix
Needs: Hilbert
Example: IHilbert(3) Þ [[9,-36,30][-36,192,-180][30,-180,180]]
Note: IHilbert(n) is usually faster than Hilbert(n)-1. For example, Hilbert(10)-1 takes about 17.2 seconds while IHilbert(10) takes about 4.0 seconds.

§         Pascal(n) returns the Pascal matrix of order n
Example: Pascal(4) Þ [[1,1,1,1][1,2,3,4][1,3,6,10][1,4,10,20]]

§         Toeplitz(row,col) returns the corresponding Toeplitz matrix
Needs: Reverse
Example: Toeplitz({a,b,c},{a,d,e}) Þ [[a,b,c][d,a,b][e,d,a]]
Note: row[1] must equal col[1].

§         Vanderm(list) returns the corresponding Vandermonde matrix
Example: Vanderm({v1,v2,v3}) Þ [[1,1,1][v1,v2,v3][v12,v22,v32]]
Note: Vandermonde matrices appear in association with fitting polynomials to data. In particular, you can create a Vandermonde matrix from the x-values, augment the y-values to the bottom, perform an rref of the transpose, and obtain an interpolating polynomial by doing the sum xk-1×mat[k,colDim(mat)], k=1, …, rowDim(mat), where mat is the row-reduced matrix.

w        Other functions

§         cCondNum(matrix) returns the column-norm condition number of matrix
Example: cCondNum(Hilbert(4)) Þ 28375
Note: Condition numbers can be useful for studying the propagation of errors. Higher values for the condition number indicate that the matrix is close to singular, which means that approximate numerical operations like inversion or linear solution will return inaccurate or unreliable results. The base-b logarithm of the matrix condition number is a worst-case estimate of how many base-b digits are lost in solving a linear system with that matrix (see MathWorld entry).

§         Cofactor(mat) returns the matrix of cofactors of mat
Needs: Minor
Example: Cofactor([[a,b][c,d]]) Þ [[d,-c][-b,a]]
Notes: The adjoint of the matrix is simply the transpose of the matrix of cofactors. The determinant of the matrix mati j can be computed from its matrix of cofactors Ci j by multiplying one of the rows or columns of mat by the corresponding row or column of C and summing the elements.

§         CompMat(poly,x) returns the companion matrix of poly
Needs: Degree, Reverse
Example: CompMat(x^2+1,x) Þ [[0,-1][1,0]], eigVl([[0,-1][1,0]]) Þ {i, -i}
Note: The eigenvalues of the companion matrix are the roots of poly, as demonstrated in the example above.

§         CondNum(matrix) returns the Turing condition number estimate
Example: CondNum(Hilbert(4)) Þ 25920
Note: Condition numbers can be useful for studying the propagation of errors. Higher values for the condition number indicate that the matrix is close to singular, which means that approximate numerical operations like inversion or linear solution will return inaccurate or unreliable results. The base-b logarithm of the matrix condition number is a worst-case estimate of how many base-b digits are lost in solving a linear system with that matrix (see MathWorld entry).

§         FastDet(mat) returns the determinant of the matrix mat
Example: det([[a,b,c][d,sin(e),f][g,Ö(h),i]]) Þ a×(sin(e)×i-f×Ö(h))-b×(d×i-f×g)+
                                                        c×(d×Ö(h)-sin(e)×g) [takes about 5.5 seconds],
                FastDet([[a,b,c][d,sin(e),f][g,Ö(h),i]]) Þ a×(sin(e)×i-f×Ö(h))-
                                      b×(d×i-f×g)+c×(d×Ö(h)-sin(e)×g) [takes about 0.6 seconds]
Note: FastDet is usually faster than the built-in det() function, especially for symbolic and exact numeric matrices. The determinant of an orthogonal matrix is ±1. One can compute the absolute value of the determinant of a matrix M from its QR factorization as product(diag(R)), and from its SVD as product(SingVals(M)).

§         Givens(matrix,i,j) does a Givens rotation to zero out the [i,j] element of matrix
Example: Givens([[1][1][1][1]],2,1) Þ [[Ö(2)][0][1][1]] (zero out 2nd element),
                Givens([[Ö(2)][0][1][1]],3,1) Þ [[Ö(3)][0][0][1]] (zero out 3rd elem),
                Givens([[Ö(3)][0][0][1]],4,1) Þ [[2][0][0][0]] (zero out 4th element).
Note: A Givens rotation is similar to a Householder reflection, but a Givens rotation is useful when you want to zero out a particular element (or small set of elements), whereas a Householder reflection is more effective when you want to zero out many elements. Givens rotations are skew-symmetric and orthogonal.

§         KProduct(a,b) returns the Kronecker/matrix/tensor product of matrices a and b
Example: KProduct([[a11,a12][a21,a22]],[[b11,b12]]) Þ [[a11×b11,a11×b12][a12×b11,a12×b12][a21×b11,a21×b12][a22×b11,a22×b12]]
Note: KProduct can be used to create Dirac matrices from quantum field theory. The Kronecker product is distributive with respect to addition as in (A+B)ÄC = AÄC+BÄC, and is associative as in AÄ(BÄC) = (AÄB)ÄC. For the transpose and inverse of a Kronecker product, (AÄB)T = ATÄBT and (AÄB)-1 = A-1ÄB-1. For the determinant and trace, if A is n´n and B is m´m, det(AÄB) = det(A)m×det(B)n, and trace(AÄB) = trace(A)×trace(B).

§         KSum(a,b) returns the Kronecker sum of square matrices a and b
Needs: KProduct
Example: KSum([[a,b][c,d]],[[e,f][g,h]]) Þ [[a+e,f,b,0][g,a+h,0,b]
                                                                        [c,0,d+e,f][0,c,g,d+h]]

§         MatFunc(str,mat) applies the function given in the string str to the matrix mat
Needs: Jordan
Example: MatFunc(“e^#”,[[1,2][3,4]]) Þ [[(1/2-Ö(33)/22)×e^(…)+…][…]]
Notes: Use “#” in the string to denote the input matrix. MatFunc can handle functions of matrices (e^A, sin(A), Ö(A), etc.) on a diagonalizable square matrix.

§         MatRank(mat) returns the matrix rank of mat
Needs: NullVecs
Examples: MatRank([a,b;c,d;e,f]) Þ 2, MatRank([[a,0,b][c,d,0][0,0,e]]) Þ 3,
                  MatRank([[a,0,b][c,d,0][0,0,0]]) Þ 2

§         Minor(mat) returns the matrix of minors for mat
Needs: DelCol, DelRow
Example: Minor([[a,b][c,d]]) Þ [[d,c][b,a]]

§         mMinPoly(mat,var) returns the minimal polynomial of mat in terms of var
Needs: NullVecs
Examples: mMinPoly([[0,0][0,0]],x) Þ x, mMinPoly([[0,1][0,0]],x) Þ x2,
                  mMinPoly([[1,0][0,1]],x) Þ x-1
Note: The minimal polynomial of a matrix M is the lowest degree polynomial such that , where n is the polynomial degree. It divides the characteristic polynomial of M.

§         NullVecs(mat) returns a matrix whose columns form a basis for the nullspace of mat
Examples: NullVecs([[4,2,-4][2,1,-2][-4,-2,4]]) Þ [[-1/2,1][1,0][0,1]],
                  NullVecs({{1,1},{1,0}}) Þ "nullspace={0}"
Note: As shown in the second example above, the null space of a nonsingular matrix mat is trivial, so that there is no nonzero vector v such that mat×v = 0.

§         SparsRnd({m,n},nz) returns an m´n random sparse array with nz nonzero values
Needs: ToSparse
Example: RandSeed 0:SparsRnd({2,3},2) Þ
                 sparse(2,{2,3},{0,1,2},{2,3},{0.733812311248,0.405809641769},0)

§         ToDense(arr) converts a sparse matrix arr to a matrix in the usual form
Examples: ToDense(ToSparse([[4,-2][0,-7][8,0]],a)) Þ [[4,-2][0,-7][8,0]],
                  ToDense(sparse(4,{3,2},{0,2,3,4},{1,2,2,1},{4,-2,-7,8},a))
                  Þ [[4,-2][a,-7][8,a]]

§         ToSparse(mat,pat) converts a matrix in the usual form to a sparse matrix with default element pat
Needs: MatchQ
Example: ToSparse([[4,-2][0,-7][8,0]],0) Þ
                sparse(4,{3,2},{0,2,3,4},{1,2,2,1},{4,-2,-7,8},0)
Note: The format used is the CSR (compressed sparse row) format.

§         VecNorm(vec,p) returns the p-norm of the vector vec
Examples: VecNorm({a,b},2) Þ Ö(a2+b2),
                  VecNorm({a,b},3) Þ (|a3|+|b3|)^(1/3),
                  VecNorm({a,b},¥) Þ max(|b|,|a|)

 

 

MathTools and its documentation copyright Ó Bhuvanesh Bhatt (bbhatt1@towson.edu)