代码之家  ›  专栏  ›  技术社区  ›  fc67

Gram Schmidt可能失去精度

  •  3
  • fc67  · 技术社区  · 10 年前

    我有一个在循环中使用Gram Schmidt的代码。我希望尽可能减少对该算法的调用次数,但问题是,尽管在调用前后得到了相同的结果,但当我使用这些值打印某些操作的结果时,它们是不同的。例如,在下面的代码中 abs(muGS[k][0]) - abs(before2) 应为0或非常接近0,因为此变量的打印值(调用前后)相同。然而,事实并非如此。 muGS 是一个双矩阵,其值通常在0和1之间。

    int k = 1;
    double before2;
    
    while(k < end) {
    
        before2 = muGS[k][0];
    
        gramSchmidt(b, muGS, cGS, k);
    
        //prints for debug
        if (abs(muGS[k][0]) - abs(before2) > 0.1) {
    
            if (abs(muGS[k][0]) - abs(before2) > 0.1)  {
                cout << "1 muGS[k] diff:" << abs(muGS[k][0]) - abs(before2) << endl;
                cout << "1 muGS[k] before:" << muGS[k][0] << endl;
                cout << "1 muGS[k] after:" << muGS[k][0] << endl;
                cout << "1 muGS[k] mult before:" << before2 * before2 << endl;
                cout << "1 muGS[k] mult after:" << muGS[k][0] * muGS[k][0] << endl;
                cout << "1 muGS[k] abs before:" << abs(before2) << endl;
                cout << "1 muGS[k] abs after:" << abs(muGS[k][0]) << endl;
            }
            getchar();
        }
    
        for (i = k-1; i >= 0; i--) {
            for (j = 0; j < i; j++) {
                muGS[k][j] -= round(muGS[k][i]) * muGS[i][j];
            }
        }
    
        //some other operations that don't change the value of muGS
        k++;
    }
    

    输出:

    1 muGS[k] diff:0.157396
    1 muGS[k] before:0.288172
    1 muGS[k] after:0.288172
    1 muGS[k] mult before:0.0171023
    1 muGS[k] mult after:0.083043
    1 muGS[k] abs before:0.130776
    1 muGS[k] abs after:0.288172
    

    发生的另一件事是 before2 之前2 . 有没有可能是我的精度降低了,或者为什么会这样?

    谢谢

    1 回复  |  直到 10 年前
        1
  •  4
  •   Iurii    10 年前

    没有精度损失。您的代码中有一个错误:

            cout << "1 muGS[k] before:" << muGS[k][0] << endl;
            cout << "1 muGS[k] after:" << muGS[k][0] << endl;
    

    前后打印相同的值。 但应该是:

            cout << "1 muGS[k] before:" << before2 << endl;
            cout << "1 muGS[k] after:" << muGS[k][0] << endl;