coin-or / qpOASES

Open-source C++ implementation of the recently proposed online active set strategy
GNU Lesser General Public License v2.1
397 stars 132 forks source link

Eigen for linear algebra #116

Open jlack1987 opened 3 years ago

jlack1987 commented 3 years ago

I'm not sure if this is a desire for anybody other than myself or if the maintainers are interested in this or not so feel free to "close won't do" if i'm alone in this, but it would really make life easier if Eigen was used for linear algebra. Especially for the more advanced features like the sparse utilizing sparse matrices.

jhallier commented 3 years ago

I would second this, and could support in implementation. Using both Eigen::SparseMatrix and qpOASES::SparseMatrix in client code becomes a hassle, when no easy convert function exists. It would be great to have an overload of all init and hotstart methods that accepts at least the Hessian H and constraint matrix A as Eigen::Matrix or Eigen::SparseMatrix.

apotschka commented 3 years ago

It would be great to have an easier interface for users of Eigen. However, a compile-time dependency of qpOASES on Eigen would probably create a lot of headache for those that need to run qpOASES on embedded devices.

My feeling is that conversion shouldn't be too hard, as the data structures for matrices in qpOASES are completely standard.

If one would like a more native support without deep copies, I think deriving Eigen-based classes from the abstract matrix classes in qpOASES should be doable (but I don't have the time to figure out the details).

Would it be a viable solution to start a separate small project that delivers the conversion routines and the derived matrix classes? In this case, current users wouldn't be bothered with the new dependencies on Eigen.

jlack1987 commented 3 years ago

yeah perhaps, it's really the more advanced feature sets where I have run into issues that have taken me great amounts of time/effort to support. For the simple init/hotstart case where no advanced features are used one can simply make these two typedefs,

typedef Eigen::Matrix<qpOASES::real_t, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> QPMatrix;
typedef Eigen::Matrix<qpOASES::real_t, Eigen::Dynamic, 1> QPVector;

and say I then have QPMatrix H, I can feed that into init/hotstart as H.data() which returns a pointer to the underlaying matrix/vector data and that all works fine and no difficult type conversions are necessary. I ran into real issues though when trying to use Eigens sparse matrix stuff and then trying to use that with the qpOASES sparse matrix API. I guess perhaps I should remake the ticket and maybe more specifically point out the sparse matrix pain points, because since Eigens .data() gives a pointer to the underlying array there's really not much reason to add a full dependency on Eigen.

jhallier commented 3 years ago

Yes, I also referred to Eigen::SparseMatrix. I understand that the dependency on Eigen would be not desired.

Currently, I'm using a naive mapping, such as let's say I have an Eigen::SparseMatrix Eigen_sparse

std::vector<qpOASES::real_t> values;
std::vector<qpOASES::sparse_int_t> row_indices;
std::vector<qpOASES::sparse_int_t> col_pointer;
qpOASES::int_t colpointer = 0;
for (int k = 0; k < Eigen_sparse.outerSize(); ++k) {
    col_pointer.push_back(colpointer);
    for (Eigen::SparseMatrix<Float64>::InnerIterator it(Eigen_sparse, k); it; ++it) {
        values.push_back(it.value());
        row_indices.push_back(it.row());
        ++colpointer;
    }
}
col_pointer.push_back(colpointer);

And then I construct the qpOASES::SparseMatrix from the three vectors values, row_indices and col_pointer. It's not very pretty though. Perhaps someone here has a better solution? What may be interesting if the constructor for qpOASES::SparseMatrix would get another overload that accepts Eigen style triplets (row, column, value). It is always then necessary when Eigen matrices are created that are used for calculations, and the result of which is a matrix for the qpOASES solver, H or A.

jhallier commented 3 years ago

I just accidentally found that it is a lot easier than that, perhaps it is useful to someone else: If you have a Eigen::SparseMatrix, you can easily transform it to a qpOASES::SparseMatrix by directly using the CCS-storage pointers of Eigen.

Example:

Eigen;:SparseMatrix<double> eigen_sparse;
qpOASES::SparseMatrix qpoases_sparse(eigen_sparse.rows(), eigen_sparse.cols(), 
        eigen_sparse.innerIndexPtr(), eigen_sparse.outerIndexPtr(), eigen_sparse.valuePtr());
qpoases_sparse.createDiagInfo();

Just make sure that eigen_sparse is still valid when qpoases_sparse is being used, e.g. by creating it as a class member via a unique pointer. What I don't really understand is why the method createDiagInfo() needs to be called. What exactly does it do?

jlack1987 commented 3 years ago

I just accidentally found that it is a lot easier than that, perhaps it is useful to someone else: If you have a Eigen::SparseMatrix, you can easily transform it to a qpOASES::SparseMatrix by directly using the CCS-storage pointers of Eigen.

Example:

Eigen;:SparseMatrix<double> eigen_sparse;
qpOASES::SparseMatrix qpoases_sparse(eigen_sparse.rows(), eigen_sparse.cols(), 
        eigen_sparse.innerIndexPtr(), eigen_sparse.outerIndexPtr(), eigen_sparse.valuePtr());
qpoases_sparse.createDiagInfo();

Just make sure that eigen_sparse is still valid when qpoases_sparse is being used, e.g. by creating it as a class member via a unique pointer. What I don't really understand is why the method createDiagInfo() needs to be called. What exactly does it do?

NIce, I will need to give that a try

apotschka commented 3 years ago

What I don't really understand is why the method createDiagInfo() needs to be called. What exactly does it do?

It allocates an index array jd, which points to the diagonal elements of a sparse matrix. It is useful for fast access and modification of the diagonal elements, e.g., in SparseMatrix::addToDiag.