# Generating Random Orthogonal Matrices in C++

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?

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

• 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. – rchilton1980 Mar 28 '20 at 1:29
• Thanks for your suggestion. Indeed I need an implementation that samples a matrix uniformly at random from the space of orthogonal matrices. – Mateus de Oliveira Mar 28 '20 at 12:29
• 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. – Daniel Shapero Mar 28 '20 at 17:51

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

• This works only under the very restrictive assumption that $n = 3$, right? – Federico Poloni Apr 27 '20 at 18:57
• @FedericoPoloni Ah yes, you're right, stupid of me :-/ – Stéphane Laurent Apr 27 '20 at 19:00

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;