c++ - 机器精度问题计算随机正交矩阵

标签 c++ matrix precision numerical-methods

我试图计算任意大小的随机正交矩阵，并遇到一个问题，因为矩阵尺寸较小，因此机器误差很大。我通过Q^T * Q = I检查最终矩阵是否正交
其中Q是正交矩阵的计算值。例如10 * 10矩阵的此操作返回

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
因此我们可以看到矩阵是正交的，但非对角元素有很大的误差
有什么解决办法吗？
我如何计算n * n正交矩阵:
  • 以大小为1的开始平方正交矩阵Q = |1|
  • 取维度为2的随机 vector y = |rand(), rand()|并对其进行规范y = y/norm(y)
  • 用 vector y构造一个Householder reflection并将其应用于右上角编号为1的矩阵Q，因此我得到大小为2的正交矩阵Q。
  • 在没有n * n矩阵的情况下重复执行，以增加维度的方式获取新的随机y。

    • 码:
    #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是具有所有重载运算符的稀疏矩阵的模板类。

    最佳答案

    我试图重现该问题，并且似乎使用基于Householder的过程，每次执行该程序时，数字错误量都会发生变化。
    这意味着由于某种原因，该错误对所使用的随机数序列很敏感。该顺序取决于通过srand(time(NULL))函数调用的程序的确切启动时间。
    如果不是Householder程序是必不可少的条件，我建议改用Gram-Schmidt based procedure。
    该过程如下所示:

    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
    不幸的是，我对DOK稀疏矩阵框架不够熟悉，无法在其中找到这样的解决方案。但这当然可以由更熟悉它的程序员来完成。
    改用Eigen开源线性代数库，代码如下所示:
    #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;
}
    关于产生伪随机数的旁注:
    DOK程序中的函数rand()提供了零和大整数之间的均匀分布。在这里，我切换到高斯(也称为“正态”)分布是因为:
  • 这是相关数学文献正在使用的
  • 这样可以确保随机列 vector 没有特权方向

    • 旋转对称性仅来自表达式exp(-x2)* exp(-y2)* exp(-z2)等于exp(-(x2 + y2 + z2))的事实。
    使用的测试程序:
    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;
}
    即使为随机种子尝试了多个值，数值误差的数量也看起来令人满意:
    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
    在非本征上下文中使用该函数:
    为了允许基于DOK的程序使用该随机unit矩阵，我们可以使用常规的std::vector对象，以使DOK类型的系统与本征类型的系统隔离。使用以下管道代码:
    // 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;
}
    最后一步，可以在事物的DOK方面重建随机unit矩阵，而不会过度暴露于特征类型。此代码可用于从中间 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程序输出:
    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
    因此，正如预期的那样，数值误差仍然很小。
    请注意，与发布中的DOK程序不同，除非在源代码中手动更改了随机种子的值，否则所使用的随机数序列始终保持不变。常见的改进包括使用诸如int seed = std::stoi(argv[1]);之类的代码将随机种子作为命令行参数传递。

    关于c++ - 机器精度问题计算随机正交矩阵，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/64804418/

    上一篇：c++ - 如何自动更新最大值？

    下一篇：c++ - 为已经转换的QGraphicsItem设置转换原点

    相关文章：

    c++ - Linux 和 Windows 上的输出差异

    c++ - 如何使用这个二叉树类？

    c++ - 在世界坐标中计算 z_far 平面的顶点。

    c++ - 使用迭代方法计算矩阵的行列式

    c - 获取二维矩阵的相邻元素(仅限深度一)

    ios - 为什么我在 GPU 上有很大的精度错误？

    c++ - 这段代码是字节序安全的吗？

    python - 如何在 Keras 中计算精度和召回率

    math - float 学坏了吗？

    c++ - 如何调用所有基类的析构函数？ (或一个共同的功能)

    ©2024 IT工具网  联系我们