Matrix determinant with the Lapack routine dspsv

April 6, 2010
By

(This article was first published on BioStatMatt » R, and kindly contributed to R-bloggers)

The Lapack routine dspsv solves the linear system of equations Ax=b, where A is a symmetric matrix in packed storage format. However, there appear to be no Lapack functions that compute the determinant of such a matrix. We need to compute the determinant, for instance, in order to compute the multivariate normal density function. The dspsv function performs the factorization A=UDU’, where U is a unitriangular matrix and D is a block diagonal matrix where the blocks are of dimension 1×1 or 2×2. In addition to the solution for x, the dspsv function also returns the matrices U and D. The matrix D may then be used to compute the determinant of A. Recall from linear algebra that det(A) = det(UDU’) = det(U)det(D)det(U’). Since U is unitriangular, det(U) = 1. The determinant of D is the product of the determinants of the diagonal blocks. If a diagonal block is of dimension 1×1, then the determinant of the block is simply the value of the single element in the block. If the diagonal block is of dimension 2×2 then the determinant of the block may be computed according to the well-known formula b11*b22-b12*b21, where bij is the value in the i’th row and j’th column of the block. The following C code snip demonstrates the procedure.

int i, q, info, *ipiv, one = 1;
double *b, *A, *D, det;

/*
** A and D are upper triangular matrices in packed storage
** A[] = a00, a01, a11, a02, a12, a22, a03, a13, a23, a33, ...
** use the following macro to address the element in the
** i'th row and j'th column for i <= j
*/
#define UMAT(i, j) (i + j * ( j + 1 ) / 2)

/*
** additional code should be here
** - set q
** - allocate ipiv...
** - allocate and fill A and b...
*/

/*
** solve Ax=B using A=UDU' factorization, D is placed in A
** where A represents a qxq matrix, b a 1xq vector
*/
F77_CALL(dspsv)("U", &q, &one, A, ipiv, b, &q, &info);
if( info > 0 ) { /*issue warning, system is singular*/ }
if( info < 0 ) { /*issue error, invalid argument*/ }

/*
** compute the determinant det = det(A)
** if ipiv[i] > 0, then D(i,i) is a 1x1 block diagonal
** if ipiv[i] = ipiv[i-1] < 0, then D(i-1,i-1),
** D(i-1,i), and D(i,i) form a 2x2 block diagonal
*/
D = A;
det = 1.0;
for( i = 0; i < q; i++ ) {
  if( ipiv[ i ] > 0 ) {
    det *= D[ UMAT(i,i) ];
  } else if( i > 0 && ipiv[ i ] < 0 && ipiv[ i-1 ] == ipiv[ i ] ) {
    det *= D[ UMAT(i,i) ] * D[ UMAT(i-1,i-1) ] -\
           D[ UMAT(i-1,i) ] * D[ UMAT(i-1,i) ];
  }
}

To leave a comment for the author, please follow the link and comment on his blog: BioStatMatt » R.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Tags: , ,

Comments are closed.