Prev Next cppad_det_minor.cpp

CppAD Speed: Gradient of Determinant Using Expansion by Minors

Operation Sequence
Note that the expansion by minors operation sequence does not depend on the matrix being factored. Hence we use the same ADFun object for all the matrices.

compute_det_minor
Routine that computes the gradient of determinant using CppAD:
# include <cppad/vector.hpp>
# include <cppad/speed/det_by_minor.hpp>
# include <cppad/speed/uniform_01.hpp>

void compute_det_minor(
     size_t                     size     , 
     size_t                     repeat   , 
     CppAD::vector<double>     &matrix   ,
     CppAD::vector<double>     &gradient )
{
     // -----------------------------------------------------
     // setup

     // object for computing determinant
     typedef CppAD::AD<double>       ADScalar; 
     typedef CppAD::vector<ADScalar> ADVector; 
     CppAD::det_by_minor<ADScalar>   Det(size);

     size_t i;               // temporary index
     size_t m = 1;           // number of dependent variables
     size_t n = size * size; // number of independent variables
     ADVector   A(n);        // AD domain space vector
     ADVector   detA(m);     // AD range space vector
     
     // vectors of reverse mode weights 
     CppAD::vector<double> w(1);
     w[0] = 1.;

     // choose a matrix
     CppAD::uniform_01(n, matrix);
     for( i = 0; i < size * size; i++)
          A[i] = matrix[i];

     // declare independent variables
     Independent(A);

     // AD computation of the determinant
     detA[0] = Det(A);

     // create function object f : A -> detA
     CppAD::ADFun<double> f(A, detA);

     // ------------------------------------------------------
     while(repeat--)
     {    // get the next matrix
          CppAD::uniform_01(n, matrix);

          // evaluate the determinant at the new matrix value
          f.Forward(0, matrix);

          // evaluate and return gradient using reverse mode
          gradient = f.Reverse(1, w);
     }
     return;
}

Input File: speed/cppad/det_minor.cpp