Machine precision problem calculating random orthogonal matrix

Question

I was trying to calculate random orthogonal matrix of any size and faced a problem that machine error is huge since small sizes of matrix. I check final matrix is orthogonal by Q^T * Q = I Where Q is calculated orthogonal matrix. For example this operation for 10*10 matrix returns

1.000001421720586184 -0.000000728640227713 0.000001136830463799 -0.000000551609342727 -0.000001177027039965 0.000000334599582398 -0.000000858589995413 0.000000954985769303 0.000032744809653293 -0.000000265053286108 
-0.000000728640227713 1.000000373167701527 -0.000000583888104495 0.000000285028920487 0.000000602479963850 -0.000000171504851561 0.000000439149502041 -0.000000489282575621 -0.000016836862737655 0.000000139458235281 
0.000001136830463799 -0.000000583888104495 1.000000903071175539 -0.000000430043020645 -0.000000944743177025 0.000000267453533747 -0.000000690730534104 0.000000764348989692 0.000025922602192780 -0.000000194784469538 
-0.000000551609342727 0.000000285028920487 -0.000000430043020645 1.000000193583126484 0.000000463290242271 -0.000000129639515781 0.000000340879278672 -0.000000371868992193 -0.000012221761904027 0.000000071060844115 
-0.000001177027039965 0.000000602479963850 -0.000000944743177025 0.000000463290242271 1.000000972303993511 -0.000000277069934225 0.000000708304641621 -0.000000790186119274 -0.000027265457679301 0.000000229727452845 
0.000000334599582398 -0.000000171504851561 0.000000267453533747 -0.000000129639515781 -0.000000277069934225 1.000000078745856667 -0.000000202136378957 0.000000224766235153 0.000007702164180343 -0.000000062098652538 
-0.000000858589995413 0.000000439149502041 -0.000000690730534104 0.000000340879278672 0.000000708304641621 -0.000000202136378957 1.000000515565399593 -0.000000576213323968 -0.000019958173806416 0.000000172131276688 
0.000000954985769303 -0.000000489282575621 0.000000764348989692 -0.000000371868992193 -0.000000790186119274 0.000000224766235153 -0.000000576213323968 1.000000641385961051 0.000022026878760393 -0.000000180133590665 
0.000032744809653293 -0.000016836862737655 0.000025922602192780 -0.000012221761904027 -0.000027265457679301 0.000007702164180343 -0.000019958173806416 0.000022026878760393 1.000742765170839780 -0.000005353869978161 
-0.000000265053286108 0.000000139458235281 -0.000000194784469538 0.000000071060844115 0.000000229727452845 -0.000000062098652538 0.000000172131276688 -0.000000180133590665 -0.000005353869978161 1.000000000000000000

So we can see that matrix is orthogonal but non-diagonal elements have big error Is there any solution of that?

How I calculate n*n orthogonal matrix:

1. Take start square orthogonal matrix of size 1 Q = |1|
2. Take random vector of dimension 2 y = |rand(), rand()| and norm it y = y/norm(y)
3. Construct a Householder reflection with vector y and apply it to matrix Q with number 1 in right corner, so I have orthogonal matrix Q of size 2.
4. Repeat while don't have n*n matrix, taking new random y with increased dimension.

Code:

#include <iostream>
#include <map>
#include <iostream>
#include <vector>
#include <iterator>
#include <cmath>
#include <utility>
using namespace std;
template<typename T>
T tolerance = T(1e-3);

template<typename T>
struct Triplet{
    int i;
    int j;
    T b;
};
template<typename T>
T Tabs(T num){
    if(num<T(0)) return -num;
    else return num;
}

template<typename T>
class DOK{
private:
    /*
     * Dictionary of Keys, pair<int, int> is coordinates of non-zero elements,
     * next int is value
     */

    int size_n;
    int size_m;
    map<pair<int, int>, T> dict;
    // int count;
public:

    DOK(vector<Triplet<T>> &matrix, int n, int m){
        this->resize(n, m);
        this->fill(matrix);
    }

    DOK(int n, int m){
        this->resize(n, m);
    }
    ~DOK() = default;

    void fill(vector<Triplet<T>> &matrix){
        //this->count=matrix.size();
        //cout<<"Input your coordinates with value in format \"i j val\" "<<endl;
        for(int k = 0; k < matrix.size(); k++){
            this->insert(matrix[k]);
        }
    }

    void insert(const Triplet<T> &Element){
        if(Element.i >= this->size_n){
            this->size_n = Element.i+1;
        }
        if(Element.j >= this->size_m){
            this->size_m = Element.j+1;
        }
        pair<int, int> coordinates = {Element.i, Element.j};
        this->dict.insert(pair(coordinates, Element.b));
    }

    void resize(int n, int m){
        this->size_n=n;
        this->size_m=m;
    }
    void print() const{
        cout<<endl;
        for(int i = 0; i < this->size_n; i++){
            for(int j = 0; j < this->size_m; j++){
                if(this->dict.find({i, j})!= this->dict.cend()) cout<< this->dict.find(pair(i, j))->second<<" "; else cout<<0<<" ";
            }
            cout<<endl;
        }
    }

    void clearZeros(){
        for(auto i = this->dict.begin(); i!=this->dict.end();){
            if(Tabs(i->second) <=  tolerance<T>){
                i = this->dict.erase(i);
            } else{
                i++;
            }
        }
    }

    [[nodiscard]] pair<int, int> getSize() const{
        return {size_n, size_m};
    }


    DOK<T> transpose(){
        DOK<T> A = DOK<T>(this->size_m, this->size_n);
        for(auto &i: this->dict){
            A.insert({i.first.second, i.first.first, i.second});
        }
        return A;
    }

    DOK<T>& operator-=(const DOK<T>& matrix){
        try{
            if(this->size_n != matrix.size_n || this->size_m != matrix.size_m) throw 1;
            for(auto j: matrix.dict){
                if(this->dict.find(j.first)!=this->dict.cend()) {
                    this->dict[j.first] -= j.second;
                }else{
                    this->dict.insert({j.first, -j.second});
                    //M.count++;
                }
            }
            this->clearZeros();
            return *this;
        }
        catch (int a) {
            cout<<"Sizes of Matrices are different."<<endl;
        }
    }

    DOK<T> operator-(const DOK<T> &matrix) const{
        DOK<T> t = *this;
        return move(t-=matrix);
    }

    DOK<T>& operator*=(const DOK<T> &matrix){
        try {
            if(this->size_m != matrix.size_n) throw 1;
            DOK<T> M = DOK(this->size_n, matrix.size_m);
            for (int i = 0; i < this->size_n; i++) {
                for (int j = 0; j < matrix.size_m; j++) {
                    T a=0;
                    for(int k = 0; k<this->size_m; k++){
                        if(this->dict.find({i,k}) != this->dict.cend() && matrix.dict.find({k, j})!=matrix.dict.cend()){
                            a+=this->dict.find({i,k})->second*matrix.dict.find({k,j})->second;
                            //cout<<a<<endl;
                        }
                    }
                    Triplet<T> m = {i, j, a};
                    M.insert(m);
                }
            }
            this->clearZeros();
            *this=M;
            return *this;
        }
        catch (int a) {
            cout<<"Wrong sizes of matrices to multiplication"<<endl;
        }
    }

    DOK<T> operator*(const DOK<T>& matrix) const{
        DOK<T> t = *this;
        return t*=matrix;
    }

    DOK<T>& operator*=(T& k){
        for(auto i: this->dict){
            this->dict[i.first]*=k;
        }
        this->clearZeros();
        return *this;
    }

    DOK<T> operator*(T& k) const{
        DOK<T> t = *this;
        return move(t*=k);
    }
    DOK<T>& operator*=(const T& k){
        for(auto i: this->dict){
            this->dict[i.first]*=k;
        }
        this->clearZeros();
        return *this;
    }


};

template<typename T>
vector<T> operator*(const DOK<T> &matrix, const vector<T> &x){
    vector<T> result;
    for(int i = 0; i < x.size(); i++){
        T temp = 0;
        for(int j = 0; j < x.size(); j++){
            temp+=matrix(i, j)*x[j];
        }
        result.push_back(temp);
    }
    return move(result);
}

template<typename T>
T operator*(const vector<T> &x, const vector<T> &b) {
    T result = 0;
    for(int i = 0; i < x.size(); i++){
        result+=x[i]*b[i];
    }
}

template<typename T>
vector<T> operator*(const vector<T> &x, const DOK<T> &matrix) {
    vector<T> result;
    for(int i = 0; i < x.size(); i++){
        T temp = 0;
        for(int j = 0; j < x.size(); j++){
            temp+=matrix(j, i)*x[j];
        }
        result.push_back(temp);
    }
    return move(result);
}

template<typename  T>
DOK<T> operator*(T& k, const DOK<T> &matrix){
    return matrix*k;
}
template<typename  T>
DOK<T> operator*(const T& k, const DOK<T> &matrix){
    return matrix*k;
}

template<typename T>
vector<T>& operator*=(const DOK<T> &matrix, const vector<T> &x){
    return matrix*x;
}

template<typename T>
vector<T>& operator*=(const vector<T> &x, const DOK<T> &matrix){
    return x*matrix;
}

template<typename T>
vector<T> operator*(const vector<T> &x, T k){
    vector<T> result = x;
    for(int i = 0; i<x.size(); i++){
        result[i]*=k;
    }
    return result;
}
template<typename T>
vector<T> operator*(T k, const vector<T> &x) {
    return x*k;
}
template<typename  T>
ostream& operator<<(ostream &os, const DOK<T> &matrix) {
    matrix.DOK<T>::print();
    return os;
}

template<typename T>
T norm(const vector<T> x){
    T result = 0;
    for(int i = 0; i < x.size(); i++){
        result+=pow(x[i],2);
    }
    return sqrt(result);
}

template<typename T>
DOK<T> Ortogonal(int n){
    srand(time(NULL));
    vector<Triplet<T>> in = {{0, 0, T(1)}};
    DOK<T> Q = DOK<T>(in, 1, 1);
    DOK<T> E = Q;
    vector<T> y;

    for(int i = 1; i<n; i++){
        y.clear();
        for(int m = 0; m<i+1; m++){
            y.emplace_back(rand());
        }

        y = (1/norm(y))*y;
        DOK<T> Y = DOK<T>(i+1, i+1);
        for(int j = 0; j<i+1; j++){
            for(int k = 0; k<i+1; k++){
                Y.insert({j, k, y[j]*y[k]});
            }
        }
        Q.insert({i, i, T(1)});
        cout<<Q;
        Y*=T(2);
        E.insert({i, i, T(1)});
        Q = (E - Y)*Q;
    }

    return Q;
}

main.cpp:

#include <iostream>
#include "DOK.h"
using namespace std;
int main() {
    DOK<long double> O = Ortogonal<long double>(10);

    cout<<O.transpose()*O;
    return 0;
}

DOK is template class of sparse matrix with all overloaded operators.

Answer 1

I have tried to reproduce the problem, and it seems that, with the Householder-based procedure, the amount of numerical error changes for each execution of the program.

This means that for some reason the error is sensitive to the random number sequence used. That sequence depends on the exact launching time of the program thru the srand(time(NULL)) function call.

If the Householder procedure is not an essential requirement, I would suggest using a Gram-Schmidt based procedure instead.

The procedure works like this:

1) create a full n*n random square matrix
2) normalize its column vectors
3) for each column vector except the leftmost one:
       substract from that vector its projections on column vectors located to its left
       renormalize the column vector

Unfortunately, I do not feel conversant enough with the DOK sparse matrix framework to work out such a solution within it. But it could certainly be done by a programmer more familiar with it.

Using the Eigen open source linear algebra library instead, the code looks like this:

#include  <vector>
#include  <random>
#include  <iostream>
#include  <Eigen/Dense>

using namespace  Eigen;

MatrixXd  makeRandomOrthogonalMatrix(int dim, int seed)
{
    std::normal_distribution<double>  dist(0.0, 1.0);
    std::mt19937                      gen(seed);
    std::function<double(void)>  rng = [&dist, &gen] () { return dist(gen); };

    MatrixXd rawMat = MatrixXd::NullaryExpr(dim, dim, rng);

    // massage rawMat into an orthogonal matrix - Gram-Schmidt process:
    for (int j=0; j < dim; j++) {
        VectorXd colj = rawMat.col(j);
        colj.normalize();
        rawMat.col(j) = colj;
    }
    for (int j=1; j < dim; j++) {
        VectorXd colj = rawMat.col(j);
        for (int k = 0; k < j; k++) {
            VectorXd  colk = rawMat.col(k);
            colj -= (colk.dot(colj)) * colk;
        }
        colj.normalize();
        rawMat.col(j) = colj;
    }
    return  rawMat;
}

Side note on the production of pseudo-random numbers:

Function rand() in the DOK program provides an uniform distribution between zero and a large integer value. Here, I switched to a Gaussian (a.k.a. “normal”) distribution because:

1. this is what the relevant mathematical literature is using
2. this ensures the random column vectors have no privileged directions

Rotational symmetry comes simply from the fact that expression exp(-x²) * exp(-y²) * exp(-z²) is equal to exp(-(x² + y² + z²)).

Test program used:

int main()
{
    int           ranSeed = 4241;
    const int     dim     = 10;

    MatrixXd  orthoMat = makeRandomOrthogonalMatrix(dim, ranSeed);

    MatrixXd  trOrthoMat  = orthoMat.transpose();
    MatrixXd  productMat2 = trOrthoMat * orthoMat;

    std::cout << "productMat2 = \n" << productMat2 << '\n' << std::endl;

    return EXIT_SUCCESS;
}

The amount of numerical error looks satisfactory, even after trying several values for the random seed:

productMat2 = 
          1  1.38778e-17            0  1.38778e-17 -9.71445e-17  1.79544e-16  1.11022e-16  1.80411e-16  4.16334e-17 -2.98372e-16
 1.38778e-17            1            0 -6.93889e-17 -1.38778e-17  5.89806e-17  5.55112e-17 -4.16334e-17 -2.77556e-17 -7.63278e-17
           0            0            1            0 -1.38778e-17  7.37257e-17            0  9.71445e-17  1.80411e-16 -1.38778e-17
 1.38778e-17 -6.93889e-17            0            1  1.38778e-17 -1.64799e-17 -5.55112e-17            0 -5.55112e-17  2.77556e-17
-9.71445e-17 -1.38778e-17 -1.38778e-17  1.38778e-17            1 -4.25007e-17  5.55112e-17 -1.38778e-17 -6.93889e-17  1.78677e-16
 1.79544e-16  5.89806e-17  7.37257e-17 -1.64799e-17 -4.25007e-17            1  7.80626e-17 -5.37764e-17  2.60209e-17 -9.88792e-17
 1.11022e-16  5.55112e-17            0 -5.55112e-17  5.55112e-17  7.80626e-17            1  8.32667e-17  8.32667e-17 -7.63278e-17
 1.80411e-16 -4.16334e-17  9.71445e-17            0 -1.38778e-17 -5.37764e-17  8.32667e-17            1  4.16334e-17  1.04083e-16
 4.16334e-17 -2.77556e-17  1.80411e-16 -5.55112e-17 -6.93889e-17  2.60209e-17  8.32667e-17  4.16334e-17            1 -7.28584e-17
-2.98372e-16 -7.63278e-17 -1.38778e-17  2.77556e-17  1.78677e-16 -9.88792e-17 -7.63278e-17  1.04083e-16 -7.28584e-17            1

Using the function in a non-Eigen context:

To allow a DOK-based program to use that random unitary matrix, we can use a regular std::vector object in order to shield the DOK type system from the Eigen type system. With this plumbing code:

// for export to non-Eigen programs:
std::vector<double>  makeRandomOrthogonalMatrixAsVector(int dim, int seed)
{
    MatrixXd randMat = makeRandomOrthogonalMatrix(dim, seed);

    std::vector<double>  matVec;
    for (int i=0; i < dim; i++)
        for (int j=0; j < dim; j++)
            matVec.push_back(randMat(i,j));
    return matVec;
}

As a last step, the random unitary matrix can be rebuilt on the DOK side of things, without any undue exposure to Eigen types. This code can be used to rebuild the matrix from the intermediate vector:

extern vector<double>  makeRandomOrthogonalMatrixAsVector(int, int);

int main(int argc, const char* argv[])
{
    int         dim = 10;
    int  randomSeed = 4241;

    vector<double> vec = makeRandomOrthogonalMatrixAsVector(dim, randomSeed);

    // create matrix in DoK format:
    vector<Triplet<double>>  trv;
    for (int i=0; i<dim; i++)
        for (int j=0; j<dim; j++) {
            double x = vec[i*dim + j];
            trv.push_back(Triplet<double>{i, j, x});
        }
    DOK<double>  orthoMat(trv, dim, dim);

    DOK<double>  prodMat = orthoMat.transpose() * orthoMat;

    cout << "\nAlmost Unit Matrix = \n" << prodMat << std::endl;

    return 0;
}

DOK program output:

Almost Unit Matrix = 

1 0 4.16334e-17 1.11022e-16 -8.32667e-17 -2.77556e-17 0 -2.77556e-17 -1.38778e-17 8.32667e-17 
0 1 -1.5786e-16 -1.38778e-16 -5.55112e-17 2.42861e-17 1.11022e-16 2.22045e-16 -5.96745e-16 -2.77556e-17 
4.16334e-17 -1.5786e-16 1 1.50487e-16 6.93889e-18 5.55112e-17 -5.72459e-17 -8.32667e-17 5.10009e-16 -1.83881e-16 
1.11022e-16 -1.38778e-16 1.50487e-16 1 -6.93889e-18 -1.21431e-17 1.73472e-17 0 1.45717e-16 2.94903e-17 
-8.32667e-17 -5.55112e-17 6.93889e-18 -6.93889e-18 1 -9.36751e-17 -2.77556e-16 -2.498e-16 5.68989e-16 -2.91434e-16 
-2.77556e-17 2.42861e-17 5.55112e-17 -1.21431e-17 -9.36751e-17 1 1.04083e-16 -3.46945e-17 1.04083e-17 -3.40006e-16 
0 1.11022e-16 -5.72459e-17 1.73472e-17 -2.77556e-16 1.04083e-16 1 -1.38778e-16 2.77556e-16 -1.04083e-16 
-2.77556e-17 2.22045e-16 -8.32667e-17 0 -2.498e-16 -3.46945e-17 -1.38778e-16 1 0 9.71445e-17 
-1.38778e-17 -5.96745e-16 5.10009e-16 1.45717e-16 5.68989e-16 1.04083e-17 2.77556e-16 0 1 2.08167e-17 
8.32667e-17 -2.77556e-17 -1.83881e-16 2.94903e-17 -2.91434e-16 -3.40006e-16 -1.04083e-16 9.71445e-17 2.08167e-17 1 

So, as expected, the numerical errors remain tiny.

Note that unlike for the DOK program in your posting, the random number sequence used always remains the same, unless the value of the random seed is changed manually within the source code. A common improvement consists in passing the random seed as a command line argument, using some code like: int seed = std::stoi(argv[1]);.