linSolNonnegdef¶
Solves a real symmetric nonnegative definite system of linear equations \(Ax = b\). Using options, computes a Cholesky factorization of the matrix A, such that \(A = R^TR = LL^T\). Computes the solution to \(Ax = b\) given the Cholesky factor.
Synopsis¶
linSolNonnegdef (a, b)
Required Arguments¶
- float
a[[]](Input) - Array of size
n×ncontaining the matrix. - float
b[](Input) - Array of size
ncontaining the right-hand side.
Return Value¶
Using required arguments, linSolNonnegdef returns a solution x of the
linear system. If no value can be computed, None is returned.
Optional Arguments¶
factor(Output)- An array of size
n×ncontaining the \(LL^T\) factorization of A. When this option is specified, the space for the factor matrix is allocated bylinSolNonnegdef. The lower-triangular part of the factor array contains L, and the upper-triangular part contains \(L^TR\). inverse(Output)- An array of size
n×ncontaining the inverse of A. The space for this array is allocated bylinSolNonnegdef. tolerance(Input)Tolerance used in determining linear dependence. See the documentation for machine (
machine(float)) in Chapter 12, “Utilities.”Default:
tolerance=100 * machine(4)factorOnly- Compute the \(LL^T\) factorization of A only. The argument
bis ignored, and the optional argumentfactoris required. solveOnly- Solve \(Ax = b\) using the factorization previously computed by this
function. The argument
ais ignored, and the optional argumentfactoris required. inverseOnly- Compute the inverse of A only. The argument
bis ignored, and either the optional argumentinverseis required.
Description¶
The function linSolNonnegdef solves a system of linear algebraic
equations having a symmetric nonnegative definite (positive semidefinite)
coefficient matrix. It first computes a Cholesky (\(LL^T\) or
\(R^TR\)) factorization of the coefficient matrix A.
The factorization algorithm is based on the work of Healy (1968) and proceeds sequentially by columns. The i-th column is declared to be linearly dependent on the first i − 1 columns if
where ε (specified in tolerance) may be set by the user. When a linear
dependence is declared, all elements in the i-th row of R (column of
L) are set to zero.
Modifications due to Farebrother and Berry (1974) and Barrett and Healy
(1978) for checking for matrices that are not nonnegative definite also are
incorporated. The function linSolNonnegdef declares A to not be
nonnegative definite and issues an error message if either of the following
conditions are satisfied:
- \(a_{ii} - \textstyle\sum_{j=1}^{i-1} r_{ji}^2 < -\varepsilon | a_{ii} |\)
- \(r_{ii} = 0 \text{ and } | a_{ik} - \sum_{j=1}^{i-1} r_{ji} r_{jk} | > \varepsilon \sqrt{a_{ii} a_{kk}}, k > i\)
Healy’s (1968) algorithm and the function linSolNonnegdef permit the
matrices A and R to occupy the same storage. Barrett and Healy (1978) in
their remark neglect this fact. The function linSolNonnegdef uses
for \(a_{ii}\) in the above condition 2 to remedy this problem.
If an inverse of the matrix A is required and the matrix is not (numerically) positive definite, then the resulting inverse is a symmetric \(g_2\) inverse of A. For a matrix G to be a \(g_2\) inverse of a matrix A, G must satisfy conditions 1 and 2 for the Moore-Penrose inverse, but generally fail conditions 3 and 4. The four conditions for G to be a Moore-Penrose inverse of A are as follows:
- AGA = A
- GAG = G
- AG is symmetric
- GA is symmetric
The solution of the linear system \(Ax = b\) is computed by solving the factored version of the linear system \(R^TRx = b\) as two successive triangular linear systems. In solving the triangular linear systems, if the elements of a row of R are all zero, the corresponding element of the solution vector is set to zero. For a detailed description of the algorithm, see Section 2 in Sallas and Lionti (1988).
Examples¶
Example 1¶
A solution to a system of four linear equations is obtained. Maindonald (1984, pp. 83-86 and 104-105) discusses the computations for the factorization and solution to this problem.
from numpy import *
from pyimsl.math.linSolNonnegdef import linSolNonnegdef
from pyimsl.math.writeMatrix import writeMatrix
a = array([[36.0, 12.0, 30.0, 6.0],
[12.0, 20.0, 2.0, 10.0],
[30.0, 2.0, 29.0, 1.0],
[6.0, 10.0, 1.0, 14.0]])
b = ([18.0, 22.0, 7.0, 20.0])
# Solve Ax = b for x
x = linSolNonnegdef(a, b)
# Print solution, x, of Ax = b
writeMatrix("Solution, x", x)
Output¶
Solution, x
1 2 3 4
0.167 0.500 0.000 1.000
Example 2¶
The symmetric nonnegative definite matrix in the initial example is used to
compute the factorization only in the first call to linSolNonnegdef. On
the second call, both the \(LL^T\) factorization and the right-hand side
vector in the first example are used as the input to compute a solution x.
It also illustrates another way to obtain the solution array x.
from numpy import *
from pyimsl.math.linSolNonnegdef import linSolNonnegdef
from pyimsl.math.writeMatrix import writeMatrix
factor = []
a = [[36.0, 12.0, 30.0, 6.0],
[12.0, 20.0, 2.0, 10.0],
[30.0, 2.0, 29.0, 1.0],
[6.0, 10.0, 1.0, 14.0]]
b = (18.0, 22.0, 7.0, 20.0)
# Factor A
linSolNonnegdef(a, b,
factor=factor,
factorOnly=True)
writeMatrix("factor", factor)
# Get the solution using
# the factorized matrix.
x = linSolNonnegdef(a, b,
factor=factor,
solveOnly=True)
writeMatrix("Solution, x, of Ax = b", x)
Output¶
factor
1 2 3 4
1 6 2 5 1
2 2 4 -2 2
3 5 -2 0 0
4 1 2 0 3
Solution, x, of Ax = b
1 2 3 4
0.167 0.500 0.000 1.000
Example 3¶
This example uses the inverse option to compute the symmetric g
inverse of the symmetric nonnegative matrix in the first example. Maindonald
(1984, p. 106) discusses the computations for this problem.
from numpy import *
from pyimsl.math.linSolNonnegdef import linSolNonnegdef
from pyimsl.math.matMulRect import matMulRect
from pyimsl.math.writeMatrix import writeMatrix
n = 4
p_inva = empty(0)
a = array([[36.0, 12.0, 30.0, 6.0],
[12.0, 20.0, 2.0, 10.0],
[30.0, 2.0, 29.0, 1.0],
[6.0, 10.0, 1.0, 14.0]])
# Get g2_inverse(a)
linSolNonnegdef(a, None, inverse=p_inva, inverseOnly=True)
# Form a*g2_inverse(a)
p_a_inva = matMulRect("A*B", aMatrix=a, bMatrix=p_inva)
# Form a*g2_inverse(a)*a
p_a_inva_a = matMulRect("A*B", aMatrix=p_a_inva, bMatrix=a)
writeMatrix("The g2 inverse of a", p_inva)
writeMatrix("a*g2_inverse(a)\nviolates condition 3 of the M-P inverse",
p_a_inva)
writeMatrix("a = a*g2_inverse(a)*a\ncondition 1 of the M-P inverse",
p_a_inva_a)
Output¶
The g2 inverse of a
1 2 3 4
1 0.0347 -0.0208 0.0000 0.0000
2 -0.0208 0.0903 0.0000 -0.0556
3 0.0000 0.0000 0.0000 0.0000
4 0.0000 -0.0556 0.0000 0.1111
a*g2_inverse(a)
violates condition 3 of the M-P inverse
1 2 3 4
1 1.0 -0.0 0.0 0.0
2 0.0 1.0 0.0 0.0
3 1.0 -0.5 0.0 0.0
4 0.0 0.0 0.0 1.0
a = a*g2_inverse(a)*a
condition 1 of the M-P inverse
1 2 3 4
1 36 12 30 6
2 12 20 2 10
3 30 2 29 1
4 6 10 1 14
Warning Errors¶
IMSL_INCONSISTENT_EQUATIONS_2 |
The linear system of equations is inconsistent. |
IMSL_NOT_NONNEG_DEFINITE |
The matrix A is not nonnegative definite. |