2
$\begingroup$

I'm looking for an open-source library for the generation of random n-dimensional orthogonal matrices in C++.

In python, it looks like such a function is available in the NumPy package. But I was not able to find yet a solution in C++. Any suggestions?

$\endgroup$

3 Answers 3

3
$\begingroup$

I've seen in your comment that you want a uniform sampling.

With the Eigen library, you can uniformly generate at random a unit quaternion:

Eigen::Quaterniond q = Eigen::Quaterniond::UnitRandom();

and then convert it to a rotation (orthogonal) matrix:

Eigen::MatrixXd M = q.toRotationMatrix();
$\endgroup$
2
  • 1
    $\begingroup$ This works only under the very restrictive assumption that $n = 3$, right? $\endgroup$ Apr 27, 2020 at 18:57
  • $\begingroup$ @FedericoPoloni Ah yes, you're right, stupid of me :-/ $\endgroup$ Apr 27, 2020 at 19:00
2
$\begingroup$

The example you cited appears to be generating random Householder vectors and multiplying them out using backwards accumulation.

Another simple thing to do would be to generate a random matrix $\mathbf A$, then compute its $\mathbf A=\mathbf Q \mathbf R$ decomposition and discard the $\mathbf R$ factor. The two LAPACK functions that you need are [geqrf] (to factor $\mathbf A = \mathbf Q \mathbf R$ using Householder reflections) and [orgqr] (to reconstruct $\mathbf Q$ from the implicit reflector form). Although this is more work (about 2x) than the numpy algorithm, it might be faster in practice for large N because it will use BLAS3 kernels (whereas one-by-one backwards accumulation is only BLAS2).

There are a lot of C++ algebraic libraries, just look for one that happens to wrap these two functions. I happen to author/maintain such a library that has them (myramath), see below for a test program to generate a random orthogonal matrix:

#include <myramath/dense/Matrix.h>
#include <myramath/dense/geqrf.h>
#include <myramath/dense/orgqr.h>

#include <myramath/dense/gemm.h>
#include <myramath/dense/frobenius.h>

#include <iostream>

myra::Matrix<double> rvs(int N)
  {
  auto A = myra::Matrix<double>::random(N,N);
  auto tau = myra::geqrf_inplace(A);
  myra::orgqr_inplace(A,tau);
  return A;
  }

int main()
  {
  // Form Q.
  int N = 10;
  auto Q = rvs(N);
  std::cout << "Q = " << Q << std::endl;
  // Check Q is orthogonal.
  auto I = myra::Matrix<double>::identity(N);
  std::cout << "|Q'Q-I| = " << myra::frobenius(myra::gemm(Q,'T',Q)-I) << std::endl;
  std::cout << "|QQ'-I| = " << myra::frobenius(myra::gemm(Q,Q,'T')-I) << std::endl;
  return 0;
  }
$\endgroup$
4
  • 1
    $\begingroup$ Hmm, upon reading the underlying code in numpy, it's possible that these two approaches are not the same in the statistical sense. Their documentation indicates that numpy.rvs() draws uniformly from the space of orthogonal matrices, I am not sure that my suggested approach (draw uniform random A, then QR it) does the same. [Granted, I am not sure that it doesn't .. I just don't know enough statistics to state with certainty either way]. If you just need an orthogonal Q that is unpredictable, this suffices. But if you need some certain/stronger statistics, maybe stick with the numpy algorithm. $\endgroup$ Mar 28, 2020 at 1:29
  • $\begingroup$ Thanks for your suggestion. Indeed I need an implementation that samples a matrix uniformly at random from the space of orthogonal matrices. $\endgroup$ Mar 28, 2020 at 12:29
  • 1
    $\begingroup$ This paper has a really nice exposition of how to generate uniformly random unitary matrices. It does involve generating an arbitrary random matrix and then taking the QR factorization, but there's one extra step to guarantee uniformity. $\endgroup$ Mar 28, 2020 at 17:51
  • $\begingroup$ That paper of Mezzadri warns that the QR decomposition is not unique for unitary case and that many software packages produce Q and R with phases chosen such that resulting R is not uniformly sampled. Ha also suggests the cure for it. AFAIK these problems do not apply to the real orthogonal case. $\endgroup$
    – Kphysics
    Feb 11 at 0:24
1
$\begingroup$

If $X$ is a $(n \times m)$-matrix whose entries are independently generated values from the standard normal distribution, then $X(X^{\top}X)^{-\frac12}$ is a uniformly generated random orthogonal matrix. Source.

Here is an implementation with Eigen:

#include<iostream>
#include<random>
#include<Eigen/Eigen>
#include<ctime>
using namespace std;
using namespace Eigen;

static default_random_engine e(time(0));
static normal_distribution<double> gaussian(0,1);

MatrixXd randomOrthogonalMatrix(const unsigned long n){
  MatrixXd X = MatrixXd::Zero(n,n).unaryExpr([](double dummy){return gaussian(e);});
  MatrixXd XtX = X.transpose() * X;
  SelfAdjointEigenSolver<MatrixXd> es(XtX);
  MatrixXd S = es.operatorInverseSqrt();
  return X * S;
}

I'm not a star with C++, not sure where to place the random generator.

$\endgroup$
1
  • $\begingroup$ This works beautifully. Note that the routine may generate improper rotation matrices, i.e. those that are the product of a proper rotation and a reflection. Improper rotation matrices have a determinant of -1, "pure" rotations have +1. If you need pure rotations, simply pre-multiply the returned matrix by an identity matrix which has -1.0 instead of 1.0 as its upper left corner element. $\endgroup$ Sep 15, 2020 at 15:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.