randomMatrices.h

Go to the documentation of this file.

/* Copyright 2020 CNRS-AIST JRL */
 
#pragma once
 
#include <Eigen/Core>
 
#include <3rd-party/effolkronium/random.hpp>
 
namespace jrl::qp::test
{
template<typename Scalar>
struct scalar_normal_random_op
{
  scalar_normal_random_op(double mean = 0, double stddev = 1) : d_(mean, stddev) {}
 
  inline const Scalar operator()() const
  {
    return effolkronium::random_static::get(d_);
  }
 
  mutable std::normal_distribution<Scalar> d_;
};
 
inline auto randnVec(Eigen::Index size, double mean = 0, double stddev = 1)
{
  return Eigen::VectorXd::NullaryExpr(size, scalar_normal_random_op<double>(mean, stddev));
}
 
inline auto randnMat(Eigen::Index rows, Eigen::Index cols, double mean = 0, double stddev = 1)
{
  return Eigen::MatrixXd::NullaryExpr(rows, cols, scalar_normal_random_op<double>(mean, stddev));
}
 
inline Eigen::VectorXd randUnitVec(Eigen::Index size)
{
  Eigen::VectorXd r = randnVec(size);
  r.normalize();
  return r;
}
 
inline Eigen::MatrixXd randOrtho(Eigen::Index size, bool special = false)
{
  assert(size >= 0);
  Eigen::VectorXd buffv(size);
  Eigen::VectorXd buffw(size);
  Eigen::MatrixXd Q = Eigen::MatrixXd::Identity(size, size);
 
  if(size == 0) return Q;
 
  bool positiveDet;
  if(effolkronium::random_static::get<bool>())
  {
    Q(size - 1, size - 1) = 1;
    positiveDet = true;
  }
  else
  {
    Q(size - 1, size - 1) = -1;
    positiveDet = false;
  }
 
  for(Eigen::Index i = 2; i <= size; ++i)
  {
    auto v = buffv.head(i);
    auto w = buffw.head(i);
 
    // Random unit vector
    v = randUnitVec(i);
 
    // compute (v+sgn(v_1)*e_1)/||v+sgn(v_1)*e_1||
    bool positive = v[0] >= 0;
    if(positive)
    {
      v[0] += 1;
      if(i % 2 == 0) positiveDet = !positiveDet;
    }
    else
    {
      v[0] -= 1;
      positiveDet = !positiveDet;
    }
    v.normalize();
 
    Eigen::MatrixXd H;
    if(positive)
      H = -(Eigen::MatrixXd::Identity(i, i) - 2 * v * v.transpose());
    else
      H = Eigen::MatrixXd::Identity(i, i) - 2 * v * v.transpose();
    // Apply the Householder transformation H = -sgn(v_1)*(I-2*vv^T) : Q = H*Q
    w = Q.bottomRightCorner(i, i).transpose() * v;
    if(positive)
    {
      Q(size - i, size - i) = -1;
      Q.bottomRightCorner(i - 1, i - 1) *= -1;
      Q.bottomRightCorner(i, i).noalias() += 2 * v * w.transpose();
    }
    else
    {
      Q.bottomRightCorner(i, i).noalias() -= 2 * v * w.transpose();
    }
  }
 
  if(special && !positiveDet) Q.col(0) *= -1;
 
  return Q;
}
 
inline Eigen::MatrixXd randn(Eigen::Index rows, Eigen::Index cols, Eigen::Index rank = -1)
{
  assert(rank <= rows && rank <= cols && "Invalid rank");
  Eigen::Index p = std::min(rows, cols);
  if(rank < 0) rank = p;
 
  if(rank == 0) return Eigen::MatrixXd::Zero(rows, cols);
 
  if(rank == p) return randnMat(rows, cols);
 
  Eigen::VectorXd s = Eigen::VectorXd::Zero(p);
  s.head(rank).setRandom();
  // set the correct variance
  s.head(rank) *= std::sqrt((3. * static_cast<double>(rows) * static_cast<double>(cols)) / static_cast<double>(rank));
 
  Eigen::MatrixXd M(rows, cols);
 
  if(rows <= cols)
  {
    M.leftCols(rows).noalias() = randOrtho(rows) * s.asDiagonal();
    M.rightCols(cols - rows).setZero();
    return M * randOrtho(cols);
  }
  else
  {
    M.topRows(cols).noalias() = s.asDiagonal() * randOrtho(cols);
    M.bottomRows(rows - cols).setZero();
    return randOrtho(rows) * M;
  }
}
 
inline std::pair<Eigen::MatrixXd, Eigen::MatrixXd> randDependent(Eigen::Index cols,
                                                                 Eigen::Index rowsA,
                                                                 Eigen::Index rankA,
                                                                 Eigen::Index rowsB,
                                                                 Eigen::Index rankB,
                                                                 Eigen::Index rankAB)
{
  assert(rankA <= rowsA && rankA <= cols);
  assert(rankB <= rowsB && rankB <= cols);
  assert(rankAB >= rankA && rankAB >= rankB && rankAB <= rankA + rankB && rankAB <= cols);
 
  Eigen::MatrixXd M = randn(rankA + rankB, cols, rankAB);
  Eigen::MatrixXd A(rowsA, cols);
  Eigen::MatrixXd B(rowsB, cols);
 
  if(rankA == rowsA)
    A = M.topRows(rankA);
  else
    A.noalias() = randOrtho(rowsA).leftCols(rankA) * M.topRows(rankA);
 
  if(rankB == rowsB)
    B = M.bottomRows(rankB);
  else
    B.noalias() = randOrtho(rowsB).leftCols(rankB) * M.bottomRows(rankB);
 
  return {A, B};
}
} // namespace jrl::qp::test