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.

• user1235831: Thank you for your review! One question, do you know if "map" method optimized by some JavaScript uglifier out there? If not then I am thinking about to write a parser for that purposes :-)
• Marc-Andre: @user1235831 You can see the reason why your edit was rejected in your user profile. It would be informative for you. You should have had a comment on this answer instead of editing directly. If you have a question of want more information you can visit the Code Review Meta.
• Flambino: @Marc-Andre Just saw the rejected edit(s). Thanks. (Turns out they were partly correct though)
• Flambino: @user1235831 What do you mean by "'map' method optimized by an uglifier"? Also: Saw your edits. You were right about the pivot.index/pivot.column confusion; fixed that. Thanks. But the row === sibling comparison remains unchanged; it works just fine, and is more direct than comparing the indices.
• Flambino: @user1235831 JS is object oriented, and we're looping through the same array, so row and sibling are just object references. So, at a certain index, row is the exact same object as sibling and vice-versa; it's two names pointing at the same thing. You could also compare the indices of course, but what really we want to know is whether it's the same row. So instead of asking "is this number the same as that number?" we just ask "is this the same thing?"