All geek questions in one place

Fixed incorrect matrix decomposition code in the negative case

I am writing my own matrix class with the methods you expect to accompany a matrix class.

Code gits were taken from here

Everything seems to work just fine, except for my Decompose method (for more information on what LU means, expand the matrix, see 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 method on this matrix"); } // swap rows so 0.0 no longer on diagonal double[] rowPtr = decomposedMatrix.GetRow(goodRow); decomposedMatrix.SetRowOfMatrix(goodRow, decomposedMatrix.GetRow(columnIndex)); decomposedMatrix.SetRowOfMatrix(columnIndex, 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

This is the Unit Test, in which it fails: (Note! The exact same Unit Test passes if the initial matrix uses positive numbers instead of negative numbers)

  [TestMethod()] public void Matrix_DecomposeTest3_DifferentSigns() { Matrix matrixA = new Matrix(9); //Set up matrix A double[] row1OfMatrixA = { 3, 0, 0, -7, 0, 0, 0, 0, 0 }; double[] row2OfMatrixA = { 0, 1, 8, 0, -4, 2, 0, 0, 0 }; double[] row3OfMatrixA = { 0, 2, 8, 0, -9, 3, 0, 0, 0 }; double[] row4OfMatrixA = { -5, 0, 0, 4, 0, 7, -1, 0, 3 }; double[] row5OfMatrixA = { 0, -3, -7, 0, 7, -8, 0, -3, 0 }; double[] row6OfMatrixA = { 0, 7, 2, 5, -3, 7, -2, 0, 4 }; double[] row7OfMatrixA = { 0, 0, 0, -2, 0, -8, 4, 0, 7 }; double[] row8OfMatrixA = { 0, 0, 0, 0, -9, 0, 0, 3, 0 }; double[] row9OfMatrixA = { 0, 0, 0, 1, 0, 4, 7, 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, 0, 0, 0, 0, 0, 0, 0 }; double[] row2OfLMatrix = { 0, 1, 0, 0, 0, 0, 0, 0, 0 }; double[] row3OfLMatrix = { 0, .143, 1, 0, 0, 0, 0, 0, 0 }; double[] row4OfLMatrix = { -.6, 0, 0, 1, 0, 0, 0, 0, 0 }; double[] row5OfLMatrix = { 0, 0, 0, 0, 1, 0, 0, 0, 0 }; double[] row6OfLMatrix = { 0, 0, 0, .435, 0, 1, 0, 0, 0 }; double[] row7OfLMatrix = { 0, 0, 0, -.217, 0, -.5, 1, 0, 0 }; double[] row8OfLMatrix = { 0, -.429, -.796, -.342, -.319, .282, -.226, 1, 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, 0.000, 4.000, 0.000, 7.000, -1.000, 0.000, 3.000 }; double[] row2OfUMatrix = { 0.000, 7.000, 2.000, 5.000, -3.000, 7.000, -2.000, 0.000, 4.000 }; double[] row3OfUMatrix = { 0.000, 0.000, 7.714, -0.714, -3.571, 1.000, 0.286, 0.000, -0.571 }; double[] row4OfUMatrix = { 0.000, 0.000, 0.000, -4.600, 0.000, 4.200, -0.600, 0.000, 1.800 }; double[] row5OfUMatrix = { 0.000, 0.000, 0.000, 0.000, -9.000, 0.000, 0.000, 3.000, 0.000 }; double[] row6OfUMatrix = { 0.000, 0.000, 0.000, 0.000, 0.000, -9.826, 4.261, 0.000, 6.217 }; double[] row7OfUMatrix = { 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 9.000, 0.000, 9.500 }; double[] row8OfUMatrix = { 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, -2.043, 2.272 }; double[] row9OfUMatrix = { 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 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(); } } }

I got the Unit Test numbers using the tool I found here .

I have Unit Checked all other methods for Matrix, and they all pass.

What am I missing in my decomposition method?

+4
source share
2 answers

If you select the largest value of the resulting column:

 maxOfColumn = decomposedMatrix.GetElement(rowIndex, columnIndex);

You really have to choose the largest absolute value found there (you do this in the if statement and seem to have forgotten after that). I would do it explicitly by adding a comment so that the next person, to read this code, can see it more clearly.

  //if a value with a larger absolute value has been found replace the currently selectedmax if (Math.Abs(decomposedMatrix.GetElement(rowIndex,columnIndex)) > maxOfColumn) { maxOfColumn = Math.Abs(decomposedMatrix.GetElement(rowIndex, columnIndex));
+1
source

it

 if (Math.Abs(decomposedMatrix.GetElement(rowIndex,columnIndex)) > maxOfColumn) { maxOfColumn = decomposedMatrix.GetElement(rowIndex, columnIndex);

it should be:

 if (Math.Abs(decomposedMatrix.GetElement(rowIndex,columnIndex)) > maxOfColumn) { maxOfColumn = Math.Abs(decomposedMatrix.GetElement(rowIndex, columnIndex));
0
source

Source: https://habr.com/ru/post/1494329/

More articles:


All Articles