Matrices

Creating matrices

Matrices may be created as follows:

B=(2, 3, 4 ; 4, 6, 6)
2  3  4
4  6  6

If the rest of values in a row are zero, you may use a `;` and begin entering the next row. Zeros will be added:

C=(1 ; 0, 1 ; 0, 0, 3)
1  0  0
0  1  0
0  0  3

The following functions also create matrices:

• `matrix(r,c,v)` returns a new matrix with `r` rows, `c` columns, and all entries are equal to value `v`.
• `matrixd(n,v)` returns a new diagonal matrix of size `n` and all diagonal entries have a value `v`.
• `matrixd(n)`. If `n` is an integer, it returns a new identity matrix of size `n`.
• `matrixd(a)`. If `a` is an array, it returns a diagonal matrix such that the diagonal terms are the entries of array `a`.

Matrices can also be created by creating any entry:

D[2,3]=-2
-2
D
0  0  0
0  0 -2

You may use functions `joinright()` and `joinbottom()` to augment matrices. A typical session follows.

A = ( 1, 0 ; 3, 2 )
1  0
3  2
B = ( 5 ; 7 )
5
7
C = joinright( A, B )
1  0  5
3  2  7
D = ( 9 8 2 )
9  8  2
E = joinbottom( C, D )
1  0  5
3  2  7
9  8  2

Use function `submatrix()` to build a smaller matrix from an existing one. Functions `row()` and `column()` are special cases of function `submatrix()`.

Also note that column matrices and row matrices are treated like arrays too.

Basic functionality

Use the `+` and `-` operators provided the matrices have the same dimensions.

Operator `*` (matrix multiplication) works as long as the inner dimensions of the matrices are the same.

The following functions can be used with matrix arguments:

• `b=inverse(a)` returns the inverse of matrix `a`.
• `b=pseudoinv(a)` returns the pseudo inverse of matrix `a`.
• `r=determinant(a)` returns the determinant of matrix `a`.
• `r=trace(a)` returns the trace of matrix `a` (sum of diagonal terms).
• `transpose(v)` returns the transpose of matrix `v`.
• `x = solve(a,b)` solves the linear set of equations A*x=B.
• `apply(a, f)` transforms matrix `a` by applying function `f` to each of the entries.

Solving sets of linear equations

Use function `solve()` to solve a set of linear equations of the form:

A*x = B

where both A and B are matrices with equal number of rows.
There are three cases to consider. Assuming matrix A has `r` rows and `c` columns, we show each case with an example.

Case r = c

Given `A` and `B`:

A=(1,2,3;5,6,8;0,-2,3)
1  2  3
5  6  8
0 -2  3
B=(3;5;7)
3
5
7

We solve this systems as follows:

X = solve(A,B)
-14/13
-19/26
24/13

Calcugator uses an LU Decomposition to solve the system of equations. Matrices `A` and `B` can be real or complex.

Case r<c

Let's create a matrix with 3 rows and 5 columns with random integers between 0 and 10:

A=randommatrix(3,5,10)
6  5  2  8  8
2  1  7  1  7
6  2  8  8  5
B=(7;5;6)

The systems is solved using the same function as above:

X = solve(A,B)
0.1559
0.1811
0.1572
0.1387
0.4667

To solve the system Calcugator uses the pseudo inverse of matrix `A` computed with a Singular Value Decomposition. Both matrices `A` and `B` must be real.

Case r>c

Given `A` and `B`:

A=(3,1 ; 6,6 ; 3,3 ; 0, 4)
3  1
6  6
3  3
0  4
B=(10 ; 9 ; 0 ; 6)
10
9
0
6

Again, use function `solve`:

X = solve(A,B)
0.9425
0.6896

In this case, Calcugator finds the least squares solution. This solution minimizes the norm of A*X-B. Both matrices `A` and `B` must be real.

Matrix decomposition

Three types of matrix decompositions are supported.

QR Decomposition

A matrix A (square or rectangular matrix with more rows than columns) that has full rank, can be decomposed into the product:

A = Q*R

where Q is an orthonormal matrix and R is upper triangular. Use function `qrgetQ` or function `qrgetR` to compute those matrices.

Example:

A=(3,1 ; 6,6 ; 3,3 ; 0, 4)
3  1
6  6
3  3
0  4
Q=qrgetQ(A)
0.4082  -0.3790
0.8164   0.1516
0.4082   0.0758
0.0      0.9097
R=qrgetR(A)
7.3484  6.5319
0.0     4.3969
Q*R
3.0000  1.0000
6.0000  6.0000
3.0000  3.0000
0.0     4.0000

For properties of matrices `Q` and `R`, consult a book on linear Algebra.

Singular Value Decomposition

A square matrix A (or rectangular with more rows than columns), can be decomposed into the product:

A = U*S*VT

where U and V are orthonormal and S is diagonal. Use functions `svdgetU()`, `svdgetS()` and `svdgetV()` to compute those matrices.

Example:

A=(3,1 ; 6,6 ; 3,3 ; 0, 4)
3  1
6  6
3  3
0  4
U=qrgetU(A)
0.2685  -0.4880
0.8228  -0.1124
0.4114  -0.0562
0.2856   0.8636
S=qrgetS(A)
10.3037  0.0
0.0     3.1358
V=qrgetV(A)
0.6771  -0.7358
0.7358   0.6771
U*S*transpose(V)
3.0000      1.0000
5.9999      5.9999
2.9999      2.9999
1.1102E-15  4.0000
normf(ans-A)
2.4825E-15

The norm returned by function `normf()` returns the Frobenius norm of the matrix.

Eigenvalue Decomposition

A square matrix A can be decomposed into the following:

A = S*L*S-1

where S is a matrix such that its columns are the eigenvectors and L is a diagonal matrix containing the eigenvalues. Use functions `eigenmatrix()` and `eigenvalues()` to compute those matrices.

Example:

A=(2,4,6 ; 4,5 ; 0, -3, 5)
2   4  6
4   5  0
0  -3  5
S=eigenmatrix(A)
-0.8436   0.0885 + 0.8650 i    0.0885 - 0.8650 i
0.4918   0.9228 + 0.5488 i    0.9228 - 0.5488 i
0.2150  -0.9174 + 0.3833 i   -0.9174 - 0.3833 i
L=eigenvalues(A)
-1.8613  0                 0
0       6.9306+2.6011 i   0
0       0                 6.9306-2.6011 i
S*L*inverse(S)
2.0000+1.1102E-16 i   4.0000   5.9999+1.102E-15 i
3.9999                4.9999  -2.1094E-15+4.4408E-16 i
-3.0531E-16           -3.0000   5.0000+4.4408E-16 i

Notice all imaginary values are small; let's take the real part of the matrix:

Re(ans)
2.0000        4.0000   5.9999
3.9999        4.9999  -2.1094E-15
-3.0531E-16  -3.0000   5.0000

LU Decomposition

Any matrix A can be decomposed as the multiplication of a lower triangular matrix L times and upper triangular matrix U. If row permutations were required, then a permutation matrix P needs to premultiply the product L*U in order to retrieve A. To obtain matrices L, U and P use functions `lugetL()`, `lugetU()` and `lugetP()` respectively.

A = ( 0, 2, 3, 0 ; 1, 1, 0, 5 ; 1, 2, 3, 1 )
0  2  3  0
1  1  0  5
1  2  3  1
L = lugetL( A )
1   0   0
0   1   0
1  1/2  1
U = lugetU( A )
1  1   0    5
0  2   3    0
0  0  3/2  -4
P = lugetP( A )
0  1  0
1  0  0
0  0  1
P*L*U
0  2  3  0
1  1  0  5
1  2  3  1

Echelon forms and Null Space

To find the echelon form of a matrix, use function `echelon()`.

A = ( 0, 2, 3, 0 ; 1, 1, 0, 5 ; 1, 2, 3, 1 )
0  2  3  0
1  1  0  5
1  2  3  1
echelon( A )
1  1   0    5
0  2   3    0
0  0  3/2  -4

To find the reduced echelon form of a matrix, use function `echelonr()`

echelonr( A )
1  0  0   1
0  1  0   4
0  0  1  -8/3

Function `nullspace` returns a matrix that spans the null space of a matrix.

A = ( 1, 2, 2, 3 ; 8, -8, -2, 6 ; -2, 4, 2, 0 )
1   2   2   3
8  -8  -2   6
-2   4   2   0
N = nullspace( A )
-1/2  -3/2
-3/4  -3/4
1     0
0     1

Function `solven()` solves exactly for a matrix equation A*X=B for the case that the number of columns of matrix A is bigger than the number of rows. Function `solven()` returns a solution such that adding a linear combination of columns of the null space, the equality still hods.

b = ( 1 ; -1 ; 1 )
1
-1
1
X = solven( A, b )
1/4
3/8
0
0
X2 = X + 2*column(N,1) - 7*column(N,2)
39/4
33/8
2
-7
A*X2
1
-1
1

Linear algebra

The following functions are implemented:

• `m=minor(a,i,j)` returns the minor of matrix `a` corresponding to the entry at row `i` and column `j`.
• `m=cofactor(a,i,j)` returns the cofactor of matrix `a` corresponding to the entry at row `i` and column `j`. The cofactor is the minor times `-1` to the power `i+j`.
• `m=minormatrix(a,i,j)` returns the matrix obtained by deleting row `i` and column `j` of matrix `a`. Note that `minor(a,i,j)` is equal to the determinant of `minormatrix(a,i,j)`.
• `L=eigenvalues(a)` returns a diagonal matrix with all eigenvalues of matrix `a`. Matrix `a` must be a real matrix.
• `V=eigenmatrix(a)` returns a square matrix with all eigenvectors of matrix `a`. The eigenvectors are the columns of the returned matrix. Matrix `a` must be a real matrix.
• `cond(a)` returns the condition number of matrix `a`.
• `rank(a)` returns the rank of matrix `a`.
• `norm1(a)` returns the norm1 of matrix `a`.
• `norm2(a)` returns the norm2 of matrix `a`.
• `normf(a)` returns the Frobenius norm of matrix `a`.
• `norminf(a)` returns the infinity norm of matrix `a`.
• `mexp(a)` returns a matrix equal to `ea`. Matrix `a` must be square.
• `msqrt(a)` returns the square root of matrix `a`. Matrix `a` must be square.
• `mlog(a)` returns the logarithm of matrix `a`. Matrix `a` must be square.
• `mpow(a, p)` returns the matrix equal to `a` to the power `p`. Matrix `a` must be square.

Special functionality

When testing for singularity while finding an inverse, Calcugator uses a default tolerance equal to `1.e-40`.

To change this value do:

• `set("matrixSingularityTolerance", tol)`

Where `tol` is the new tolerance to be used in the check for singularity.

Example:

set("matrixSingularityTolerance", 1e-20)
Property "matrixSingularityTolerance" was set to 1.0E-20.

sets the tolerance to `1e-20`.

Other functions to use with matrices are:

• `c=row(a,i)` returns a row matrix (a matrix with only one row) equal to the `ith` row of matrix `a`.
• `c=column(a,i)` returns a column matrix (a matrix with only one column) equal to the `ith` column of matrix `a`.
• `m=submatrix(a, r1,c1, r2,c2)` returns a matrix obtained from row `r1` to row `r2` and from column `c1` to column `c2` of matrix `a`.
• `isHermitian(a)` returns `true` if matrix `a` is Hermitian.
• `rowadd(a,i,j)` is an educational function. It adds row `i` to row `j` of matrix `a`.
• `rowtimes(a,i,c)` is an educational function. It multiplies row `i` of matrix `a` times values `c`.

Random matrices

Matrices with random values can be created using the `randommatrix()` function.

• `a=randommatrix( rws, cls )` returns a matrix with `rws` rows and `cls` columns with real random values between `0.` and `1.0`.
• `a=randommatrix( rws, cls, value )` returns a matrix with `rws` rows and `cls` columns with random values. If `value` is real, all values are real between `0.` and `value`. If `value` is integer, all values are integers between `0` and `value` inclusive.
• `a=randommatrix( rws, cls, value1, value2 )` returns a matrix with `rws` rows and `cls` columns with random values. If any of the values is real, all values are real between `value1` and `value2`. If both values are integer, all values are integer between `value1` and `value2` inclusive.

Example:

A random matrix with real values between 0. and 5.:

a=randommatrix(4,4,5.)
4.8761  0.5769  0.6039  4.5528
2.2761  4.1217  1.3126  0.0923
0.6550  4.8586  1.3032  4.7012
3.2179  0.6328  1.2388  2.8118

A random matrix with integer values between 3 and 7:

a=randommatrix(4,4,3,7)
7  6  3  5
6  7  5  6
4  4  4  7
3  7  3  5