Generating Random Orthogonal Matrices in C++

Question

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?

Answer 1

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();

Answer 2

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;
  }

Answer 3

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.