Reduced Row Echelon Implementation in Java

Question

Via this code:

import java.util.Arrays;

public class RREFS {

    public static void main(String[] args){

        double[][] rref;
        double[][] matrix;
        matrix=new double[][]{
                new double[]{  8.385956100787163E-8,   0.9103664774626016, 0.41380294430118253,                 0.0, 0.0, 0.0  },
                new double[]{ 2.0779736498353668E-8,  0.22558161873553792, -0.4962795612181877,  0.8383433249008051, 0.0, 0.0  },
                new double[]{ 1.0081125874457642E-7, -0.34690893617920376,  0.7631996595942279,  0.5451426139958772, 0.0, 0.0  },
                new double[]{                   0.0,                  0.0,                 0.0,                 0.0, 1.0, 0.0  },
                new double[]{                   0.0,                  0.0,                 0.0,                 0.0, 0.0, 1.0  }
        };
        rref=RREFS.rref(matrix);
        print(rref);

//      this
//      [ 8.385956100787163E-8,   0.9103664774626016, 0.41380294430118253,                 0.0, 0.0, 0.0]
//      [2.0779736498353668E-8,  0.22558161873553792, -0.4962795612181877,  0.8383433249008051, 0.0, 0.0]
//      [1.0081125874457642E-7, -0.34690893617920376,  0.7631996595942279,  0.5451426139958772, 0.0, 0.0]
//      [                  0.0,                  0.0,                 0.0,                 0.0, 1.0, 0.0]
//      [                  0.0,                  0.0,                 0.0,                 0.0, 0.0, 1.0]

//       produces
//       [1.0, 0.0, 0.0,                 0.0, 0.0, 0.0]
//       [0.0, 1.0, 0.0,                 0.0, 0.0, 0.0]
//       [0.0, 0.0, 1.0, -1.3999999999999997, 0.0, 0.0]
//       [0.0, 0.0, 0.0,                 0.0, 1.0, 0.0]
//       [0.0, 0.0, 0.0,                 0.0, 0.0, 1.0]

//       insteada
//       [1.0, 0.0, 0.0,  1.38165312352941182995E7, 0.0, 0.0]
//       [0.0, 1.0, 0.0,        -0.636363636363636, 0.0, 0.0]
//       [0.0, 0.0, 1.0,                      -1.4, 0.0, 0.0]
//       [0.0, 0.0, 0.0,                         0, 1.0, 0.0]
//       [0.0, 0.0, 0.0,                         0, 0.0, 1.0]

    }


    /**
     * Puts a matrix into reduced row echelon form
     *
     * @param matrix input matrix
     *
     * @return 2D result matrix
     */
    public static double[][] rref(double[][] matrix){
        int columnIndex = 0;
        int cursor;

        // number of rows and columns in matrix
        int getRowSize = matrix.length;
        int getColumnSize = matrix[0].length;


        loop:
        for(int rowIndex = 0; rowIndex < getRowSize; rowIndex++){
            if(getColumnSize <= columnIndex){
                break loop;
            }
            cursor = rowIndex;
            while(matrix[cursor][columnIndex] == 0){
                cursor++;
                if(getRowSize == cursor){
                    cursor = rowIndex;
                    columnIndex++;
                    if(getColumnSize == columnIndex){
                        break loop;
                    }
                }

            }

            matrix = rowSwap(matrix, cursor, rowIndex);
            if(matrix[rowIndex][columnIndex] != 0){
                matrix = rowScale(matrix, rowIndex, (1/matrix[rowIndex][columnIndex]));
            }

            for(cursor = 0; cursor < getRowSize; cursor++){
                if(cursor != rowIndex){
                    matrix = rowAddScale(matrix, rowIndex, cursor,((-1) * matrix[cursor][columnIndex]));
                }
            }columnIndex++;
        }return matrix;
    }

    /**
     * Swap positions of 2 rows
     *
     * @param matrix matrix before row additon
     * @param rowIndex1 int index of row to swap
     * @param rowIndex2 int index of row to swap
     *
     * @return matrix after row swap
     */
    private static double[][] rowSwap(double[][] matrix, int rowIndex1,
                                      int rowIndex2){
        // number of columns in matrix
        int numColumns = matrix[0].length;

        // holds number to be swapped
        double hold;

        for(int k = 0; k < numColumns; k++){
            hold = matrix[rowIndex2][k];
            matrix[rowIndex2][k] = matrix[rowIndex1][k];
            matrix[rowIndex1][k] = hold;
        }

        return matrix;
    }

    /**
     * Adds 2 rows together row2 = row2 + row1
     *
     * @param matrix matrix before row additon
     * @param rowIndex1 int index of row to be added
     * @param rowIndex2 int index or row that row1 is added to
     *
     * @return matrix after row addition
     */
    private static double[][] rowAdd(double[][] matrix, int rowIndex1,
                                     int rowIndex2){
        // number of columns in matrix
        int numColumns = matrix[0].length;

        for(int k = 0; k < numColumns; k++){
            matrix[rowIndex2][k] += matrix[rowIndex1][k];
        }

        return matrix;
    }

    /**
     * Multiplies a row by a scalar
     *
     * @param matrix matrix before row additon
     * @param rowIndex int index of row to be scaled
     * @param scalar double to scale row by
     *
     * @return matrix after row scaling
     */
    private static double[][] rowScale(double[][] matrix, int rowIndex,
                                       double scalar){
        // number of columns in matrix
        int numColumns = matrix[0].length;

        for(int k = 0; k < numColumns; k++){
            matrix[rowIndex][k] *= scalar;
        }

        return matrix;
    }

    /**
     * Adds a row by the scalar of another row
     * row2 = row2 + (row1 * scalar)
     * @param matrix matrix before row additon
     * @param rowIndex1 int index of row to be added
     * @param rowIndex2 int index or row that row1 is added to
     * @param scalar double to scale row by
     *
     * @return matrix after row addition
     */
    private static double[][] rowAddScale(double[][] matrix, int rowIndex1,
                                          int rowIndex2, double scalar){
        // number of columns in matrix
        int numColumns = matrix[0].length;

        for(int k = 0; k < numColumns; k++){
            matrix[rowIndex2][k] += (matrix[rowIndex1][k] * scalar);
        }

        return matrix;
    }

    public static void print(double[][] matrix) {
        print("",matrix);
    }
    public static void print(String message, double[][] matrix) {
        System.out.println(message);
        for (int i = 0; i < matrix.length; i++) {
            System.out.println(Arrays.toString(matrix[i]));
        }System.out.println();
    }
}

I am trying to build a custom rref calculator thus I have edie-zhou's implementation asper.

It works with every simple matrix I have tested so far, However the result of

 [ 8.385956100787163E-8,   0.9103664774626016, 0.41380294430118253,                 0.0, 0.0, 0.0]
 [2.0779736498353668E-8,  0.22558161873553792, -0.4962795612181877,  0.8383433249008051, 0.0, 0.0]
 [1.0081125874457642E-7, -0.34690893617920376,  0.7631996595942279,  0.5451426139958772, 0.0, 0.0]
 [                  0.0,                  0.0,                 0.0,                 0.0, 1.0, 0.0]
 [                  0.0,                  0.0,                 0.0,                 0.0, 0.0, 1.0]

is this

[1.0, 0.0, 0.0,                 0.0, 0.0, 0.0]
[0.0, 1.0, 0.0,                 0.0, 0.0, 0.0]
[0.0, 0.0, 1.0, -1.3999999999999997, 0.0, 0.0]
[0.0, 0.0, 0.0,                 0.0, 1.0, 0.0]
[0.0, 0.0, 0.0,                 0.0, 0.0, 1.0]

instead of

[1.0, 0.0, 0.0,  1.38165312352941182995E7, 0.0, 0.0]
[0.0, 1.0, 0.0,        -0.636363636363636, 0.0, 0.0]
[0.0, 0.0, 1.0,                      -1.4, 0.0, 0.0]
[0.0, 0.0, 0.0,                         0, 1.0, 0.0]
[0.0, 0.0, 0.0,                         0, 0.0, 1.0]

according to this online calculator and some others I tested, please what seems to be the problem with the code.

Please don't link to external code. Include all relevant code in the question itself. If you decide to take down the gist in future, this question will become useless for future readers. — Michael, Jan 31 '23 at 10:12
Happy to remove it when the question contains all the relevant code. If Stack Overflow complains that the ratio of text to code is too low, then you'll either need to add more explanation, or figure out how you can distill your problem to use less code. Once again, a gist is a poor solution for the reason I already described. — Michael, Feb 01 '23 at 12:08
@Michael Tbh there’s no other explanation I can possibly give other than how I have simply explained the problem. Since you have decided to stay unhappy until your wish is granted, I can only promise that my gist would not be taken down anytime in the future just the same way concerning some of your posts on Stack Overflow with external links you have no problem in believing associated websites or services where they are hosted won’t be taken down in the future. I only want you to be happy, that’s why I typed this much. I just really hope you can be happy though… after all I have typed — linker, Feb 01 '23 at 16:08
It's not about me being unhappy, it's about no one person being above the rules. The rules exist to make things better for everyone. "*just the same way concerning some of your posts on Stack Overflow with external links*". No, this is entirely different. Your question is *entirely predicated* on the existence of that link. My answers may be supplemented, but even if the links go down, the answers should stand on their own. If some don't stand on their own, then feel free to downvote them and comment on how they can be improved. — Michael, Feb 01 '23 at 16:14
OK, the code is included in the question now. Perhaps someone is nice enough to go through it and spot the bug, but I think you'll have more luck by either commenting under the Gist to see if the original author knows what's going wrong, or by manually calculating the matrix operation and then stepping through the code to see if you spot the error yourself. Anyway, best of luck! — Bart Kiers, Feb 02 '23 at 19:28
If you print the matrix after each pass, after pass two it has some alarming numbers in it (E16, E23). If this gives correct results for some matrices but not others, I would suspect loss of precision due to numerical instability of the chosen method. — teapot418, Feb 02 '23 at 20:39
@teapot418 I switched from working with floats initially to double because after each pass I noticed some funny numbers were showing up occasionally. Please what do you recommend in order to get better precision: BigDecimal? — linker, Feb 03 '23 at 13:59
I think LU decomposition is usually preferred as being less explody. I think you can get rref from LU, but I won't point you to a specific formula, this topic is complex and my memory of it fuzzy. Or you can cheat and [use fractions](https://rosettacode.org/wiki/Reduced_row_echelon_form#Java). — teapot418, Feb 03 '23 at 14:05
The above probably rather with BigFraction, Fraction is int-based and can overflow. — teapot418, Feb 03 '23 at 14:19
As I understand RREF is an ill-posed problem so if you start from the double-matrix as input you might work a lot to give better results, while ignoring that the input matrix may have been rounded to double precision and that rounding may cause differences as big as the one we see. — Hans Olsson, Feb 03 '23 at 14:30

score 3 · Accepted Answer · answered Feb 03 '23 at 14:35

3

You may actually have a "do not compare floats for equality" problem in there. If I do:

public static boolean isZero(double foo) {
    // demonstration only, replace with less silly implementation
    return Math.abs(foo)<0.0000001;
}

  ...
  while(isZero(matrix[cursor][columnIndex])){
  //instead of
  while(matrix[cursor][columnIndex] == 0){

I get

[1.0, 0.0, 0.0, 1.3816531235294117E7, 0.0, 0.0]
[0.0, 1.0, 0.0, -0.6363636363636362, 0.0, 0.0]
[0.0, 0.0, 1.0, -1.3999999999999997, 0.0, 0.0]
[0.0, 0.0, 0.0, 0.0, 1.0, 0.0]
[0.0, 0.0, 0.0, 0.0, 0.0, 1.0]

which is the result you wanted I believe.

More discussion e.g. in How should I do floating point comparison?

answered Feb 03 '23 at 14:35

teapot418

1,239
1
3
9

you just made me realize that java is more than just a tea. big thanks @teapot. You’re the real OG – linker Feb 03 '23 at 15:42
1

In addition to this, if you're doing math you may find the [BigDecimal](https://docs.oracle.com/en/java/javase/17/docs/api/java.base/java/math/BigDecimal.html) class to be preferable over floating point implementations, which can suffer from inaccuracies due to rounding that may compound on each other. If you use a BigDecimal, you can check if an instance of it is 0 with `foo.signum() == 0` – kagof Feb 08 '23 at 18:16

Reduced Row Echelon Implementation in Java

1 Answers1