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

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

转载作者：行者123 更新时间：2023-12-03 07:03:21

我试图计算任意大小的随机正交矩阵，并遇到一个问题，因为矩阵尺寸较小，因此机器误差很大。我通过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++ - 为已经转换的QGraphicsItem设置转换原点

文章推荐： c++ - 如何自动更新最大值？

文章推荐： c++ - 如何用C++中的两个字符串 vector 创建一个字符串？

服务器端的 Firebird 计算(计算)字段
SQL 和一般开发的新手，我有一个表(COUNTRIES)，其中包含字段(INDEX、NAME、POPULATION、AREA) 通常我添加一个客户端(Delphi)计算字段(DENSITY)和 On
jquery - 计算(百分比)计算(像素)
我想使用 calc(100%-100px)，但在我的 demo 中不起作用由于高度只接受像素，因此如何将此百分比值转换为像素。最佳答案以下将为您提供高度: $(window).height();
MySql 计算
我正在尝试在 MySQL 中添加列并动态填充其他列。例如我有一张表“数字”并具有第 1 列、第 2 列、第 3 列，这些总数应填充在第 4 列中最佳答案除非我误解了你的问题，否则你不只是在寻找:
mysql - 计算
我想返回简单计算的结果，但我不确定如何执行此操作。我的表格如下: SELECT COUNT(fb.engineer_id) AS `total_feedback`, SUM(fb.ra
嵌套for循环中的c++计算
我一直在尝试做这个程序，但我被卡住了，我仍然是一个初学者，任何帮助将不胜感激。我需要程序来做打印一个 10 X 10 的表格，其中表格中的每个条目都是行号和列号的总和包含一个累加器，用于计算所有表
c - 计算
这个计算背后一定有一些逻辑。但我无法得到它。普通数学不会导致这种行为。谁能帮我解释一下原因 printf ("float %f\n", 2/7 * 100.0); 结果打印 1.000000 为什么会
计算 AND 的算法
我想计算从 0 到 (n)^{1/2} - 1 的数字的 AND每个数字从 0 到 (n)^{1/2} - 1 .我想在 O(n) 中执行此操作时间，不能使用 XOR、OR、AND 运算。具体来说，
Excel - 在数字格式中使用公式/计算
如何在 Excel 中将公式放入自定义数字格式？例如(出于说明目的随机示例)，假设我有以下数据: 输入输出在不编辑单元格中的实际数据的情况下，我想显示单元格中的值除以 2，并保留两位小数: 有没
Flutter:隔离内存泄漏(计算)
每次我在 Flutter 应用程序中调用计算()时，我都会看到内存泄漏，据我所知，这基本上只是一种生成隔离的便捷方法。我的应用程序内存占用增加并且在 GC 之后永远不会减少。我已将我的代码简化为仅调
R中的RMSE(均方根偏差)计算
我有数字特征观察 V1通过 V12用于目标变量 Wavelength .我想计算 Vx 之间的 RMSE列。数据格式如下。每个变量“Vx”以 5 分钟的间隔进行测量。我想计算所有 Vx 变量的观测值
计算 C 文件中未知数量的字符
我正在寻找一种使用 C 语言计算文件中未知字符数的简单方法。谢谢你的帮助最佳答案 POSIX 方式(可能是您想要的方式): off_t get_file_length( FILE *file ) {
sql - 计算/派生连续日期跨度中的第一个开始日期
我正在使用 Postgres，并且我正试图围绕如何在连续日期跨度中得出第一个开始日期的问题进行思考。例如 :- ID | Start Date | End Date =================
jquery - 计算，用逗号替换点
我有一个订单表格，我在其中使用 jQuery 计算插件来汇总总数。此求和工作正常，但生成的“总和”存在问题。总之，我希望用逗号替换任何点。代码的基础是； function ($this) {
Delphi错误的 double 计算
我在使用 double 变量计算简单算术方程时遇到问题。我有一个具有 double 属性 Value 的组件，我将此属性设置为 100。然后我做一个简单的减法来检查这个值是否真的是 100: va
openssl CRC32 计算
我在这里看到了一些关于 CRC 32 计算的其他问题。但没有一个让我满意，因此是这样。 openssl 库是否有任何用于计算 CRC32 的 api 支持？我已经在为 SHA1 使用 openssl，
php - 计算-1个月时的PHP天错误
当我在PHP日期计算中遇到问题时，我感到惊讶。 $add = '- 30 days'; echo date('Y-m-01', strtotime($add)); // result is 2017-
持有变量的 JavaScript 计算
我正在使用 javascript 进行练习，我编写了这个脚本来计算 2 个变量的总和，然后在第三个方程中使用这个总和!关于如何完成这项工作的任何想法都将非常有用! First Number:
audio - sample 计算
我有一个来自EAC的提示单和一个包含完整专辑的FLAC文件。我正在尝试制作一些python脚本来播放文件，因为我需要能够设置在flac文件中开始的位置。如何从CueSheet格式MM:SS:FF转
javascript - 计算 for 循环中输入值的总和
这个问题已经有答案了: Adding two numbers concatenates them instead of calculating the sum (24 个回答) 已关闭去年。我有一个
使用输入和跨度字段的 Javascript 计算
4000 我需要上面字段 name="quantity" 和 id="price" 中的值，并使用 javascript 函数进行计算，并将其显示在字段 id= 中仅当我单击计算按钮时才显示“总

行者123

个人简介

我是一名优秀的程序员,十分优秀！

作者热门文章

滴滴打车优惠券免费领取

全站热门文章

首页

博学

6Ren·AI

商城

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