numpy/scipy.linalg.solve with MX RHS

bransfor
2012-06-29
2013-04-11
  • bransfor

    bransfor - 2012-06-29

    I would like to multiply the inverse of a (potentially large) numpy array by a casadi MX object using the python interface. The goal is to obtain a transformed MX object. Is there some way to do this using numpy/scipy.linalg routines? The following code fails when b is converted to a numpy array by nump.linalg.solve.

    import numpy
    import casadi
    A = numpy.eye(3)
    b = casadi.MX('b', 3)
    x = numpy.linalg.solve(A, b)
    Traceback (most recent call last):
      File "<stdin>", line 1, in <module>
      File "C:\Python27\lib\site-packages\numpy\linalg\linalg.py", line 307, in solve
        b, wrap = _makearray(b)
      File "C:\Python27\lib\site-packages\numpy\linalg\linalg.py", line 67, in _makearray
        new = asarray(a)
      File "C:\Python27\lib\site-packages\numpy\core\numeric.py", line 235, in asarray
        return array(a, dtype, copy=False, order=order)
      File "C:\Python27\lib\site-packages\casadi\casadi.py", line 10197, in __array__
        return self.__array_custom__(*args,**kwargs)
      File "C:\Python27\lib\site-packages\casadi\casadi.py", line 10203, in __array_custom__
        raise Exception("MX cannot be converted to an array. MX.__array__ purely exists to allow ufunc/numpy goodies")
    Exception: MX cannot be converted to an array. MX.__array__ purely exists to allow ufunc/numpy goodies
    

    I actually need to do this operation many times with the same matrix but different right hands sides, so I would prefer to use scipy.linalg.cho_factor to factorize the matrix and then run scipy.linalg.cho_solve to evaluate the operation.

    Thanks!

     
  • Joel Andersson

    Joel Andersson - 2012-06-29

    Hello!

    The MX datatype cannot be used for this purpose. MX is not a matrix where the entries are expressions, but an expression which is matrix-valued. This is different from SXMatrix, which is indeed a matrix with entries that are expressions.

    I would first of all advice you to use the SXMatrix datatype instead of MX. Just replace "b = casadi.MX('b', 3)" (which is equivalent to "b = casadi.msym('b', 3)") with "b = casadi.ssym('b', 3)". MX will eventually support "backslash nodes" and use the matrix chain rule for these in the automatic differentiation, but this is not currently implemented. Until then, I think it's better to use SXMatrix.

    In any case, I am very skeptical that scipy.linalg.cho_factor will work well with any symbolic datatype (SX in this case). The reason is that operations such as Cholesky factorization relies on the ability to compare if one value is larger than another (for example in the partial pivoting), which is always possible with a numeric datatype (e.g. float), but not with a symbolic datatype. If you do get it to work, I would be very interested.

    What is much more likely to work is CasAD's native QR-factorization, usage  = casadi.qr(A). This contains a simple implementation of Gram-Schmidt orgonalization _without_ partial pivoting, but with a preparation step which ensures that structural zeros do not appear along the diagonal (using a Dulmage-Mendelsohn decomposition). Since there is no pivoting, it might fail under certain circumstances. I guess it should be easy to implement Cholesky factorization (which should be more efficient) the same way, though I imagine that QR is a safer bet if you cannot pivot.

    Note that you can form an SXFunction out of the expression that you got during the symbolic QR factorization and include calls to this function in MX expression graphs.

    Regards,
    Joel

     
  • bransfor

    bransfor - 2012-06-29

    Thanks, Joel. I didn't understand the difference between MX and SXMatrix. You comments made it clear.

     
  • Anonymous - 2012-06-29

    I'd also be interested in getting this to work - are there any examples of solving Ax=b using the symbolic QR factorization?

     
  • Joel Andersson

    Joel Andersson - 2012-06-30

    The function "solve" does this, usage: x = casadi.solve(A,b). You can look into the implementation of this function.

    Best,
    Joel

     
  • bransfor

    bransfor - 2012-07-01

    Can the elements of the RHS SXMatrix point to views from an MX object? The following works when line 3 starts with "if 1" but fails when line 3 starts with "if 0"

    import numpy, casadi
    b = casadi.SXMatrix(3, 1)
    if 1:
        y = [casadi.SX(str(i)) for i in range(10)]
    else:    
        y = casadi.MX('y', 10)
    b[0, 0] = y[0]
    b[1, 0] = y[1]
    b[2, 0] = y[2]
    A = numpy.eye(3)
    x = casadi.solve(A, b)
    
     
  • Joel Andersson

    Joel Andersson - 2012-07-02

    You cannot mix SXMatrix and MX like this. The relation between SX and MX is described in the users guide. I suggest you to read that.

    Futhermore, in your code, you create variables named "0" and "1", etc. Is this really what you want to do?

    I think that the code you want to write looks as follows:

    import numpy, casadi
    if 1:
        b = casadi.SXMatrix.zeros(3, 1)
        y = casadi.ssym('y',10)
    else:    
        b = casadi.MX.zeros(3, 1)
        y = casadi.msym('y',10)
    b[0, 0] = y[0]
    b[1, 0] = y[1]
    b[2, 0] = y[2]
    A = numpy.eye(3)
    x = casadi.solve(A, b)
    

    This is valid syntax, but will not work with "if 0" since "x = casadi.solve(A, b)" is not yet supported for the MX datatype. Note that MX does support index access and assignments, but the way it works internally is that an "assignment expression" is created.

    But why would you want to work with the MX-datatype in this case? SX will do the job and should be much faster.

    Best,
    Joel

     
  • bransfor

    bransfor - 2012-07-02

    Thanks for the comments, Joel. I will review the documentation.

    The code I posted is just for testing purposes. The actual application is a direct collocation NLP problem solved with IPOPT. What I wish to do is place a multivariate normal prior on a subset of the decision variables. Referring to the example code, assume that the first three elements of the decision variable vector ‘y’ are parameters in the model. The vector b is created to point to the first three elements of y. To the objective function I want to add

    prior = 0.5 * casadi.mul(b.T, casadi.solve(A, b)) + numpy.log(numpy.linalg.det(A))
    

    where A is the parameter covariance matrix. Up until now the decision variable vector in the code was an MX, but now I have updated it to use an ssym.

    Thanks

     
  • Joel Andersson

    Joel Andersson - 2012-07-02

    Hello!

    The covariance matrix A, is it known? I mean, do you have it's numerical values? In that case, I guess it makes the most sense that you factorize this matrix using numpy/scipy. casadi.solve, really only makes sense if "A" is contains symbolic expressions.

    I guess you can use whatever numerical routine to factorize it, and then just perform the backsolve in CasADi (also using casadi.solve). This should be a lot more efficient and stable.

    Joel

     
  • bransfor

    bransfor - 2012-07-02

    Yes, A is know. It's just a numpy array.

     
  • Joel Andersson

    Joel Andersson - 2012-07-02

    Oh, OK. Then it's probably much smarter to use scipy routines to factorize the matrix. Maybe even solving the linear system will work then. Did you try?

     
  • bransfor

    bransfor - 2012-07-02

    Yes, I am using scipy.linalg.cho_factor to factorize the matrix and casadi.solve to solve the linear systems. Here is an example of the relevant steps

    import numpy as np
    import casadi
    from scipy import linalg
    class Prior(object):
        def __init__(self, mu, covariance):
            self._mu = mu
            self._cho_upper, _ = linalg.cho_factor(covariance)
            self._log_det = 2.0 * np.sum(np.log(np.diag(self._cho_upper)))  
        def negative_log_likelihood(self, x):   
            'x is a casadi.SXMatrix'
            dx = x - self._mu
            y = casadi.solve(self._cho_upper.T, dx)    
            y = casadi.solve(self._cho_upper, y)    
            like = 0.5 * casadi.mul(dx.T, y) + self._log_det
            return like.toScalar()
    
     
  • Joel Andersson

    Joel Andersson - 2012-07-02

    OK. Nice! Looks like a good approach.

     

Log in to post a comment.

Get latest updates about Open Source Projects, Conferences and News.

Sign up for the SourceForge newsletter:





No, thanks