[CppAD] LuSolve nan

[Brad Bell]

This is a subtitle issue related to the theory of Algorithmic 
Differentiation; namely a precise definition of what gets taped and what 
does not get taped.

For the case below, Float is AD<double> and logdet is computing the log 
of the absolute value of the determinant and signdet is the sign of the 
determinant. This is accomplished int LuSolve code by the statements:
                // update the determinant
                if( LeqZero ( pivot ) )
                {       logdet += log( - pivot );
                        signdet = - signdet;
                }
                else    logdet += log( pivot );

Now comes the actual problem. During the taping of the function LeqZero 
returns false, and the evaluation of LeqZero is not taped. It is 
possible to use CondExp to get the evaluation the LeqZero comparison 
into the tape. But there is a more general problem here; namely that the 
choice of pivot elements is not being taped. Even more generally, one 
needs to consider if there are other conditional statements in ones 
algorithm that have different choices depending on the value of the 
independent variables.

The example function
    http://www.coin-or.org/CppAD/Doc/luvecad.xml
actually tapes all of the pivots as well as the evaluation of the 
determinant so that it will work for all different arguments.

Thus you have two choices, one is the use LuVecAD and the other is to 
retape when you need to (By retape I mean create a new function argument 
corresponding to the values for the independent variables). I have been 
considering adding an argument to LuSolve that indicates when one should 
retape, but have not yet done so.

PS.
The file
    Speed/LuSolve.cpp
contains a speed comparison between retaping and using LuVecAD.

Kasper Kristensen wrote:

>Hello all,
>
>I have attached a small code example which tapes the mapping
>A->log(abs(det(A))) for a 1x1 matrix A=1 using LuSolve.
>I was expecting the tape to be valid for all non-zero matrices A.
>However if I evaluate the tape with a matrix for which det(A)<0 the
>evaluation gives nan (but the Jacobian is correct!).
>
>Best regards
>Kasper
>
>PROGRAM OUTPUT:
>
>A=2
>{ 0.693147 }{ 0.5 }
>A=-2
>{ nan }{ -0.5 }
>
>CODE:
>
># include <CppAD/CppAD.h>
>using namespace CppAD;
>int main()
>{
>   // dimension of the matrix
>   size_t n = 1;
>
>   // Dummy arguments to LuSolve
>   vector<AD<double> > B(0);
>   vector<AD<double> > X(0);
>
>   // construct ADFun object F:A->log(abs(det(A)))
>   vector<AD<double> > A(n);
>   vector<AD<double> > y(n);
>   A[0]=1;
>   Independent(A);
>   AD<double> logdet;
>   LuSolve(n, 0, A, B, X, logdet);
>   y[0]=logdet;
>   ADFun<double> F(A,y);
>
>   // evaluate the tape
>   vector<double> A0(n);
>   A0[0]=2;
>   std::cout << "A=2\n";
>   std::cout << F.Forward(0,A0) ;
>   std::cout << F.Jacobian(A0) ;
>   A0[0]=-2;
>   std::cout << "\nA=-2\n";
>   std::cout << F.Forward(0,A0) ;
>   std::cout << F.Jacobian(A0) << "\n";
>
>}
>
>
>_______________________________________________
>CppAD mailing list
>CppAD at list.coin-or.org
>http://list.coin-or.org/mailman/listinfo/cppad
>
>
>  
>