Solves a set of positive definite linear equations.
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
x (NxK matrix or M-dimensional array) – where the last two dimensions are NxK, the solutions for the system of equations, \(Ax = b\).
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
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\).
If the A matrix is not positive definite:
return scalar error code 30.
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
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);