numeric::matlinsolve

Solve a linear matrix equation

Use only in the MuPAD Notebook Interface.

This functionality does not run in MATLAB.

Syntax

numeric::matlinsolve(A, B, options)

Description

numeric::matlinsolve(A, B) returns the matrix solution X of the matrix equation AX = B together with the kernel of the matrix A.

numeric::matlinsolve is a fast numerical linear solver for both sparse and dense systems. It is also a recommended solver for linear systems with exact or symbolic coefficients (use option Symbolic).

If no return type is specified via the option ReturnType = d, the domain type of the return data depends on the type of the input matrix A:

  • The special solution X as well as the kernel of an array A are returned as arrays.

  • The special solution X as well as the kernel of an hfarray A are returned as hfarrays.

  • For a dense matrixA of type Dom::DenseMatrix(), both the special solution X as well as the kernel of A are returned as matrices of type Dom::DenseMatrix() over the ring of MuPAD® expressions.

  • For all other matrices of category Cat::Matrix, both the special solution X as well as the kernel of A are returned as matrices of type Dom::Matrix() over the ring of MuPAD expressions. This includes input matrices A of type Dom::Matrix(...), Dom::SquareMatrix(...), Dom::MatrixGroup(...) etc.

Without Symbolic, exact numerical input data such as PI + sqrt(2), sin(3) etc. are converted to floating-point numbers. Floating point arithmetic is used. Its precision is given the environment variable DIGITS. If symbolic data are found that cannot be converted to floating-point numbers, numeric::matlinsolve automatically switches to its symbolic mode, issuing a warning. This warning may be suppressed via NoWarning.

With Symbolic, symbolic parameters in the coefficient matrix A as well as in the right hand side B are accepted and processed without a warning.

With SofwareFloats, the right hand side B may contain symbolic parameters that cannot be converted to floating-point numbers. All entries of the coefficient matrix A, however, must be convertible to floating-point numbers.

With HardwareFloats, neither A nor B must contain symbolic parameters that cannot be converted to floating-point numbers.

X is a special solution of the equation AX = B. If A has a non-trivial kernel, the solution X is not unique.

    Note:   The result computed with HardwareFloats may differ from the solution computed with SoftwareFloats or Symbolic! In particular, this is the case for systems with a non-trivial kernel.

Cf. Example 9.

The n×d matrix KernelBasis is the most general solution of AX = 0. Its columns span the d-dimensional kernel of A.

Thus, the kernel of A may be computed via numeric::matlinsolve(A, [0, ..., 0])[2].

If the kernel is 0-dimensional, the return value of KernelBasis is the integer 0. If KernelBasis is returned as an array, the dimension d of the kernel is d = op(KernelBasis, [0, 3, 2]]). If KernelBasis is returned as a matrix of type Dom::Matrix() or Dom::DenseMatrix(), the dimension d of the kernel is d = KernelBasis::dom::matdim(KernelBasis)[2].

    Note:   Due to round-off errors, some or all basis vectors in the kernel of A may be missed in the numerical modes.

The special solution X in conjunction with KernelBasis provides the general solution of AX = B. Solutions generated without the option ShowAssumptions are valid for arbitrary complex values of the symbolic parameters which may be present in A and B. If no such solution exists, then [FAIL,NIL] is returned. Solutions that are valid only for special values of the symbolic parameters may be obtained with ShowAssumptions. See Example 3, Example 4, and Example 7.

numeric::matlinsolve internally uses a sparse representation of the matrices. It is suitable for solving large sparse systems. See Example 5.

    Note:   numeric::matlinsolve does not react to any assumptions on symbolic parameters in A,B that are set via assume.

    Note:   Gaussian elimination with partial pivoting is used. Without the option Symbolic, the pivoting strategy takes care of numerical stabilization. With Symbolic, exact data are assumed. The symbolic pivoting strategy tries do maximize speed and does not take care of numerical stabilization! Do not use Symbolic for linear systems involving floating-point entries! See Example 6.

    Note:   Apart from matrices of type Dom::Matrix(...), Cat::Matrix objects A from matrix domains such as Dom::DenseMatrix(...) or Dom::SquareMatrix(...) are internally converted to arrays over expressions via expr(A). Note that the option Symbolic should be used if the entries cannot be converted to numerical expressions.

    Note that linalg::matlinsolve must be used, when the solution is to be computed over the component domain. See . Example 8.

Environment Interactions

Without the option Symbolic, the function is sensitive to the environment variable DIGITS, which determines the numerical working precision.

Examples

Example 1

We use equivalent input formats (B1, B2) to represent a vector with components [a, π]. First, this vector is defined as a 2-dimensional array:

A := array(1..2, 1..3, [[1, 2, 3],[1, 1, 2]]):
B1 := array(1..2, 1..1, [[a], [PI]]):
numeric::matlinsolve(A, B1)

Next, we use a 1-dimensional array and compute an exact solution:

B2 := array(1..2, [a, PI]):
numeric::matlinsolve(A, B2, Symbolic)

Now, a list is used to specify the vector. No internal assumptions were used by numeric::matlinsolve to obtain the solution:

B3 := [a, PI]:
numeric::matlinsolve(A, B3, ShowAssumptions)

Finally, we use Dom::Matrix objects to specify the system. Note that the results are returned as corresponding matrix objects:

A := matrix([[1, 2, 3],[1, 1, 2]]):
B4 := matrix([a, PI]):
numeric::matlinsolve(A, B4)

delete A, B1, B2, B3, B4:

Example 2

We invert a matrix by solving AX = 1:

A := hfarray(1..3, 1..3, [[1, 1, 0], [0, 1, 1], [0, 0, 1]]):
B := matrix::identity(3, 3):
InverseOfA := numeric::matlinsolve(A, B, Symbolic)[1]

delete A, B, InverseOfA:

Example 3

We solve an equation with a symbolic parameter x:

A := matrix([[2, 2, 3], [1, 1, 2], [3, 3, 5]]):
B := matrix([sin(x)^2, cos(x)^2, 0]):
[X, Kernel, Constraints, Pivots] :=
 numeric::matlinsolve(A, B, Symbolic, ShowAssumptions)

This solution holds subject to the constraint sin(x)2 + cos(x)2 = 0 on the parameter x. numeric::matlinsolve does not investigate the Constraints and does not realize that they cannot be satisfied. We check the consistency of the "result" by inserting the solution into the original system. Since the input matrix A was of type Dom::Matrix(), the results X and Kernel were returned as corresponding matrices. The overloaded operators * and - for matrix multiplication and subtraction can be used:

A*X - B, A*Kernel

delete A, B, X, Kernel, Constraints, Pivots:

Example 4

We give a further demonstration of the option ShowAssumptions. The following system does not have a solution for all values of the parameter a:

A := array(1..2, 1..2, [[1, 1], [1, 1]]):
B := array(1..2, 1..1, [[1], [a]]):
numeric::matlinsolve(A, B, Symbolic)

With ShowAssumptions, numeric::matlinsolve investigates under which conditions (on the parameter a a) there is a solution:

numeric::matlinsolve(A, B, Symbolic, ShowAssumptions)

We conclude that there is a 1-dimensional solution space for a = 1.

delete A, B:

Example 5

We solve a sparse system with 3 diagonal bands:

n := 100: 
A := matrix(n, n, [1, -2, 1], Banded):
B := array(1..n, [1 $ n]): 
numeric::matlinsolve(A, B)

[matrix([[-50.0], [-99.0], [-147.0], [-194.0], [-240.0], [-285.0], [-329.0], [-372.0], [-414.0], [-455.0], [-495.0], [-534.0], [-572.0], [-609.0], [-645.0], [-680.0], [-714.0], [-747.0], [-779.0], [-810.0], [-840.0], [-869.0], [-897.0], [-924.0], [-950.0], [-975.0], [-999.0], [-1022.0], [-1044.0], [-1065.0], [-1085.0], [-1104.0], [-1122.0], [-1139.0], [-1155.0], [-1170.0], [-1184.0], [-1197.0], [-1209.0], [-1220.0], [-1230.0], [-1239.0], [-1247.0], [-1254.0], [-1260.0], [-1265.0], [-1269.0], [-1272.0], [-1274.0], [-1275.0], [-1275.0], [-1274.0], [-1272.0], [-1269.0], [-1265.0], [-1260.0], [-1254.0], [-1247.0], [-1239.0], [-1230.0], [-1220.0], [-1209.0], [-1197.0], [-1184.0], [-1170.0], [-1155.0], [-1139.0], [-1122.0], [-1104.0], [-1085.0], [-1065.0], [-1044.0], [-1022.0], [-999.0], [-975.0], [-950.0], [-924.0], [-897.0], [-869.0], [-840.0], [-810.0], [-779.0], [-747.0], [-714.0], [-680.0], [-645.0], [-609.0], [-572.0], [-534.0], [-495.0], [-455.0], [-414.0], [-372.0], [-329.0], [-285.0], [-240.0], [-194.0], [-147.0], [-99.0], [-50.0]]), 0]
delete n, A, B:

Example 6

The option Symbolic should not be used for equations with floating-point coefficients, because the symbolic pivoting strategy favors efficiency instead of numerical stability.

A := array(1..2, 1..2, [[1, 10^20], [1, 1]]):
B := array(1..2, 1..1, [[10^20], [0]]):

The float approximation of the exact solution is:

map(numeric::matlinsolve(A, B, Symbolic)[1], float)

We now convert the exact input data to floating-point approximations:

A := map(A, float): B := map(B, float):

The default pivoting strategy of the floating-point algorithm stabilizes floating-point operations. Consequently, one gets a correct result:

numeric::matlinsolve(A, B)[1]

With the option Symbolic, however, the pivoting strategy optimizes speed, assuming exact arithmetic. Numerical instabilities may occur if floating-point coefficients are involved. The following result is caused by internal round-off effects ("cancellation"):

numeric::matlinsolve(A, B, Symbolic)[1]

We need to increase DIGITS to obtain a better result:

DIGITS := 20:
numeric::matlinsolve(A, B, Symbolic)[1]

delete A, B, DIGITS:

Example 7

We demonstrate how a complete solution of the following linear system with symbolic parameters may be found:

A := array(1..3, 1..2, [[1, 1], [a, b], [1, c]]):
B := array(1..3, 1..1, [[1], [1], [1]]):
numeric::matlinsolve(A, B, Symbolic, ShowAssumptions)

This is the general solution assuming ab. We now set b = a to investigate further solution branches:

A := subs(A, b = a):
numeric::matlinsolve(A, B, Symbolic, ShowAssumptions)

This is the general solution for a = b, assuming c ≠ 1. We finally set c = 1 to obtain the last solution branch:

A := subs(A, c = 1):
numeric::matlinsolve(A, B, Symbolic, ShowAssumptions)

From the constraint on a, we conclude that there is a 1-dimensional solution space for the special values a = b = c = 1 of the symbolic parameters.

delete A, B:

Example 8

Matrices from a domain such as Dom::Matrix(...) are converted to arrays with numbers or expressions. Hence, numeric::matlinsolve finds no solution for the following system:

M := Dom::Matrix(Dom::IntegerMod(7)):
A := M([[1, 4], [6, 3], [3, 2]]): 
B := M([[9], [5], [0]]):
numeric::matlinsolve(A, B)

Use linalg::matlinsolve to solve the system over the coefficient field of the matrices. A solution does exist over the field Dom::IntegerMod(7):

linalg::matlinsolve(A, B)

delete M, A, B:

Example 9

We demonstrate the difference between Symbolic, HardwareFloats, and SoftwareFloats. The following matrix A has a 1-dimensional kernel. Due to round-off, a further spurious kernel vector appears with SoftwareFloats. No kernel vector is detected with HardwareFloats:

A := matrix([[2*10^14 + 2, 2*10^(-9), 2*10^(-4)], 
             [3*10^15 + 3, 3*10^(-8), 3*10^(-3)],
             [4*10^16 + 4, 4*10^(-7), 4*10^(-2)]
            ]): 
b := matrix([2*10^(-9), 3*10^(-8), 4*10^(-7)]):
float(numeric::matlinsolve(A, b, Symbolic));
numeric::matlinsolve(A, b, SoftwareFloats);
numeric::matlinsolve(A, b, HardwareFloats)

delete A, b:

Parameters

A

An m×n matrix of domain type DOM_ARRAY, DOM_HFARRAY, or of category Cat::Matrix

B

An m×p matrix of domain type DOM_ARRAY, DOM_HFARRAY, or of category Cat::Matrix. Column vectors B may also be represented by a 1-dimensional array(1..m, [B1, B2, …] )), a 1-dimensional hfarray(1..m, [B1, B2, …] )), or by a list [B1, B2, …].

Options

Hard, HardwareFloats, Soft, SoftwareFloats

With Hard (or HardwareFloats), computations are done using fast hardware float arithmetic from within a MuPAD session. Hard and HardwareFloats are equivalent. With this option, the input data are converted to hardware floats and processed by compiled C code. The result is reconverted to MuPAD floats and returned to the MuPAD session.

With Soft (or SoftwareFloats) computations are dome using software float arithmetic provided by the MuPAD kernel. Soft and SoftwareFloats are equivalent. SoftwareFloats is used by default if the current value of DIGITS is larger than 15 and the input matrix A is not of domain type DOM_HFARRAY.

Compared to the SoftwareFloats used by the MuPAD kernel, the computation with HardwareFloats may be many times faster. Note, however, that the precision of hardware arithmetic is limited to about 15 digits. Further, the size of floating-point numbers may not be larger than approximately 10308 and not smaller than approximately 10- 308.

If no HardwareFloats or SoftwareFloats are requested explicitly, the following strategy is used: If the current value of DIGITS is smaller than 16 or if the matrix A is a hardware float array of domain type DOM_HFARRAY, then hardware arithmetic is tried. If this is successful, the result is returned.

If the result cannot be computed with hardware floats, software arithmetic by the MuPAD kernel is tried.

If the current value of DIGITS is larger than 15 and the input matrix A is not of domain type DOM_HFARRAY, or if one of the options Soft, SoftwareFloats or Symbolic is specified, MuPAD computes the result with its software arithmetic without trying to use hardware floats first.

There may be several reasons for hardware arithmetic to fail:

  • The current value of DIGITS is larger than 15.

  • The data contains symbolic objects.

  • The data contains numbers larger than 10308 or smaller than 10- 308 that cannot be represented by hardware floats.

If neither HardwareFloats nor SoftwareFloats is specified, the user is not informed whether hardware floats or software floats are used.

If HardwareFloats are specified but fail due to one of the reasons above, a warning is issued that the (much slower) software floating-point arithmetic of the MuPAD kernel is used.

Note that HardwareFloats can only be used if all input data can be converted to floating-point numbers.

The trailing digits in floating-point results computed with HardwareFloats and SoftwareFloats may differ.

    Note:   For ill-conditioned matrices, the results returned with HardwareFloats and SoftwareFloats may differ significantly! See Example 9.

Symbolic

Prevents the conversion of the input data to floating-point numbers. Exact arithmetic is used. This option overrides HardwareFloats and SoftwareFloats.

This option must be used if the matrix A contains symbolic parameters that cannot be converted to floating-point numbers.

    Note:   This option should not be used for matrices with floating-point entries! Numerical instabilities may occur in floating-point operations. Cf. Example 6.

ShowAssumptions

Returns information on internal assumptions on symbolic parameters in A and B. With this option, either exact arithmetic or SoftwareFloats are used.

This option is only useful if the matrices contain symbolic parameters. Consequently, it should only be used in conjunction with the option Symbolic.

    Note:   This option changes the format of the return value to [X, KernelBasis, Constraints, Pivots].

X and KernelBasis represent the general solution subject to Constraints and Pivots.

Constraints is a list of equations for symbolic parameters in B which are necessary and sufficient for AX = B to be solvable.

Such constraints arise if Gaussian elimination leads to equations of the form 0 = c, where c is some expression involving symbolic parameters contained in B. All such equations are collected in Constraints; numeric::matlinsolve assumes that these equations are satisfied and returns a special solution X.

If no such constraints arise, the return value of Constraints is the empty list.

Pivots is a list of inequalities involving symbolic parameters in A. Internally, division by pivot elements occurs in the Gaussian elimination. The expressions collected in Pivots are the numerators of those pivot elements that involve symbolic parameters contained in A. If only numerical pivot elements are used, then the return value of Pivots is the empty list.

    Note:   Constraints usually is a list of non-linear equations for the symbolic parameters. It is not investigated by numeric::matlinsolve, i.e., solutions may be returned, even if the Constraints cannot be satisfied. See Example 3.

    Note:   This option changes the return strategy for "unsolvable" systems. Without the option ShowAssumptions, the result [FAIL,NIL] is returned, whenever Gaussian elimination produces an equation 0 = c with non-zero c. With ShowAssumptions, such equations are returned via Constraints, provided c involves symbolic parameters. If c is a purely numerical value, then [FAIL, NIL, [], []] is returned.

See Example 3, Example 4, and Example 7.

NoWarning

Suppresses warnings

If symbolic coefficients are found, numeric::matlinsolve automatically switches to the Symbolic mode with a warning. With this option, this warning is suppressed; numeric::matlinsolve still uses the symbolic mode for symbolic coefficients, i.e., exact arithmetic without floating-point conversions is used.

ReturnType

Option, specified as ReturnType = d

Return the (special) solution and the kernel as matrices of domain type d. The following return types d are available: DOM_ARRAY, DOM_HFARRAY, Dom::Matrix(), or Dom::DenseMatrix().

Sparse

Use a sparse internal representation for matrices.

This option only has an effect when used in conjunction with HardwareFloats. With the Sparse option, the linear solver uses a sparse representation of the matrices to save memory and increase efficiency. However, if the coefficient matrix is not sparse, this option will cost some additional memory and runtime.

Return Values

Without the option ShowAssumptions, a list [X, KernelBasis] is returned. The (special) solution X is an n×p matrix. KernelBasis is an n×d matrix (d is the dimension of the kernel of A). Its columns span the kernel of A. If the kernel is trivial, KernelBasis is the integer 0.

[FAIL, NIL] is returned if the system is not solvable.

With ShowAssumptions, a list [X, KernelBasis, Constraints, Pivots] is returned. The lists Constraints and Pivots contain equations and inequalities involving symbolic parameters in A and B. Internally these were assumed to hold true when solving the system. [FAIL, NIL, [], []] is returned if the system is not solvable.

Was this topic helpful?