| Up | Next | Prev | PrevTail | Tail |
This package provides a selection of functions that are useful in the world of linear algebra.
This package provides a selection of functions that are useful in the world of linear algebra. These functions are described alphabetically in subsection 20.32.3 and are labelled 20.32.3.1 to 20.32.3.53. They can be classified into four sections(n.b: the numbers after the dots signify the function label in section 20.32.3).
Contributions to this package have been made by Walter Tietze (ZIB).
| add_columns | … | 20.32.3.1 | add_rows | … | 20.32.3.2 |
| add_to_columns | … | 20.32.3.3 | add_to_rows | … | 20.32.3.4 |
| augment_columns | … | 20.32.3.5 | char_poly | … | 20.32.3.9 |
| column_dim | … | 20.32.3.12 | copy_into | … | 20.32.3.14 |
| diagonal | … | 20.32.3.15 | extend | … | 20.32.3.16 |
| find_companion | … | 20.32.3.17 | get_columns | … | 20.32.3.18 |
| get_rows | … | 20.32.3.19 | hermitian_tp | … | 20.32.3.21 |
| matrix_augment | … | 20.32.3.28 | matrix_stack | … | 20.32.3.30 |
| minor | … | 20.32.3.31 | mult_columns | … | 20.32.3.32 |
| mult_rows | … | 20.32.3.33 | pivot | … | 20.32.3.34 |
| remove_columns | … | 20.32.3.37 | remove_rows | … | 20.32.3.38 |
| row_dim | … | 20.32.3.39 | rows_pivot | … | 20.32.3.40 |
| stack_rows | … | 20.32.3.43 | sub_matrix | … | 20.32.3.44 |
| swap_columns | … | 20.32.3.46 | swap_entries | … | 20.32.3.47 |
| swap_rows | … | 20.32.3.48 |
Functions that create matrices.
| band_matrix | … | 20.32.3.6 | block_matrix | … | 20.32.3.7 |
| char_matrix | … | 20.32.3.8 | coeff_matrix | … | 20.32.3.11 |
| companion | … | 20.32.3.13 | hessian | … | 20.32.3.22 |
| hilbert | … | 20.32.3.23 | mat_jacobian | … | 20.32.3.24 |
| jordan_block | … | 20.32.3.25 | make_identity | … | 20.32.3.27 |
| random_matrix | … | 20.32.3.36 | toeplitz | … | 20.32.3.50 |
| Vandermonde | … | 20.32.3.52 | Kronecker_Product | … | 20.32.3.53 |
| char_poly | … | 20.32.3.9 | cholesky | … | 20.32.3.10 |
| gram_schmidt | … | 20.32.3.20 | lu_decom | … | 20.32.3.26 |
| pseudo_inverse | … | 20.32.3.35 | simplex | … | 20.32.3.41 |
| svd | … | 20.32.3.45 | triang_adjoint | … | 20.32.3.51 |
There is a separate NORMFORM package described in section 20.37 for computing the following matrix normal forms in REDUCE:
smithex, smithex_int, frobenius, ratjordan, jordansymbolic, jordan.
| matrixp | … | 20.32.3.29 | squarep | … | 20.32.3.42 |
| symmetricp | … | 20.32.3.49 |
In the examples the matrix 𝒜 will be

Throughout ℐ is used to indicate the identity matrix and 𝒜T to indicate the transpose of the matrix 𝒜.
If you have not used matrices within REDUCE before then the following may be helpful.
Initialisation of matrices takes the following syntax:
mat1 := mat((a,b,c),(d,e,f),(g,h,i));
will produce

The (i,j)th entry can be accessed by:
mat1(i,j);
The package is loaded by:
load_package linalg;
add_columns(𝒜,c1,c2,expr);
| 𝒜 | :- | a matrix. |
| c1, c2 | :- | positive integers. |
| expr | :- | a scalar expression. |
add_columns replaces column c2 of 𝒜 by
expr∗column(𝒜,c1) + column(𝒜,c2).
add_rows performs the equivalent task on the rows of 𝒜.
add_columns
add_rows
add_to_columns, add_to_rows, mult_columns, mult_rows.
See: add_columns.
add_to_columns(𝒜,column_list,expr);
| 𝒜 | :- | a matrix. |
| column_list | :- | a positive integer or a list of positive integers. |
| expr | :- | a scalar expression. |
add_to_columns adds expr to each column specified in column_list of
𝒜.
add_to_rows performs the equivalent task on the rows of 𝒜.
add_to_columns
add_to_rows
add_columns, add_rows, mult_rows, mult_columns.
See: add_to_columns.
augment_columns(𝒜,column_list);
| 𝒜 | :- | a matrix. |
| column_list | :- | either a positive integer or a list of positive integers. |
augment_columns gets hold of the columns of 𝒜 specified in column_list and
sticks them together.
stack_rows performs the same task on rows of 𝒜.
augment_columns
stack_rows
get_columns, get_rows, sub_matrix.
band_matrix(expr_list,square_size);
| expr_list | :- | either a single scalar expression or a list of an odd number of scalar expressions. |
| square_size | :- | a positive integer. |
band_matrix creates a square matrix of dimension square_size. The
diagonal consists of the middle expr of the expr_list. The expressions to the
left of this fill the required number of sub-diagonals and the expressions to the
right the super-diagonals.
band_matrix
diagonal.
block_matrix(r,c,matrix_list);
| r,c | :- | positive integers. |
| matrix_list | :- | a list of matrices. |
block_matrix creates a matrix that consists of r ×c matrices filled from the
matrix_list row-wise.
block_matrix
char_matrix(𝒜,λ);
| 𝒜 | :- | a square matrix. |
| λ | :- | a symbol or algebraic expression. |
char_matrix creates the characteristic matrix 𝒞 of 𝒜. This is 𝒞 = λℐ−𝒜.
char_matrix
char_poly.
char_poly(𝒜,λ);
| 𝒜 | :- | a square matrix. |
| λ | :- | a symbol or algebraic expression. |
char_poly finds the characteristic polynomial of 𝒜.
This is the determinant of λℐ−𝒜.
char_poly(𝒜,x) = x3 − 15 ∗x2 − 18 ∗x
char_matrix.
cholesky(𝒜);
| 𝒜 | :- | a positive definite matrix containing numeric entries. |
cholesky computes the cholesky decomposition of 𝒜.
It returns {ℒ,𝒰} where ℒ is a lower matrix, 𝒰 is an upper matrix,
𝒜 = ℒ𝒰, and 𝒰 = ℒT .
lu_decom.
coeff_matrix({lin_eqn1,lin_eqn2, …,lin_eqnn});
27
| lin_eqn1,lin_eqn2, …,lin_eqnn | :- | linear equations. Can be of the form equation = number or just equation which is equivalent to equation = 0. |
coeff_matrix creates the coefficient matrix 𝒞 of the linear equations. It returns
{𝒞,𝒳,ℬ} such that 𝒞𝒳 = ℬ.
coeff_matrix

column_dim(𝒜);
| 𝒜 | :- | a matrix. |
column_dim finds the column dimension of 𝒜.
row_dim finds the row dimension of 𝒜.
column_dim(𝒜) = 3
companion(poly,x);
| poly | :- | a monic univariate polynomial in x. |
| x | :- | the variable. |
companion creates the companion matrix 𝒞 of poly.
This is the square matrix of dimension n, where n is the degree of poly w.r.t. x. The entries of 𝒞 are: 𝒞(i,n) = −coeffn(poly,x,i − 1) for i = 1,…,n, 𝒞(i,i − 1) = 1 for i = 2,…,n and the rest are 0.

find_companion.
copy_into(𝒜,ℬ,r,c);
| 𝒜,ℬ | :- | matrices. |
| r,c | :- | positive integers. |
copy_into copies matrix 𝒜 into ℬ with 𝒜(1,1) at ℬ(r,c).
copy_into
augment_columns, extend, matrix_augment, matrix_stack,
stack_rows, sub_matrix.
diagonal({mat1,mat2,
…,matn});28
| mat1,mat2, …,matn | :- | each can be either a scalar expr or a square matrix. |
diagonal creates a matrix that contains the input on the diagonal.
jordan_block.
extend(𝒜,r,c,expr);
| 𝒜 | :- | a matrix. |
| r, c | :- | positive integers. |
| expr | :- | algebraic expression or symbol. |
extend returns a copy of 𝒜 that has been extended by r rows and c columns. The
new entries are made equal to expr.

copy_into, matrix_augment, matrix_stack, remove_columns,
remove_rows.
find_companion(𝒜,x);
| 𝒜 | :- | a matrix. |
| x | :- | the variable. |
Given a companion matrix, find_companion finds the polynomial from which
it was made.
find_companion
companion.
get_columns(𝒜,column_list);
| 𝒜 | :- | a matrix. |
| c | :- | either a positive integer or a list of positive integers. |
get_columns removes the columns of 𝒜 specified in column_list and
returns them as a list of column matrices.
get_rows performs the same task on the rows of 𝒜.
get_columns
get_rows
augment_columns, stack_rows, sub_matrix.
gram_schmidt({vec1,vec2, …,vecn});
29
| vec1,vec2, …,vecn | :- | linearly-independent vectors. Each vector must be written as a list, eg:{1,0,0}. |
gram_schmidt performs the Gram-Schmidt orthonormalisation on the input
vectors. It returns a list of orthogonal normalised vectors.
gram_schmidt({{1,0,0},{1,1,0},{1,1,1}}) = {{1,0,0},{0,1,0},{0,0,1}}
gram_schmidt({{1,2},{3,4}}) = {{
,
},{
,
}}
hermitian_tp(𝒜);
| 𝒜 | :- | a matrix. |
hermitian_tp computes the hermitian transpose of 𝒜.
This is a matrix in which the (i,j)th entry is the conjugate of the (j,i)th entry of 𝒜.
hermitian_tp
tp30 .
hessian(expr,variable_list);
| expr | :- | a scalar expression. |
| variable_list | :- | either a single variable or a list of variables. |
hessian computes the hessian matrix of expr w.r.t. the varibles in
variable_list.
This is an n×n matrix where n is the number of variables and the (i,j)th entry is df(expr,variable_list(i),variable_list(j)).

df31 .
hilbert(square_size,expr);
| square_size | :- | a positive integer. |
| expr | :- | an algebraic expression. |
hilbert computes the square hilbert matrix of dimension square_size.
This is the symmetric matrix in which the (i,j)th entry is 1∕(i + j −expr).

mat_jacobian(expr_list,variable_list);
| expr_list | :- | either a single algebraic expression or a list of algebraic expressions. |
| variable_list | :- | either a single variable or a list of variables. |
mat_jacobian computes the jacobian matrix of expr_list w.r.t.
variable_list.
This is a matrix whose (i,j)th entry is df(expr_list(i),variable_list(j)). The matrix is n ×m where n is the number of variables and m the number of expressions.
mat_jacobian
hessian, df32 .
NOTE: The function mat_jacobian used to be called just "jacobian" however us of that name was in conflict with another Reduce package.
jordan_block(expr,square_size);
| expr | :- | an algebraic expression or symbol. |
| square_size | :- | a positive integer. |
jordan_block computes the square jordan block matrix 𝒥 of dimension
square_size.
The entries of 𝒥 are: 𝒥 (i,i) = expr for i = 1,…,n, 𝒥 (i,i + 1) = 1 for i = 1,…,n − 1, and all other entries are 0.
jordan_block(x,5)
diagonal, companion.
lu_decom(𝒜);
| 𝒜 | :- | a matrix containing either numeric entries or imaginary entries with numeric coefficients. |
lu_decom performs LU decomposition on 𝒜, ie: it returns {ℒ,𝒰} where ℒ is a
lower diagonal matrix, 𝒰 an upper diagonal matrix and 𝒜 = ℒ𝒰.
Caution: The algorithm used can swap the rows of 𝒜 during the calculation. This means that ℒ𝒰 does not equal 𝒜 but a row equivalent of it. Due to this, lu_decom returns {ℒ,𝒰,vec}. The call convert(𝒜,vec) will return the matrix that has been decomposed, ie: ℒ𝒰 = convert(𝒜,vec).
lu := lu_decom
| lu := lu_decom(𝒫) = | ![]() | ||
![]() |
cholesky.
make_identity(square_size);
| square_size | :- | a positive integer. |
make_identity creates the identity matrix of dimension square_size.
make_identity(4) 
diagonal.
matrix_augment({mat1,mat2,
…,matn});33
| mat1,mat2, …,matn | :- | matrices. |
matrix_augment sticks the matrices in matrix_list together
horizontally.
matrix_stack sticks the matrices in matrix_list together vertically.
matrix_augment
matrix_stack
augment_columns, stack_rows, sub_matrix.
matrixp(test_input);
| test_input | :- | anything you like. |
matrixp is a boolean function that returns t if the input is a matrix and nil
otherwise.
matrixp(𝒜) = t
matrixp(doodlesackbanana) = nil
squarep, symmetricp.
minor(𝒜,r,c);
| 𝒜 | :- | a matrix. |
| r, c | :- | positive integers. |
minor computes the (r,c)th minor of 𝒜.
This is created by removing the rth row and the cth column from 𝒜.

remove_columns, remove_rows.
mult_columns(𝒜,column_list,expr);
| 𝒜 | :- | a matrix. |
| column_list | :- | a positive integer or a list of positive integers. |
| expr | :- | an algebraic expression. |
mult_columns returns a copy of 𝒜 in which the columns specified in
column_list have been multiplied by expr.
mult_rows performs the same task on the rows of 𝒜.
mult_columns
mult_rows
add_to_columns, add_to_rows.
pivot(𝒜,r,c);
| 𝒜 | :- | a matrix. |
| r,c | :- | positive integers such that 𝒜(r,c)≠0. |
pivot pivots 𝒜 about its (r,c)th entry.
To do this, multiples of the r’th row are added to every other row in the matrix.
This means that the c’th column will be 0 except for the (r,c)’th entry.

rows_pivot.
pseudo_inverse(𝒜);
| 𝒜 | :- | a matrix containing only real numeric entries. |
pseudo_inverse, also known as the Moore-Penrose inverse, computes the
pseudo inverse of 𝒜.
Given the singular value decomposition of 𝒜, i.e: 𝒜 = 𝒰Σ𝒱T , then the pseudo inverse 𝒜† is defined by 𝒜† = 𝒱Σ†𝒰T . For the diagonal matrix Σ, the pseudoinverse Σ† is computed by taking the reciprocal of only the nonzero diagonal elements.
If 𝒜 is square and non-singular, then 𝒜† = 𝒜. In general, however, 𝒜𝒜†𝒜 = 𝒜, and 𝒜†𝒜𝒜† = 𝒜†.
Perhaps more importantly, 𝒜† solves the following least-squares problem: given a rectangular matrix 𝒜 and a vector b, find the x minimizing ∥𝒜x −b∥2, and which, in addition, has minimum ℓ2 (euclidean) Norm, ∥x∥2. This x is 𝒜†b.

svd.
random_matrix(r,c,limit);
| r,c, limit | :- | positive integers. |
random_matrix creates an r ×c matrix with random entries in the range
−limit < entry < limit.
| imaginary | :- | if on, then matrix entries are x + iy where −limit < x,y < limit. |
| not_negative | :- | if on then 0 < entry < limit. In the imaginary case we have 0 < x,y < limit. |
| only_integer | :- | if on then each entry is an integer. In the imaginary case x,y are integers. |
| symmetric | :- | if on then the matrix is symmetric. |
| upper_matrix | :- | if on then the matrix is upper triangular. |
| lower_matrix | :- | if on then the matrix is lower triangular. |
random_matrix
on only_integer, not_negative, upper_matrix, imaginary;
random_matrix
remove_columns(𝒜,column_list);
| 𝒜 | :- | a matrix. |
| column_list | :- | either a positive integer or a list of positive integers. |
remove_columns removes the columns specified in column_list from
𝒜.
remove_rows performs the same task on the rows of 𝒜.
remove_columns
remove_rows
minor.
See: column_dim.
rows_pivot(𝒜,r,c,{row_list});
| 𝒜 | :- | a matrix. |
| r,c | :- | positive integers such that 𝒜(r,c) neq 0. |
| row_list | :- | positive integer or a list of positive integers. |
rows_pivot performs the same task as pivot but applies the pivot only to the
rows specified in row_list.
rows_pivot
pivot.
simplex(max/min,objective function,{linear
inequalities},[{bounds}]);
| max/min | :- | either max or min (signifying maximise and minimise). |
| objective function | :- | the function you are maximising or minimising. |
| linear inequalities | :- | the constraint inequalities. Each one must be of the form sum of variables (<=,=,>=) number. |
| bounds | :- | bounds on the variables as specified for the LP file format. Each bound is of one of the forms l ≤ v, v ≤ u, or l ≤ v ≤ u, where v is a variable and l, u are numbers or infinity or -infinity |
simplex applies the revised simplex algorithm to find the optimal(either
maximum or minimum) value of the objective function under the linear inequality
constraints.
It returns {optimal value,{ values of variables at this optimal}}.
The {bounds} argument is optional and admissible only when the switch fastsimplex is on, which is the default.
Without a {bounds} argument, the algorithm implies that all the variables are non-negative.
simplex(max,x+y,{x>=10,y>=20,x+y<=25});
***** Error in simplex: Problem has no feasible solution.
simplex(max,10x+5y+5.5z,{5x+3z<=200,x+0.1y+0.5z<=12,
0.1x+0.2y+0.3z<=9, 30x+10y+50z<=1500});
{525.0,{x=40.0,y=25.0,z=0}}
squarep(𝒜);
| 𝒜 | :- | a matrix. |
squarep is a boolean function that returns t if the matrix is square and nil
otherwise.
matrixp, symmetricp.
See: augment_columns.
sub_matrix(𝒜,row_list,column_list);
| 𝒜 | :- | a matrix. |
| row_list, column_list | :- | either a positive integer or a list of positive integers. |
sub_matrix produces the matrix consisting of the intersection of the rows
specified in row_list and the columns specified in column_list.
sub_matrix
augment_columns, stack_rows.
svd(𝒜);
| 𝒜 | :- | a matrix containing only real numeric entries. |
svd computes the singular value decomposition of 𝒜. If A is an m ×n real
matrix of (column) rank r, svd returns the 3-element list {𝒰,Σ,𝒱} where
𝒜 = 𝒰Σ𝒱T .
Let k = min(m,n). Then U is m ×k, V is n ×k, and and Σ = diag(σ1,…,σk), where σi ≥ 0 are the singular values of 𝒜; only r of these are non-zero. The singular values are the non-negative square roots of the eigenvalues of 𝒜T 𝒜.
𝒰 and 𝒱 are such that 𝒰𝒰T = 𝒱𝒱T = 𝒱T 𝒱 = ℐk.
Note: there are a number of different definitions of SVD in the literature, in some of which Σ is square and U and V rectangular, as here, but in others U and V are square, and Σ is rectangular.

swap_columns(𝒜,c1,c2);
| 𝒜 | :- | a matrix. |
| c1,c1 | :- | positive integers. |
swap_columns swaps column c1 of 𝒜 with column c2.
swap_rows performs the same task on 2 rows of 𝒜.
swap_columns
swap_entries.
swap_entries(𝒜,{r1,c1},{r2,c2});
| 𝒜 | :- | a matrix. |
| r1,c1,r2,c2 | :- | positive integers. |
swap_entries swaps 𝒜(r1,c1) with 𝒜(r2,c2).
swap_entries
swap_columns, swap_rows.
symmetricp(𝒜);
| 𝒜 | :- | a matrix. |
symmetricp is a boolean function that returns t if the matrix is symmetric and
nil otherwise.
matrixp, squarep.
toeplitz({expr1,expr2, …,exprn});
34
| expr1,expr2, …,exprn | :- | algebraic expressions. |
toeplitz creates the toeplitz matrix from the expression list.
This is a square symmetric matrix in which the first expression is placed on the diagonal and the i’th expression is placed on the (i-1)’th sub and super diagonals.
It has dimension n where n is the number of expressions.

triang_adjoint(𝒜);
| 𝒜 | :- | a matrix. |
triang_adjoint computes the triangularizing adjoint ℱ of matrix 𝒜 due to the algorithm of Arne Storjohann. ℱ is lower triangular matrix and the resulting matrix 𝒯 of ℱ∗𝒜 = 𝒯 is upper triangular with the property that the i-th entry in the diagonal of 𝒯 is the determinant of the principal i-th submatrix of the matrix 𝒜.
triang_adjoint

vandermonde({expr1,expr2, …,exprn});
35
| expr1,expr2, …,exprn | :- | algebraic expressions. |
Vandermonde creates the Vandermonde matrix from the expression list. This is
the square matrix in which the (i,j)th entry is expri(j−1). It has dimension n,
where n is the number of expressions.

kronecker_product(M1,M2)
| M1,M2 | :- | Matrices |
kronecker_product creates a matrix containing the Kronecker product (also
called direct product or tensor product) of its arguments.
a1 := mat((1,2),(3,4),(5,6))$
a2 := mat((1,1,1),(2,z,2),(3,3,3))$
kronecker_product(a1,a2);

By turning the fast_la switch on, the speed of the following functions will be increased:
| add_columns | add_rows | augment_columns | column_dim |
| copy_into | make_identity | matrix_augment | matrix_stack |
| minor | mult_column | mult_row | pivot |
| remove_columns | remove_rows | rows_pivot | squarep |
| stack_rows | sub_matrix | swap_columns | swap_entries |
| swap_rows | symmetricp |
The increase in speed will be insignificant unless you are making a significant number(i.e: thousands) of calls. When using this switch, error checking is minimised. This means that illegal input may give strange error messages. Beware.
Many of the ideas for this package came from the Maple[3] Linalg package [4].
The algorithms for cholesky, lu_decom, and svd are taken from the book Linear Algebra - J.H. Wilkinson & C. Reinsch[5].
The gram_schmidt code comes from Karin Gatermann’s Symmetry package[6] for REDUCE.
| Up | Next | Prev | PrevTail | Front |