0

提前抱歉。这将是一个漫长的过程。

一些背景。我已经在 J​​ava 中实现了几个与矩阵相关的算法来获得乐趣。我正在研究简化的行梯形式（rref）。我已经实现了使用和不使用部分旋转的 rref。我的理解是，部分旋转应该比不旋转在数值上更稳定。但是，我看到一些矩阵与我的实现相反，这让我认为我做错了什么。

关于我的实现的一些注意事项：

  • 我已经实现了自己的Number课程。此类允许使用复数或实数。
  • 类中的compareTo()方法Number使用复数的大小。实际值的比较是通过另一种方法完成的。
  • 我的Matrix类本质上包装了这些Number对象的二维数组，并提供了几个构造函数/方法。

代码：
没有部分旋转：

public static Matrix rrefNoPivot(Matrix A) {
    Number m, scale;
    
    int pivotRow = 0,
        pivotCol = 0;
    
    while(pivotRow<A.m && pivotCol<A.n) {
        if(!A.entries[pivotRow][pivotCol].equals(Number.ZERO)) { 
            scale = Number.divide(Number.ONE, A.entries[pivotRow][pivotCol]);
             
            for(int k=pivotCol; k<A.n; k++) { // scale the whole row
                A.entries[pivotRow][k] = Number.multiply(A.entries[pivotRow][k], scale); 
            }   
        }
            
        for(int i=0; i<A.m; i++) {
            m = A.entries[i][pivotCol];
            
            if(pivotRow != i) {
                for(int k=pivotCol; k<A.n; k++) {
                    A.entries[i][k] = Number.subtract(A.entries[i][k], 
                            Number.multiply(A.entries[pivotRow][k], m));
                }
            }
        }
        
        pivotRow++;
        pivotCol++;
    }
    
    return A;
}

部分旋转：

public static Matrix rref(Matrix A) {
    Number mult, scale, currentMax;
    int maxIndex;
    
    int pivotRow = 0,
        pivotCol = 0;
    
    while(pivotRow<A.m && pivotCol<A.n) {   
        maxIndex = pivotRow;
        currentMax = A.entries[pivotRow][pivotCol];
        
        for(int i=pivotRow; i<A.m; i++) { // find the maximum entry in the pivot column (at or below the pivot ).
            if(A.entries[i][pivotCol].compareTo(currentMax) > 0) {
                maxIndex = i;
                currentMax = A.entries[i][pivotCol];
            }
        }
        
        if(!A.entries[maxIndex][pivotCol].equals(Number.ZERO)) { // Check that the maximum absolute value is not zero.
            if(pivotRow != maxIndex) {
                A = A.swapRows(pivotRow, maxIndex); // Make the row with the largest value in the pivot column the pivot for this row.
            }
            
            scale = Number.divide(Number.ONE, A.entries[pivotRow][pivotCol]);
             
            for(int k=pivotCol; k<A.n; k++) { // scale the whole row
                A.entries[pivotRow][k] = Number.multiply(A.entries[pivotRow][k], scale); 
            }
            
            for(int i=0; i<A.m; i++) {
                mult = A.entries[i][pivotCol];
                
                if(pivotRow != i) {
                    for(int k=pivotCol; k<A.n; k++) {
                        A.entries[i][k] = Number.subtract(A.entries[i][k], 
                                Number.multiply(A.entries[pivotRow][k], mult));
                    }
                }
            }
            
            pivotRow++;
            pivotCol++;
        }
        else { // Then we do not have a pivot for this column (i.e. the column is all zeros).
            pivotCol++;
        }
    }
    
    return A;
}

现在来看看结果。以下是不同矩阵的三个结果。

    E: rref(E): rrefNoPivot(E):
    [ [1 2 3] [ [1 0 0] [ [1 0 -1]
      [4 5 6] [0 1 0] [0 1 2]
      [7 8 9] [0 0 1] [0 0 0]
      [10 11 12] [0 0 0] [0 0 0]
      [13 14 15] ] [0 0 0] ] [0 0 0 ] ]
    
    
    C: rref(C): rrefNoPivot(C):
    [ [1 2 3 4 5 ] [ [1 0 0 -0.9333333333333338 -1] [ [1 0 -1 -2 -3]
      [6 7 8 9 10] [0 1 0 0.8666666666666671 0 ] [0 1 2 3 4]
      [11 12 13 14 15] ] [0 0 1 1.0666666666666667 2 ] ] [0 0 0 0 0 ] ]
    
    
    B: rref(B): rrefNoPivot(B):
    [ [1 2 3 1 2 3 4 5 6] [ [1 0 0 0 0 0 0 0 0] [ [1 0 0 0 0 0 0 2.220446049250313E-16 0]
      [4 5 6 3 4 5 7 8 2] [0 1 0 0 0 0 0 0 0] [0 1 0 0 0 0 0 -4.440892098500626E-16 0]
      [1 5 5 2 6 7 9 0 1] [0 0 1 0 0 0 0 0 0] [0 0 1 0 0 0 0 4.440892098500626E-16 0]
      [3 4 5 2 6 7 8 9 2] [0 0 0 1 0 0 0 0 0] [0 0 0 1 0 0 0 -1.3877787807814457E-17 0]
      [1 1 1 3 4 7 8 9 1] [0 0 0 0 1 0 0 0 0] [0 0 0 0 1 0 0 -4.440892098500626E-16 0]
      [3 4 7 8 3 1 2 3 4] [0 0 0 0 0 1 0 0 0] [0 0 0 0 0 1 0 8.881784197001252E-16 0]
      [3 5 6 8 1 3 5 9 1] [0 0 0 0 0 0 1 0 0] [0 0 0 0 0 0 1 -4.440892098500626E-16 0]
      [0 4 5 3 2 0 1 2 3] [0 0 0 0 0 0 0 1 0] [0 0 0 0 0 0 0 0.9999999999999999 0]
      [1 4 5 1 7 0 5 3 1] ] [0 0 0 0 0 0 0 0 1] ] [0 0 0 0 0 0 0 7.105427357601002E-15 1] ]

对于前两个，rref()甚至没有得到正确的答案。我查看了矩阵，在 while 循环的倒数第二次迭代中，这就是它的样子……

    [[1 0 -1.0000000000000002]
      [0 1 2]
      [0 0 1.7763568394002505E-15]
      [0 0 0 ]
      [0 0 8.881784197001252E-16]]

所以有非常小的非零值（应该是零）导致了这个问题。这似乎是浮点算术错误的结果。C 类似。然而，对于矩阵 B，部分旋转确实产生了正确的答案，而不旋转方法有一些错误。

所以这是我的正式问题，最后是嘘。首先，我的部分枢轴执行是否存在问题？无论是引入数值不稳定性还是不正确的逻辑。其次，部分旋转是否保证在数值上更稳定？或者，使用部分旋转时的两个矩阵 E 和 C 示例是否不太稳定？这些类型的错误是不可避免的吗？我应该将非常小的数字四舍五入为零以便算法有效吗？
提前谢谢大家。

4

0 回答 0