R devel - Apr 2010 - Lapack, determinant, multivariate normal density, solution to linear system, C language

r-devel list,

I have recently written an R package that solves a linear least squares
problem, and computes the multivariate normal density function. The bulk
of the code is written in C, with interfacing code to the BLAS and
Lapack libraries. The motivation here is speed. I ran into a problem
computing the determinant of a symmetric matrix in packed storage.
Apparently, there are no explicit routines for this as part of Lapack.
While there IS an explicit routine for this in Linpack, I did not want
to use the older library. Also, right before I needed the determinant of
the matrix A, I had used the Lapack routine dspsv to solve the linear
system Ax=b, which does much of the work of computing a determinant
also. In fact, the solution I came up with involves the output of this
routine (which might be obvious to Lapack designers, but not me)

My modest Googleing turned up very little unique material (as is typical
with BLAS/Lapack/Linpack queries). Hence, I am writing the r-devel list
partly to document the solution I've come up with, but mainly to elicit
additional wisdom from seasoned R programmers.

My solution to the problem is illustrated in the appended discussion and
C code. Thanks for your input.

-Matt Shotwell

--------------

The Lapack routine dspsv solves the linear system of equations Ax=b,
where A is a symmetric matrix in packed storage format. 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 1x1 or 2x2. 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 1x1, then the
determinant of the block is simply the value of the single element in
the block. If the diagonal block is of dimension 2x2 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.

  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 where 
  ** A represents a qxq matrix, b a 1xq vector.
  ** dspsv outputs the elements of the matrix D
  ** is upper triangular packed storage
  ** in the memory addressed by A. That is, A is
  ** replaced by D when dspsv returns.
  */
  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 the upper triangle 
  ** of 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) ];
    }
  }

On Mon, Apr 12, 2010 at 10:27 PM, shotwelm <shotwelm at musc.edu>
wrote:
> r-devel list,
>
> I have recently written an R package that solves a linear least squares
> problem, and computes the multivariate normal density function.
For both of those applications you can use a Cholesky decomposition of
the symmetric matrix.  If the Cholesky decomposition fails then you
have a singular least squares problem or a singular
variance-covariance matrix for your multivariate normal density
function.

Have you tried comparing the speed of your code to

prod(diag(chol(mm))^2

or, probably better, is to use the logarithm of the determinant

2 * sum(log(diag(chol(mm)))

If you use the Matrix package class dpoMatrix to solve the linear
system it will cache the results of the Cholesky decomposition when
solving the system so later evaluation of the determinant will be very
fast - although I suspect you would need to be working with matrices
of sizes in the hundreds or doing the same operation thousands of
times before you would notice a difference.

If you really insist on doing this in compiled code you just need to
call F77_CALL(dpotrf) then accumulate the product of the diagonal
elements of the resulting factor.

You could use packed storage but the slight advantage in memory usage
(at best, 1/2 of the full storage usage) is not worth the pain of
writing code to navigate the packed storage locations.
> The bulk
> of the code is written in C, with interfacing code to the BLAS and
> Lapack libraries. The motivation here is speed. I ran into a problem
> computing the determinant of a symmetric matrix in packed storage.
> Apparently, there are no explicit routines for this as part of Lapack.
> While there IS an explicit routine for this in Linpack, I did not want
> to use the older library. Also, right before I needed the determinant of
> the matrix A, I had used the Lapack routine dspsv to solve the linear
> system Ax=b, which does much of the work of computing a determinant
> also. In fact, the solution I came up with involves the output of this
> routine (which might be obvious to Lapack designers, but not me)
>
> My modest Googleing turned up very little unique material (as is typical
> with BLAS/Lapack/Linpack queries). Hence, I am writing the r-devel list
> partly to document the solution I've come up with, but mainly to elicit
> additional wisdom from seasoned R programmers.
>
> My solution to the problem is illustrated in the appended discussion and
> C code. Thanks for your input.
>
> -Matt Shotwell
>
> --------------
>
> The Lapack routine dspsv solves the linear system of equations Ax=b,
> where A is a symmetric matrix in packed storage format. 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 1x1 or 2x2. 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 1x1, then the
> determinant of the block is simply the value of the single element in
> the block. If the diagonal block is of dimension 2x2 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.
>
> ?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 where
> ?** A represents a qxq matrix, b a 1xq vector.
> ?** dspsv outputs the elements of the matrix D
> ?** is upper triangular packed storage
> ?** in the memory addressed by A. That is, A is
> ?** replaced by D when dspsv returns.
> ?*/
> ?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 the upper triangle
> ?** of 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) ];
> ? ?}
> ?}
>
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>