c# – 在负面情况下几乎正确的矩阵分解代码失败

代码的内容取自here

除了我的Decompose方法之外,一切看起来都很完美(有关LU分解矩阵意义的更多信息,请查看here)：

public Matrix Decompose(out int[] permutationArray,out int toggle)
    {
        // Doolittle LUP decomposition with partial pivoting.
        // rerturns: result is L (with 1s on diagonal) and U; perm holds row permutations; toggle is +1 or -1 (even or odd)

        if (!this.IsSquare())
        {
            throw new Exception("LU-Decomposition can only be found for a square matrix");
        }

        Matrix decomposedMatrix = new Matrix(this); // copy this matrix before messing with it

        permutationArray = new int[NumberOfRows]; // set up row permutation result
        for (int i = 0; i < NumberOfRows; ++i)
        {
            permutationArray[i] = i;
        }

        toggle = 1; // toggle tracks row swaps. +1 -> even,-1 -> odd. used by MatrixDeterminant

        for (int columnIndex = 0; columnIndex < NumberOfRows - 1; columnIndex++) // each column
        {
            // find largest value in col j
            double maxOfColumn = Math.Abs(decomposedMatrix.GetElement(columnIndex,columnIndex));

            int pivotRowIndex = columnIndex;

            for (int rowIndex = columnIndex + 1; rowIndex < NumberOfRows; rowIndex++)
            {
                if (Math.Abs(decomposedMatrix.GetElement(rowIndex,columnIndex)) > maxOfColumn)
                {
                    maxOfColumn = decomposedMatrix.GetElement(rowIndex,columnIndex);
                    pivotRowIndex = rowIndex;
                }
            }

            if (pivotRowIndex != columnIndex) // if largest value not on pivot,swap rows
            {
                double[] rowPtr = decomposedMatrix.GetRow(pivotRowIndex);
                decomposedMatrix.SetRowOfMatrix(pivotRowIndex,decomposedMatrix.GetRow(columnIndex));
                decomposedMatrix.SetRowOfMatrix(columnIndex,rowPtr);

                int tmp = permutationArray[pivotRowIndex]; // and swap permutation info
                permutationArray[pivotRowIndex] = permutationArray[columnIndex];
                permutationArray[columnIndex] = tmp;

                toggle = -toggle; // adjust the row-swap toggle
            }

            if (decomposedMatrix.GetElement(columnIndex,columnIndex) == 0.0)
            {
                // find a good row to swap
                int goodRow = -1;
                for (int row = columnIndex + 1; row < decomposedMatrix.NumberOfRows; row++)
                {
                    if (decomposedMatrix.GetElement(row,columnIndex) != 0.0)
                    {
                        goodRow = row;
                    }
                }

                if (goodRow == -1)
                {
                    throw new Exception("Cannot use Doolittle's method on this matrix");
                }

                // swap rows so 0.0 no longer on diagonal
                double[] rowPtr = decomposedMatrix.GetRow(goodRow);
                decomposedMatrix.SetRowOfMatrix(goodRow,rowPtr);

                int tmp = permutationArray[goodRow]; // and swap perm info
                permutationArray[goodRow] = permutationArray[columnIndex];
                permutationArray[columnIndex] = tmp;

                toggle = -toggle; // adjust the row-swap toggle

            }

            for (int rowIndex = columnIndex + 1; rowIndex < this.NumberOfRows; ++rowIndex)
            {
                double valueToStore = decomposedMatrix.GetElement(rowIndex,columnIndex) / decomposedMatrix.GetElement(columnIndex,columnIndex);
                decomposedMatrix.SetElement(rowIndex,columnIndex,valueToStore);
                for (int nextColumnIndex = columnIndex + 1; nextColumnIndex < NumberOfRows; ++nextColumnIndex)
                {
                    double valueToStore2 = decomposedMatrix.GetElement(rowIndex,nextColumnIndex) -
                                            decomposedMatrix.GetElement(rowIndex,columnIndex) * decomposedMatrix.GetElement(columnIndex,nextColumnIndex);
                    decomposedMatrix.SetElement(rowIndex,nextColumnIndex,valueToStore2);
                }
            }

        } // main j column loop

        return decomposedMatrix;
    } // MatrixDecompose

这是它失败的单元测试:(请注意！如果开始矩阵使用正数而不是负数,则完全相同的单元测试通过)

[TestMethod()]
    public void Matrix_DecomposeTest3_DifferentSigns()
    {
        Matrix matrixA = new Matrix(9);

        //Set up matrix A
        double[] row1OfMatrixA = { 3,-7,0 };
        double[] row2OfMatrixA = { 0,1,8,-4,2,0 };
        double[] row3OfMatrixA = { 0,-9,3,0 };
        double[] row4OfMatrixA = { -5,4,7,-1,3 };
        double[] row5OfMatrixA = { 0,-3,-8,0 };
        double[] row6OfMatrixA = { 0,5,-2,4 };
        double[] row7OfMatrixA = { 0,7 };
        double[] row8OfMatrixA = { 0,0 };
        double[] row9OfMatrixA = { 0,6 };

        matrixA.SetRowOfMatrix(0,row1OfMatrixA);
        matrixA.SetRowOfMatrix(1,row2OfMatrixA);
        matrixA.SetRowOfMatrix(2,row3OfMatrixA);
        matrixA.SetRowOfMatrix(3,row4OfMatrixA);
        matrixA.SetRowOfMatrix(4,row5OfMatrixA);
        matrixA.SetRowOfMatrix(5,row6OfMatrixA);
        matrixA.SetRowOfMatrix(6,row7OfMatrixA);
        matrixA.SetRowOfMatrix(7,row8OfMatrixA);
        matrixA.SetRowOfMatrix(8,row9OfMatrixA);

        Console.Out.Write(matrixA);

        //The LUP Decomposition

        //Correct L Part
        Matrix correctLPartOfLUPDecomposition = new Matrix(9);

        double[] row1OfLMatrix = { 1,0 };
        double[] row2OfLMatrix = { 0,0 };
        double[] row3OfLMatrix = { 0,.143,0 };
        double[] row4OfLMatrix = { -.6,0 };
        double[] row5OfLMatrix = { 0,0 };
        double[] row6OfLMatrix = { 0,.435,0 };
        double[] row7OfLMatrix = { 0,-.217,-.5,0 };
        double[] row8OfLMatrix = { 0,-.429,-.796,-.342,-.319,.282,-.226,0 };
        double[] row9OfLMatrix = { 0.000,0.286,0.963,0.161,0.523,0.065,0.013,0.767,1.000 };

        correctLPartOfLUPDecomposition.SetRowOfMatrix(0,row1OfLMatrix);
        correctLPartOfLUPDecomposition.SetRowOfMatrix(1,row2OfLMatrix);
        correctLPartOfLUPDecomposition.SetRowOfMatrix(2,row3OfLMatrix);
        correctLPartOfLUPDecomposition.SetRowOfMatrix(3,row4OfLMatrix);
        correctLPartOfLUPDecomposition.SetRowOfMatrix(4,row5OfLMatrix);
        correctLPartOfLUPDecomposition.SetRowOfMatrix(5,row6OfLMatrix);
        correctLPartOfLUPDecomposition.SetRowOfMatrix(6,row7OfLMatrix);
        correctLPartOfLUPDecomposition.SetRowOfMatrix(7,row8OfLMatrix);
        correctLPartOfLUPDecomposition.SetRowOfMatrix(8,row9OfLMatrix);

        //Correct U Part
        Matrix correctUPartOfLUPDecomposition = new Matrix(9);

        double[] row1OfUMatrix = { -5.000,0.000,4.000,7.000,-1.000,3.000 };
        double[] row2OfUMatrix = { 0.000,2.000,5.000,-3.000,-2.000,4.000 };
        double[] row3OfUMatrix = { 0.000,7.714,-0.714,-3.571,1.000,-0.571 };
        double[] row4OfUMatrix = { 0.000,-4.600,4.200,-0.600,1.800 };
        double[] row5OfUMatrix = { 0.000,-9.000,3.000,0.000 };
        double[] row6OfUMatrix = { 0.000,-9.826,4.261,6.217 };
        double[] row7OfUMatrix = { 0.000,9.000,9.500 };
        double[] row8OfUMatrix = { 0.000,-2.043,2.272 };
        double[] row9OfUMatrix = { 0.000,-3.153 };

        correctUPartOfLUPDecomposition.SetRowOfMatrix(0,row1OfUMatrix);
        correctUPartOfLUPDecomposition.SetRowOfMatrix(1,row2OfUMatrix);
        correctUPartOfLUPDecomposition.SetRowOfMatrix(2,row3OfUMatrix);
        correctUPartOfLUPDecomposition.SetRowOfMatrix(3,row4OfUMatrix);
        correctUPartOfLUPDecomposition.SetRowOfMatrix(4,row5OfUMatrix);
        correctUPartOfLUPDecomposition.SetRowOfMatrix(5,row6OfUMatrix);
        correctUPartOfLUPDecomposition.SetRowOfMatrix(6,row7OfUMatrix);
        correctUPartOfLUPDecomposition.SetRowOfMatrix(7,row8OfUMatrix);
        correctUPartOfLUPDecomposition.SetRowOfMatrix(8,row9OfUMatrix);

        //Calculate values for the above
        int[] permutationArray;
        int toggleValue;
        Matrix decomposeResult = matrixA.Decompose(out permutationArray,out toggleValue);

        Matrix calculatedLPartOfLUPDecomposition = decomposeResult.ExtractLower();
        Matrix calculatedUPartOfLUPDecomposition = decomposeResult.ExtractUpper();

        //Compare the two matrices
        double tolerance = .001;
        Matrix LDifferenceMatrix = calculatedLPartOfLUPDecomposition - correctLPartOfLUPDecomposition;
        Matrix UDifferenceMatrix = calculatedUPartOfLUPDecomposition - correctUPartOfLUPDecomposition;

        for (int i = 0; i < LDifferenceMatrix.NumberOfRows; i++)
        {
            for (int j = 0; j < LDifferenceMatrix.NumberOfColumns; j++)
            {
                (Math.Abs(LDifferenceMatrix.GetElement(i,j)) < tolerance).Should().BeTrue();
            }
        }

        for (int i = 0; i < UDifferenceMatrix.NumberOfRows; i++)
        {
            for (int j = 0; j < UDifferenceMatrix.NumberOfColumns; j++)
            {
                (Math.Abs(UDifferenceMatrix.GetElement(i,j)) < tolerance).Should().BeTrue();
            }
        }
    }

我使用找到的工具here得到了单元测试的数字.

我有单元测试Matrix的所有其他方法,他们都通过.

我对分解方法有什么看法？