Linear algebra, reduced row echelon form - function 1

Question

This is a derivative post from here.

This is just a general review, so the question is the same:

Is there something...

1. That you would consider as a bad practice and why?
2. That is just bad in some other way?

This question applies to the function toRowEchelonForm(matrix, returnMethod).

(Just to be clear, this is not a production code. There are too few tests and I do not guarantee it does what it is supposed to do. WolframAlpha limits its inputs and I haven't found alternative yet.)

/*jslint browser: true, indent: 8 */
/*global console */

function toRowEchelonForm(matrix, returnMethod) {
        'use strict';
        var f, pivotRow, toReduceRow, scaledRow, len, pos, new_matrix;

        function sortRows(matrix) {
                var pivots, c;
                // map rows to arrays with the structure
                // [<column index>, <pivot coefficient>, <row index>]
                // which makes them easily sortable
                pivots = matrix.map(function (row, r) {
                        var len = row.length;

                        for (c = 0; c < len; c += 1) {
                                if (row[c]) {
                                        return [c, row[c], r]; // found a non-zero coefficient at matrix[r][c]
                                }
                        }

                        return [Infinity, null, r]; // row is all-zero
                });

                // sort the pivots array, and use it to (re)build
                // the matrix with the rows in the correct order
                return pivots.sort().map(function (pivot) {
                        return matrix[pivot[2]];
                });
        }

        function leadingPivotsToTop(matrix) {
                var r, c, len, has_pivot, irrelevant, positions, new_matrix, count;

                len = {
                        row: matrix.length,
                        col: matrix[0].length
                };
                positions  = [];
                has_pivot  = [];
                irrelevant = [];
                new_matrix = [];
                count      = 0;

                // Find pivot positions
                for (c = 0; c < len.col; c += 1) {
                        for (r = 0; r < len.row; r += 1) {
                                if (matrix[r][c] === 1 && has_pivot[r] !== r) {
                                        has_pivot[r] = r;
                                        positions[positions.length] = r;
                                        break;
                                }
                        }
                }

                // Find irrelevant vectors positions
                for (r = 0; r < len.row; r += 1) {
                        if (has_pivot[r] === undefined) {
                                irrelevant[irrelevant.length] = r;
                                count += 1;
                        }
                }

                count = 0;

                // Sort positions
                for (r = 0; r < len.row; r += 1) {
                        if (matrix[positions[r]] !== undefined) {
                                new_matrix[r] = matrix[positions[r]];
                        } else {
                                new_matrix[r] = matrix[irrelevant[count]];
                                count += 1;
                        }
                }

                return new_matrix;
        }

        new_matrix = matrix;
        pos = {
                pivot: 0,
                reserved: []
        };
        len = {
                row: matrix.length,
                col: matrix[0].length
        };
        f = {
                getPivotPosition: function (pivotRow) {
                        var c;

                        for (c = 0; c < len.col; c += 1) {
                                if (new_matrix[pivotRow][c] !== 0
                                && pos.reserved[c] === undefined) {
                                        pos.reserved[c] = 1;
                                        return c;
                                }
                        }

                        return null;
                },
                reducePivotRow: function (pivotRow) {
                        var c, pivot;

                        if (new_matrix[pivotRow][pos.pivot] !== 1) {
                                pivot = new_matrix[pivotRow][pos.pivot];
                                for (c = 0; c < len.col; c += 1) {
                                        new_matrix[pivotRow][c] /= pivot;
                                }
                        }
                },
                scaleRow: function (pivotRow, toReduceRow) {
                        var c, row;

                        row = [];

                        for (c = 0; c < len.col; c += 1) {
                                row[c]  = new_matrix[pivotRow][c];
                                row[c] *= new_matrix[toReduceRow][pos.pivot];
                        }

                        return row;
                },
                rowReduction: function (toReduceRow, scaledRow) {
                        var c;

                        for (c = 0; c < len.col; c += 1) {
                                new_matrix[toReduceRow][c] -= scaledRow[c];
                        }
                }
        };

        for (pivotRow = 0; pivotRow < len.row; pivotRow += 1) {

                pos.pivot = null;

                pos.pivot = f.getPivotPosition(pivotRow);

                if (pos.pivot !== null) {
                        f.reducePivotRow(pivotRow);

                        for (
                                toReduceRow = 0;
                                toReduceRow < len.row;
                                toReduceRow += 1
                        ) {
                                if (toReduceRow !== pivotRow
                                && new_matrix[toReduceRow][pos.pivot] !== 0) {
                                        scaledRow = f.scaleRow(
                                                pivotRow,
                                                toReduceRow
                                        );
                                        f.rowReduction(
                                                toReduceRow,
                                                scaledRow
                                        );
                                }
                        }
                }
        }

        if (returnMethod === "raw") {
                return new_matrix;
        } else if (returnMethod === "lead to top") {
                return leadingPivotsToTop(new_matrix);
        } else if (returnMethod === "lead to top and zeroes to bottom") {
                return sortRows(new_matrix);
        } else {
                return leadingPivotsToTop(new_matrix);
        }
}

function compereMatrices(matrixA, matrixB) {
        'use strict';
        var i, j, len;

        if (matrixA.length !== matrixB.length) {
                return false;
        }

        len = {};

        len.i = matrixA.length;

        for (i = 0; i < len.i; i += 1) {
                if (matrixA[i].length !== matrixB[i].length) {
                        return false;
                }

                len.j = matrixA[i].length;

                for (j = 0; j < len.j; j += 1) {
                        if (matrixA[i][j] !== matrixB[i][j]) {
                                return false;
                        }
                }
        }

        return true;
}


// Tests.
var m = [
        [5, -7, -8, -4],
        [2, 8, -22, -55],
        [-3, 0, -36, 12]
];
// answer: http://www.wolframalpha.com/input/?i=solve+row+echelon+form+{{5%2C+-7%2C+-8%2C+-4}%2C{2%2C+8%2C+-22%2C+-55}%2C+{-3%2C+0%2C+-36%2C+12}}
var mr = [
        [1, 0, 0, -6.785219399538105],
        [0, 1, 0, -4.54041570438799],
        [0, 0, 1, 0.23210161662817538]
];

console.log(" ");
console.log("toRowEcholonForm test:         " + compereMatrices(toRowEchelonForm(m), mr));

// answer: http://www.wolframalpha.com/input/?i=solve+row+echelon+form+{{5%2C+-23%2C+2%2C+4%2C+5%2C+11}%2C{4%2C+-3%2C+6%2C+4%2C+5%2C+2}%2C{3%2C+7%2C+-18%2C+7%2C+9%2C+-6}%2C{4%2C+87%2C+-12%2C+7%2C+12%2C+6}%2C{5%2C+4%2C+7%2C+11%2C+7%2C+-7}}
var m = [
        [5, -23, 2, 4, 5, 11],
        [4, -3, 6, 4, 5, 2],
        [3, 7, -18, 7, 9, -6],
        [4, 87, -12, 7, 12, 6],
        [5, 4, 7, 11, 7, -7]
];

var mr = [
        [1, 0, 0, 0, 0, 10.784116921993304],
        [0, 1, 0, 0, 0, 0.3085998347488045],
        [0, 0, 1, 0, 0, -1.0969699432456959],
        [0, 0, 0, 1, 0, -1.369593780053366],
        [0, 0, 0, 0, 1, -5.630094680807834]
];

console.log(" ");
console.log("toRowEcholonForm test:         " + compereMatrices(toRowEchelonForm(m), mr));

// answer: http://www.wolframalpha.com/input/?i=solve+row+echelon+form+{{1%2C+2%2C+2%2C+2}%2C{1%2C+3%2C+3%2C+3}%2C+{1%2C+4%2C+16%2C+5}}
m = [
        [1, 2, 2, 2],
        [1, 3, 3, 3],
        [1, 4, 16, 5]
];

mr = [
        [1, 0, 0, 0],
        [0, 1, 0, 0.9166666666666666],
        [0, 0, 1, 0.08333333333333333]
];

console.log(" ");
console.log("toRowEcholonForm test:         " + compereMatrices(toRowEchelonForm(m), mr));

// answer: http://www.wolframalpha.com/input/?i=solve+row+echelon+form+{{0%2C+2%2C+-1%2C+-6}%2C{0%2C+3%2C+-2%2C+-16}%2C+{0%2C+0%2C+-3%2C+11}}
m = [
        [0, 2, -1, -6],
        [0, 3, -2, -16],
        [0, 0, -3, 11]
];

mr = [
        [0, 1, 0, 0],
        [0, 0, 1, 0],
        [0, 0, 0, 1]
];

console.log(" ");
console.log("toRowEcholonForm test:         " + compereMatrices(toRowEchelonForm(m), mr));

Answer 1

There's a lot of meat here; I'll try to do it justice. I've focussed on the core reduction logic, and haven't looked at the "post-processing" stuff that sorts the reduced matrix update: I lied - see bottom of answer.

Minor stuff

• There are a few constructions like this:

  pos.pivot = null;
  pos.pivot = f.getPivotPosition(pivotRow);
  

  but since f.getPivotPosition() defaults to returning null anyway, the first line is irrelevant. Not a huge concern by any stretch, but it made me go "huh?"

• There are a few strict comparisons against undefined. While possible, I'd prefer explicitly setting a boolean value to compare against, or using type coercion. For one, many browsers allow undefined to be set to anything, which can quickly break your code (the solution is to use typeof x === 'undefined'). But I'd also find it more straight-forward to read

  if(!pos.reserved[c])
  

  rather than

  if(pos.reserved[c] === undefined)
  

  because, semantically speaking, the condition being invoked is "if this position isn't reserved..." - not "if this position something we've never heard of...".
  Of course you want to be careful with type coercion, since zero is false'y - same as undefined, null or false itself - and you're working with indices that might well be zero. But in that case I'd explicitly set false or null (or -1, which is what indexOf() uses to mean "no index found"), and strictly check for that.

Bigger stuff

While you've extracted quite a few functions, they're very tied to an external state. You provide some of the values they need as arguments, but scaleRow (for instance), also relies on pos.pivot. So it's sort of "halfway extracted". It's not reusable, although the act of scaling a vector is actually a generic operation.

What you want is just to multiply every element in an array (aka row/vector) with a value, producing a new, scaled vector. Basically, \$\mathbf{V}\times c\$. You can achieve this like so:

 function scale(vector, factor) {
   // map returns a new array, produced by invoking the
   // iterator function on each element
   return vector.map(function (value) {
     return value * factor;
   });
 }

This function is completely independent, stateless and reusable. Of course, achieving these qualities isn't necessary for an inlined function (JS has closures for a reason: They're awesome!), but if you're going to extract a function, aim for maximum independence.

My take

I've left out the sorting functions, but this seems to work well. I've commented the code pretty heavily; hope it's clear

function reducedRowEchelonForm(matrix) {
  var knownPivotColumns = []; // this is our one piece of iffy state-keeping :(

  // Copy the input matrix (reusing the variable name) so we have a local copy to work on
  matrix = matrix.map(function (row) { return row.slice() });

  // Now, go through the matrix and do the reduction.
  // We're using forEach here, because we want to update the matrix
  // in-place, whereas `map()` will work on a separate instance
  matrix.forEach(function (row, rowIndex) {

    // Find the row's pivot
    // This is wrapped in an IIFE just for structure
    var pivot = (function () {
      // using a regular for-loop, since we may want to break out of it early
      for(var i = 0, l = row.length ; i < l ; i++ ) {
        if(!row[i] || knownPivotColumns[i]) { continue } // skip column if it's zero or its index is taken
        knownPivotColumns[i] = true;                     // "reserve" the column
        return { index: i, value: row[i] };              // return the pivot data
      }
      return null; // no pivot found
    }());

    // if there's no pivot, there's nothing else to do for this row
    if(!pivot) { return }

    // scale the row, if necessary
    if(pivot.value !== 1) {
      // using forEach as a "map in place" here
      row.forEach(function (_, i) { row[i] /= pivot.value });
    }

    // now reduce the other rows (calling them "siblings" here)
    matrix.forEach(function (sibling) {
      var siblingPivotValue = sibling[pivot.index];

      if(sibling === row || siblingPivotValue === 0) { return } // skip if possible

      // subtract the sibling-pivot-scaled row from the sibling
      // (another "forEach as map-in-place")
      sibling.forEach(function (_, i) { sibling[i] -= row[i] * siblingPivotValue });
    });
  });

  return matrix;
}

(jslint will complain about missing semi-colons,but I've only omitted them where it's "safe" - i.e. immediately before a })

You'll note that I haven't actually extracted any functions here. Everything is inlined, and there's even an IIFE. This was mostly to keep things self-contained for the purposes of this answer. Ideally, one would have a little library of generic vector/matrix functions (like scale() above) which could handle a lot of this stuff. A generic transform() function that alters an array in-place might also be handy (an extraction of the forEach-use above).

As for the sorting/post-processing stuff, I'd leave that up to the user. Reducing the matrix is one operation. Sorting its rows is completely separate. So I'd prefer to write code like:

 var reduced = reducedRowEchelonForm(matrix);
 var sorted = sortRows(reduced);

rather than passing more arguments to the reduce function. The sortRows() function is a nice, generic function (if I do say so myself) that doesn't care where its input came from. It just sorts it. If you don't want to sort it, well, don't.

---

OK, one quick look at the leadingPivotsToTop() function:

Basically, you can achieve this (I think) by

1. dividing the rows into those that have a pivot and the irrelevant ones
2. sorting the ones with a pivot using the existing sortRows
3. gluing (concatenating) the irrelevant ones back on

Something like this (untested!)

function leadingPivotsToTop(matrix) {
  var pivotIndices = [],
      pivotedRows = [],
      irrelevantRows = [];

  matrix.forEach(function (row) {
    var pivotIndex = row.findIndex(function (value) { return value === 1 });
    if(pivotIndex === -1 || pivotIndices[pivotIndex]) {
      irrelevantRows.push(row);
    } else {
      pivotIndices[pivotIndex] = true;
      pivotedRows.push(row);
    }
  });

  return sortRows(pivotedRows).concat(irrelevantRows);
}

Note that findIndex() isn't widely supported, but it's easy enough to roll your own - or just find a working one. The linked MDN page has an implementation you can borrow.