solpd#

Purpose#

Solves a set of positive definite linear equations.

Format#

x = solpd(b, A)#

Parameters:

b (matrix or array) – NxK matrix or M-dimensional array where the last two dimensions are NxK
A (matrix or array) – NxN symmetric positive definite matrix or M-dimensional array where the NxN 2-dimensional arrays described by the last two dimensions are symmetric and positive definite

Returns:

x (NxK matrix or M-dimensional array) – where the last two dimensions are NxK, the solutions for the system of equations, \(Ax = b\).

Examples#

n = 5;
format /lo 16,8;

A = rndn(n, n);
A = A'A;
x = rndn(n, 1);
b = A*x;

x2 = solpd(b, A);

print " X solpd(b, A) Difference";
print x~x2~x-x2;

X solpd(b, A) Difference

 0.32547881   0.32547881  -4.9960036e-16
 1.5190182    1.5190182   -1.7763568e-15
 0.88099266   0.88099266   1.5543122e-15
 1.8192784    1.8192784   -2.2204460e-16
-0.060848175 -0.060848175 -1.4710455e-15

Remarks#

b can have more than one column. If so, the system of equations is solved for each column, i.e., \(A*x[., i] = b[., i]\).

This function uses the Cholesky decomposition to solve the system directly. Therefore it is more efficient than using \(inv(A)*b\).

If b and A are M-dimensional arrays, the sizes of their corresponding M-2 leading dimensions must be the same. The resulting array will contain the solutions for the system of equations given by each of the corresponding 2-dimensional arrays described by the two trailing dimensions of b and A. In other words, for a 10x4x2 array b and a 10x4x4 array A, the resulting array x will contain the solutions for each of the 10 corresponding 4x2 arrays contained in b and 4x4 arrays contained in A. Therefore, \(A[n,.,.]*x[n,.,.] = b[n,.,.]\), for \(1 ≤ n ≤ 10\).

solpd() does not check to see that the matrix A is symmetric. solpd() will look only at the upper half of the matrix including the principal diagonal.

If the A matrix is not positive definite:

`trap 1`	return scalar error code 30.
`trap 0`	terminate with an error message.

One obvious use for this function is to solve for least squares coefficients. The effect of this function is thus similar to that of the / operator.

If X is a matrix of independent variables, and Y is a vector containing the dependent variable, then the following code will compute the least squares coefficients of the regression of Y on X:

b = solpd(X'Y, X'X);