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.