Log Determinant Evaluation in casadi

[sanket...@gmail.com]

Hi,

I a currenlty using casadi v3.0.0-rc2 from python 2.7. I wanted to evaluate the log determinant for an MX expression (for computing likelihood in a bayesian learning problem), but I end up getting an error:

" RuntimeError: MXNode::eval not defined for class N6casadi11DeterminantE "

If I convert all my symbolic matrices to SX class symbols, the expression evaluates fine for small matrices, but as the matrix size increases the time taken to declare the log likelihood function increases manyfold and the code gets stuck at the line where the expression is being declared as "import casadi as cs;  loglikehood = cs.log(cs.det(Matrix)) + ...;

This declaration does not finish execution (the code is stuck) even for a 10 x 10 matrix when declared through a SX symbolics.

Is there a better way to do this in casadi.

Best,

Sanket

[Joel Andersson]

Hi Sanket!

Calculating determinants has indeed not been implemented for MX and for SX, it will use minor expansion, which scales poorly, as you noticed. Indeed, the main use of the "det" function is to create a big expression to test the AD framework..

For SX, a more efficient way should be to use a QR decomposition. CasADi implements a QR factorization (with element reordering but without partial pivoting) that often works well:

So if you want to calculate log(det(A)), you can do:

[Q, R] = qr(A);

det(A) = det(Q)*det(R), which can be calculated cheaply. det(Q)=1 (it's orthogonal) and since R is triangular,  det(R) is just the multiplication of the diagonal elements. So:

log(det(A)) = log(det(R)) = trace(log(R))

This all assumes that A is positive definite etc.

I'm not sure if this is the smartest way forward though, there might be reformulations that you can use.

Joel

[sanket...@gmail.com]

Hi Joel,

Thanks for your quick reply, the QR decomposition worked well for my case. I will look into some reformulations next.

Best,

Sanket

[Gregor Reich]

Dear all

Has this been integrated into Casadi by now, or is there any plan to do it?

Best, Gregor.

[Gregor Reich]

Or asked differently: Is there any way to compute the determinant (or the log thereof) of an MX matrix? Currently, it looks like det itself cannot be evaluated, and neither qr nor chol can be done on MX (to apply the trace trick).

Thank and best regards, Gregor.

[Joel Andersson]

That trick also works for MX. But you need to encapsulate the function into an Function instance.

A = SX.sym('A',n, n)
log_det = ...
F = Function('F',[A],[log_det])

It has not yet been implemented in MX and I'm not aware of plans to do so. That being said, it's not a very hard thing to do.

Joel

[Gregor Reich]

Hi Joel

Thanks for the reply!

However, I'm sorry but I still don't get it: If I use your code below using MX (with and without qr), I get the same error:

F(magic(2))
Error using casadi.Function/call (line 735)
.../casadi/core/mx_node.cpp:291: 'eval' not defined for class N6casadi11DeterminantE

Error in casadi.Function/paren (line 1434)
        res = self.call(varargin);

Error in casadi.Function/subsref (line 1414)
        [varargout{1:nargout}]= paren(self, s.subs{:});

Am I doing it wrong?

Thanks again and best, Gregor.

[Joel Andersson]

You will only get that error message if you use the "det" operator. You should be able to avoid it by using SX and QR, as described in this thread.

Joel

[Gregor Reich]

Thanks for the quick reply, Joel!

I see now: By encapsulating it into a Function object, I can still call it with an MX argument (or whatever I want) and thus incorporate it into my MX-based code, although the definition of this particular function uses SX.

Great, problem solved!

Best, Gregor.