Program Listing for File random_qp_problems.hpp
↰ Return to documentation for file (/tmp/ws/src/proxsuite/include/proxsuite/proxqp/utils/random_qp_problems.hpp
)
#ifndef PROXSUITE_PROXQP_UTILS_RANDOM_QP_PROBLEMS_HPP
#define PROXSUITE_PROXQP_UTILS_RANDOM_QP_PROBLEMS_HPP
#include <Eigen/Core>
#include <Eigen/SparseCore>
#include <Eigen/Cholesky>
#include <Eigen/Eigenvalues>
#include <Eigen/QR>
#include <utility>
#include <proxsuite/proxqp/dense/views.hpp>
#include <proxsuite/proxqp/dense/model.hpp>
#include <proxsuite/proxqp/sparse/model.hpp>
#include <map>
#include <random>
namespace proxsuite {
namespace proxqp {
namespace utils {
using c_int = long long;
using c_float = double;
namespace proxqp = proxsuite::proxqp;
template<typename T, proxqp::Layout L>
using Mat =
Eigen::Matrix<T,
Eigen::Dynamic,
Eigen::Dynamic,
(L == proxqp::colmajor) ? Eigen::ColMajor : Eigen::RowMajor>;
template<typename T>
using Vec = Eigen::Matrix<T, Eigen::Dynamic, 1>;
template<typename Scalar>
using SparseMat = Eigen::SparseMatrix<Scalar, Eigen::ColMajor, c_int>;
using namespace proxqp;
namespace eigen {
template<typename T>
void
llt_compute( //
Eigen::LLT<T>& out,
T const& mat)
{
out.compute(mat);
}
template<typename T>
void
ldlt_compute( //
Eigen::LDLT<T>& out,
T const& mat)
{
out.compute(mat);
}
LDLT_EXPLICIT_TPL_DECL(2, llt_compute<Mat<f32, colmajor>>);
LDLT_EXPLICIT_TPL_DECL(2, ldlt_compute<Mat<f32, colmajor>>);
LDLT_EXPLICIT_TPL_DECL(2, llt_compute<Mat<f32, rowmajor>>);
LDLT_EXPLICIT_TPL_DECL(2, ldlt_compute<Mat<f32, rowmajor>>);
LDLT_EXPLICIT_TPL_DECL(2, llt_compute<Mat<f64, colmajor>>);
LDLT_EXPLICIT_TPL_DECL(2, ldlt_compute<Mat<f64, colmajor>>);
LDLT_EXPLICIT_TPL_DECL(2, llt_compute<Mat<f64, rowmajor>>);
LDLT_EXPLICIT_TPL_DECL(2, ldlt_compute<Mat<f64, rowmajor>>);
} // namespace eigen
namespace rand {
using proxqp::u32;
using proxqp::u64;
#ifdef _MSC_VER
/* Using the MSCV compiler on Windows causes problems because the type uint128
is not available. Therefore, we use a random number generator from the stdlib
instead of our custom Lehmer random number generator. The necessary lehmer
functions used in in our code are remplaced with calls to the stdlib.*/
std::mt19937 gen(1234);
std::uniform_real_distribution<> uniform_dist(0.0, 1.0);
std::normal_distribution<double> normal_dist;
using u128 = u64;
inline auto
uniform_rand() -> double
{
double output = double(uniform_dist(gen));
return output;
}
inline auto
lehmer_global() -> u128&
{
static u64 output = gen();
return output;
}
inline void
set_seed(u64 seed)
{
gen.seed(seed);
}
inline auto
normal_rand() -> double
{
return normal_dist(gen);
}
#else
using u128 = __uint128_t;
constexpr u128 lehmer64_constant(0xda942042e4dd58b5);
inline auto
lehmer_global() -> u128&
{
static u128 g_lehmer64_state = lehmer64_constant * lehmer64_constant;
return g_lehmer64_state;
}
inline auto
lehmer64() -> u64
{ // [0, 2^64)
lehmer_global() *= lehmer64_constant;
return u64(lehmer_global() >> u128(64U));
}
inline void
set_seed(u64 seed)
{
lehmer_global() = u128(seed) + 1;
lehmer64();
lehmer64();
}
inline auto
uniform_rand() -> double
{ // [0, 2^53]
u64 a = lehmer64() / (1U << 11U);
return double(a) / double(u64(1) << 53U);
}
inline auto
normal_rand() -> double
{
static const double pi2 = std::atan(static_cast<double>(1)) * 8;
double u1 = uniform_rand();
double u2 = uniform_rand();
double ln = std::log(u1);
double sqrt = std::sqrt(-2 * ln);
return sqrt * std::cos(pi2 * u2);
}
#endif
template<typename Scalar>
auto
vector_rand(isize nrows) -> Vec<Scalar>
{
auto v = Vec<Scalar>(nrows);
for (isize i = 0; i < nrows; ++i) {
v(i) = Scalar(rand::normal_rand());
}
return v;
}
template<typename Scalar>
auto
matrix_rand(isize nrows, isize ncols) -> Mat<Scalar, colmajor>
{
auto v = Mat<Scalar, colmajor>(nrows, ncols);
for (isize i = 0; i < nrows; ++i) {
for (isize j = 0; j < ncols; ++j) {
v(i, j) = Scalar(rand::normal_rand());
}
}
return v;
}
namespace detail {
template<typename Scalar>
auto
orthonormal_rand_impl(isize n) -> Mat<Scalar, colmajor>
{
auto mat = rand::matrix_rand<Scalar>(n, n);
auto qr = mat.householderQr();
Mat<Scalar, colmajor> q = qr.householderQ();
return q;
}
using Input = std::pair<u128, isize>;
} // namespace detail
template<typename Scalar>
auto
orthonormal_rand(isize n) -> Mat<Scalar, colmajor> const&
{
static auto cache = std::map<detail::Input, Mat<Scalar, colmajor>>{};
auto input = detail::Input{ lehmer_global(), n };
auto it = cache.find(input);
if (it == cache.end()) {
auto res = cache.insert({
input,
detail::orthonormal_rand_impl<Scalar>(n),
});
it = res.first;
}
return (*it).second;
}
template<typename Scalar>
auto
positive_definite_rand(isize n, Scalar cond) -> Mat<Scalar, colmajor>
{
auto const& q = rand::orthonormal_rand<Scalar>(n);
auto d = Vec<Scalar>(n);
{
using std::exp;
using std::log;
Scalar diff = log(cond);
for (isize i = 0; i < n; ++i) {
d(i) = exp(Scalar(i) / Scalar(n) * diff);
}
}
return q * d.asDiagonal() * q.transpose();
}
template<typename Scalar>
auto
sparse_positive_definite_rand(isize n, Scalar cond, Scalar p)
-> SparseMat<Scalar>
{
auto H = SparseMat<Scalar>(n, n);
for (isize i = 0; i < n; ++i) {
auto random = Scalar(rand::normal_rand());
H.insert(i, i) = random;
}
for (isize i = 0; i < n; ++i) {
for (isize j = i + 1; j < n; ++j) {
auto urandom = rand::uniform_rand();
if (urandom < p / 2) {
auto random = Scalar(rand::normal_rand());
H.insert(i, j) = random;
}
}
}
Mat<Scalar, colmajor> H_dense = H.toDense();
Vec<Scalar> eigh =
H_dense.template selfadjointView<Eigen::Upper>().eigenvalues();
Scalar min = eigh.minCoeff();
Scalar max = eigh.maxCoeff();
// new_min = min + rho
// new_max = max + rho
//
// (max + rho)/(min + rho) = cond
// 1 + (max - min) / (min + rho) = cond
// (max - min) / (min + rho) = cond - 1
// min + rho = (max - min) / (cond - 1)
// rho = (max - min)/(cond - 1) - min
Scalar rho = (max - min) / (cond - 1) - min;
if (max == min) {
rho += 1;
}
for (isize i = 0; i < n; ++i) {
H.coeffRef(i, i) += rho;
}
H.makeCompressed();
return H;
}
template<typename Scalar>
auto
sparse_positive_definite_rand_compressed(isize n, Scalar rho, Scalar p)
-> SparseMat<Scalar>
{
auto H = SparseMat<Scalar>(n, n);
H.setZero();
for (isize i = 0; i < n; ++i) {
for (isize j = i + 1; j < n; ++j) {
auto urandom = rand::uniform_rand();
if (urandom < p) {
auto random = Scalar(rand::normal_rand());
H.insert(i, j) = random;
}
}
}
Mat<Scalar, colmajor> H_dense = H.toDense();
Vec<Scalar> eigh =
H_dense.template selfadjointView<Eigen::Upper>().eigenvalues();
Scalar min = eigh.minCoeff();
for (isize i = 0; i < n; ++i) {
H.coeffRef(i, i) += (rho + abs(min));
}
H.makeCompressed();
return H;
}
template<typename Scalar>
auto
sparse_positive_definite_rand_not_compressed(isize n, Scalar rho, Scalar p)
-> Mat<Scalar, colmajor>
{
auto H = Mat<Scalar, colmajor>(n, n);
H.setZero();
for (isize i = 0; i < n; ++i) {
for (isize j = 0; j < n; ++j) {
auto urandom = rand::uniform_rand();
if (urandom < p / 2) {
auto random = Scalar(rand::normal_rand());
H(i, j) = random;
}
}
}
H = ((H + H.transpose()) * 0.5)
.eval(); // safe no aliasing :
// https://eigen.tuxfamily.org/dox/group__TopicAliasing.html
// H.array() /= 2.;
Vec<Scalar> eigh = H.template selfadjointView<Eigen::Upper>().eigenvalues();
Scalar min = eigh.minCoeff();
H.diagonal().array() += (rho + abs(min));
return H;
}
template<typename Scalar>
auto
sparse_matrix_rand(isize nrows, isize ncols, Scalar p) -> SparseMat<Scalar>
{
auto A = SparseMat<Scalar>(nrows, ncols);
for (isize i = 0; i < nrows; ++i) {
for (isize j = 0; j < ncols; ++j) {
if (rand::uniform_rand() < p) {
A.insert(i, j) = Scalar(rand::normal_rand());
}
}
}
A.makeCompressed();
return A;
}
template<typename Scalar>
auto
sparse_matrix_rand_not_compressed(isize nrows, isize ncols, Scalar p)
-> Mat<Scalar, colmajor>
{
auto A = Mat<Scalar, colmajor>(nrows, ncols);
A.setZero();
for (isize i = 0; i < nrows; ++i) {
for (isize j = 0; j < ncols; ++j) {
if (rand::uniform_rand() < p) {
A(i, j) = Scalar(rand::normal_rand());
}
}
}
return A;
}
} // namespace rand
using proxqp::usize;
namespace osqp {
auto
to_sparse(Mat<c_float, colmajor> const& mat) -> SparseMat<c_float>;
auto
to_sparse_sym(Mat<c_float, colmajor> const& mat) -> SparseMat<c_float>;
} // namespace osqp
template<typename T>
auto
matmul_impl( //
Mat<T, proxqp::colmajor> const& lhs,
Mat<T, proxqp::colmajor> const& rhs) -> Mat<T, proxqp::colmajor>
{
return lhs.operator*(rhs);
}
template<typename To, typename From>
auto
mat_cast(Mat<From, proxqp::colmajor> const& from) -> Mat<To, proxqp::colmajor>
{
return from.template cast<To>();
}
LDLT_EXPLICIT_TPL_DECL(2, matmul_impl<long double>);
LDLT_EXPLICIT_TPL_DECL(1, mat_cast<proxqp::f64, long double>);
LDLT_EXPLICIT_TPL_DECL(1, mat_cast<proxqp::f32, long double>);
template<typename MatLhs, typename MatRhs, typename T = typename MatLhs::Scalar>
auto
matmul(MatLhs const& a, MatRhs const& b) -> Mat<T, proxqp::colmajor>
{
using Upscaled = typename std::
conditional<std::is_floating_point<T>::value, long double, T>::type;
return mat_cast<T, Upscaled>(matmul_impl<Upscaled>(
Mat<T, proxqp::colmajor>(a).template cast<Upscaled>(),
Mat<T, proxqp::colmajor>(b).template cast<Upscaled>()));
}
template<typename MatLhs,
typename MatMid,
typename MatRhs,
typename T = typename MatLhs::Scalar>
auto
matmul3(MatLhs const& a, MatMid const& b, MatRhs const& c)
-> Mat<T, proxqp::colmajor>
{
return matmul(matmul(a, b), c);
}
VEG_TAG(from_data, FromData);
struct EigenNoAlloc
{
EigenNoAlloc(EigenNoAlloc&&) = delete;
EigenNoAlloc(EigenNoAlloc const&) = delete;
auto operator=(EigenNoAlloc&&) -> EigenNoAlloc& = delete;
auto operator=(EigenNoAlloc const&) -> EigenNoAlloc& = delete;
#if defined(EIGEN_RUNTIME_NO_MALLOC)
EigenNoAlloc() noexcept { Eigen::internal::set_is_malloc_allowed(false); }
~EigenNoAlloc() noexcept { Eigen::internal::set_is_malloc_allowed(true); }
#else
EigenNoAlloc() = default;
#endif
};
template<typename Scalar>
proxsuite::proxqp::dense::Model<Scalar>
dense_unconstrained_qp(proxqp::isize dim,
Scalar sparsity_factor,
Scalar strong_convexity_factor = Scalar(1e-2))
{
Mat<Scalar, colmajor> H =
rand::sparse_positive_definite_rand_not_compressed<Scalar>(
dim, strong_convexity_factor, sparsity_factor);
Vec<Scalar> g = rand::vector_rand<Scalar>(dim);
Mat<Scalar, colmajor> A =
rand::sparse_matrix_rand_not_compressed<Scalar>(0, dim, sparsity_factor);
Vec<Scalar> b = rand::vector_rand<Scalar>(0);
Mat<Scalar, colmajor> C =
rand::sparse_matrix_rand_not_compressed<Scalar>(0, dim, sparsity_factor);
Vec<Scalar> u = rand::vector_rand<Scalar>(0);
Vec<Scalar> l = rand::vector_rand<Scalar>(0);
proxsuite::proxqp::dense::Model<Scalar> model(dim, 0, 0);
model.H = H;
model.g = g;
return model;
}
template<typename Scalar>
proxsuite::proxqp::dense::Model<Scalar>
dense_strongly_convex_qp(proxqp::isize dim,
proxqp::isize n_eq,
proxqp::isize n_in,
Scalar sparsity_factor,
Scalar strong_convexity_factor = Scalar(1e-2))
{
Mat<Scalar, colmajor> H =
rand::sparse_positive_definite_rand_not_compressed<Scalar>(
dim, strong_convexity_factor, sparsity_factor);
Vec<Scalar> g = rand::vector_rand<Scalar>(dim);
Mat<Scalar, colmajor> A =
rand::sparse_matrix_rand_not_compressed<Scalar>(n_eq, dim, sparsity_factor);
Mat<Scalar, colmajor> C =
rand::sparse_matrix_rand_not_compressed<Scalar>(n_in, dim, sparsity_factor);
Vec<Scalar> x_sol = rand::vector_rand<Scalar>(dim);
auto delta = Vec<Scalar>(n_in);
for (proxqp::isize i = 0; i < n_in; ++i) {
delta(i) = rand::uniform_rand();
}
Vec<Scalar> u = C * x_sol + delta;
Vec<Scalar> b = A * x_sol;
Vec<Scalar> l(n_in);
l.setZero();
l.array() -= 1.e20;
proxsuite::proxqp::dense::Model<Scalar> model(dim, n_eq, n_in);
model.H = H;
model.g = g;
model.A = A;
model.b = b;
model.C = C;
model.u = u;
model.l = l;
return model;
}
template<typename Scalar>
proxsuite::proxqp::dense::Model<Scalar>
dense_not_strongly_convex_qp(proxqp::isize dim,
proxqp::isize n_eq,
proxqp::isize n_in,
Scalar sparsity_factor)
{
Mat<Scalar, colmajor> H =
rand::sparse_positive_definite_rand_not_compressed<Scalar>(
dim, Scalar(0), sparsity_factor);
Mat<Scalar, colmajor> A =
rand::sparse_matrix_rand_not_compressed<Scalar>(n_eq, dim, sparsity_factor);
Mat<Scalar, colmajor> C =
rand::sparse_matrix_rand_not_compressed<Scalar>(n_in, dim, sparsity_factor);
Vec<Scalar> x_sol = rand::vector_rand<Scalar>(dim);
Vec<Scalar> y_sol = rand::vector_rand<Scalar>(n_eq);
Vec<Scalar> z_sol = rand::vector_rand<Scalar>(n_in);
auto delta = Vec<Scalar>(n_in);
for (proxqp::isize i = 0; i < n_in; ++i) {
delta(i) = rand::uniform_rand();
}
auto Cx = C * x_sol;
Vec<Scalar> u = Cx + delta;
Vec<Scalar> b = A * x_sol;
Vec<Scalar> l = Cx - delta;
Vec<Scalar> g = -(H * x_sol + C.transpose() * z_sol + A.transpose() * y_sol);
proxsuite::proxqp::dense::Model<Scalar> model(dim, n_eq, n_in);
model.H = H;
model.g = g;
model.A = A;
model.b = b;
model.C = C;
model.u = u;
model.l = l;
return model;
}
template<typename Scalar>
proxsuite::proxqp::dense::Model<Scalar>
dense_degenerate_qp(proxqp::isize dim,
proxqp::isize n_eq,
proxqp::isize n_in,
Scalar sparsity_factor,
Scalar strong_convexity_factor = Scalar(1e-2))
{
Mat<Scalar, colmajor> H =
rand::sparse_positive_definite_rand_not_compressed<Scalar>(
dim, strong_convexity_factor, sparsity_factor);
Vec<Scalar> g = rand::vector_rand<Scalar>(dim);
Mat<Scalar, colmajor> A =
rand::sparse_matrix_rand_not_compressed<Scalar>(n_eq, dim, sparsity_factor);
Mat<Scalar, colmajor> C = Mat<Scalar, proxqp::colmajor>(2 * n_in, dim);
Vec<Scalar> x_sol = rand::vector_rand<Scalar>(dim);
auto delta = Vec<Scalar>(2 * n_in);
for (proxqp::isize i = 0; i < 2 * n_in; ++i) {
delta(i) = rand::uniform_rand();
}
Vec<Scalar> b = A * x_sol;
auto C_ =
rand::sparse_matrix_rand_not_compressed<Scalar>(n_in, dim, sparsity_factor);
C.setZero();
C.block(0, 0, n_in, dim) = C_;
C.block(n_in, 0, n_in, dim) = C_;
Vec<Scalar> u = C * x_sol + delta;
Vec<Scalar> l(2 * n_in);
l.setZero();
l.array() -= 1.e20;
proxsuite::proxqp::dense::Model<Scalar> model(dim, n_eq, n_in);
model.H = H;
model.g = g;
model.A = A;
model.b = b;
model.C = C;
model.u = u;
model.l = l;
return model;
}
template<typename Scalar>
proxsuite::proxqp::dense::Model<Scalar>
dense_box_constrained_qp(proxqp::isize dim,
proxqp::isize n_eq,
proxqp::isize n_in,
Scalar sparsity_factor,
Scalar strong_convexity_factor = Scalar(1e-2))
{
Mat<Scalar, colmajor> H =
rand::sparse_positive_definite_rand_not_compressed<Scalar>(
dim, strong_convexity_factor, sparsity_factor);
Vec<Scalar> g = rand::vector_rand<Scalar>(dim);
Mat<Scalar, colmajor> A =
rand::sparse_matrix_rand_not_compressed<Scalar>(n_eq, dim, sparsity_factor);
Mat<Scalar, colmajor> C = Mat<Scalar, proxqp::colmajor>(n_in, dim);
Vec<Scalar> x_sol = rand::vector_rand<Scalar>(dim);
auto delta = Vec<Scalar>(n_in);
for (proxqp::isize i = 0; i < n_in; ++i) {
delta(i) = rand::uniform_rand();
}
Vec<Scalar> b = A * x_sol;
C.setZero();
C.diagonal().array() += 1;
Vec<Scalar> u = x_sol + delta;
Vec<Scalar> l = x_sol - delta;
proxsuite::proxqp::dense::Model<Scalar> model(dim, n_eq, n_in);
model.H = H;
model.g = g;
model.A = A;
model.b = b;
model.C = C;
model.u = u;
model.l = l;
return model;
}
template<typename Scalar>
proxsuite::proxqp::sparse::SparseModel<Scalar>
sparse_strongly_convex_qp(proxqp::isize dim,
proxqp::isize n_eq,
proxqp::isize n_in,
Scalar sparsity_factor,
Scalar strong_convexity_factor = Scalar(1e-2))
{
SparseMat<Scalar> H = rand::sparse_positive_definite_rand_compressed<Scalar>(
dim, strong_convexity_factor, sparsity_factor);
Vec<Scalar> g = rand::vector_rand<Scalar>(dim);
SparseMat<Scalar> A =
rand::sparse_matrix_rand<Scalar>(n_eq, dim, sparsity_factor);
SparseMat<Scalar> C =
rand::sparse_matrix_rand<Scalar>(n_in, dim, sparsity_factor);
Vec<Scalar> x_sol = rand::vector_rand<Scalar>(dim);
auto delta = Vec<Scalar>(n_in);
for (proxqp::isize i = 0; i < n_in; ++i) {
delta(i) = rand::uniform_rand();
}
Vec<Scalar> u = C * x_sol + delta;
Vec<Scalar> b = A * x_sol;
Vec<Scalar> l(n_in);
l.setZero();
l.array() -= 1.e20;
proxsuite::proxqp::sparse::SparseModel<Scalar> res{ H, g, A, b, C, u, l };
return res;
}
} // namespace utils
} // namespace proxqp
} // namespace proxsuite
#endif /* end of include guard PROXSUITE_PROXQP_UTILS_RANDOM_QP_PROBLEMS_HPP \
*/