Current Index implementation

[rjbaucells]

Current Index implementation is missing the feature to allow the user to pick individual items from a Matrix and create a SubMatrix as a result. See the following example:

// 4x4 Matrix
A = [1, 2, 3, 4; 5, 6, 7, 8; 9, 10, 11, 12; 13, 14, 15, 16];

// B, same as A - using ranges
B = A[:,:];

// C, using a subset of A
C = A[1:2, 1:2]; // C = [1, 2; 5, 6]

// Using permutation vectors
p = [1, 2, 3];
q = [3, 2, 1];

// missing feature! (it should create a matrix using the specific indeces in the given order)
D = A[p, q]; // D = [3, 2, 1; 7, 6, 5; 11, 10, 9]

// use row 1, 2, 3, all columns
E = A[p,:]; // E = [1, 2, 3, 4; 5, 6, 7, 8; 9, 10, 11, 12];

// use columns 3, 2, 1, all rows
F = A[:,q]; // F = [3, 2, 1; 7, 6, 5; 11, 10, 9; 15, 14, 13];

The same syntax exists in Matlab and Scilab (replace square brackets by parentheses).

Implementing this feature in V2 will introduce a breaking change since creating an Index with one of the parameters being an Array will use the array items to create a Range (array code and matrix code).

The change I am proposing is:

• Allow Index to have a Range or a Set for each one of the dimensions
• Create a Set if one of the dimensions parameter is an Array (breaking change)
• Modify Index consumers to use the new feature (DenseMatrix, SparseMatrix)

This functionality is needed to correctly implement Matrix LU and QR factorization using full pivoting. Current partial pivoting (LU) will benefit form it as well since the row permutation is implemented by creating a Matrix out of the permutation vector and doing a Matrix x Matrix operation.

What do you think?

[josdejong]

This functionality has been there in the past, before we even had an Index type, see #59 (comment) and #61. We removed this feature then as the API was too complicated.

It's ok with me to reintroduce this functionality again, as long as we can come up with a clear API which works conveniently in both JavaScript as well as the expression parser.

If I understand you correctly you would like to replace the support for specifying a range as an Array [start, end, step] with an Array of arbitrary index values, a "set"? How to specify a range in the JS API then?

[rjbaucells]

Creating an Index with an Array parameter will no longer create a range for the given dimension, it will use the array as the actual indices for the dimension:

Actual behavior:

a = new Index([1, 2]); // dimension 0 of Index contains a Range '1:2'

New behavior:

a = new Index([1, 2]); // dimension 0 of Index contains a Set with the values [1, 2]

This allows the following (not possible today):

A = [.....]; // MxN matrix
p = [2, 7, 4, 5];
B = A[:,p]; // B contains the same number of rows as A but only columns 2, 7, 4 and 5. In that order

The same can be expressed in JavaScript as:

A = math.matrix([....]);  // MxN matrix
p = [2, 7, 4, 5];
B = A.subset(math.index(new Range(0, M, 1), p));

[josdejong]

Yes exactly. For the expression parser it will work like a charm. I'm worried though about the JavaScript notation A.subset(math.index(new Range(0, M, 1), p));. The notation is currently already quite complicated (see also #260) and would become even more complicated. And a pitfall here would be to call A.subset(math.index(math.range(0, M, 1), p));, which will expand the whole range into an array... We have to find a solution for that.

[rjbaucells]

Yes; calling math.range() will expand the whole range into an array, but the good part is that the result will be the same and no performance penalties.

math.range(0, M, 1) will return an array and the Index will use a Set for the given dimension.

[josdejong]

No performance penalties? I would expect that expanding a range into an Array would take more time/memory than just passing three values and running a simple for loop with that?

[rjbaucells]

It is a tradeoff, by having the Range as it is you need to run a range.forEach(), that will invoke a callback function with the actual index value, thing that is prohibited in SparseMatrix since seeking a value is too costly. Expanding the Range into an Array range.toArray() will give you a permutation vector that can be used in a fast algorithm in Sparse Matrices.

The Set contains the same API operations as a Range so the consuming code is not aware that the actual Index dimension is represented as a Set or a Range.

[josdejong]

I have to say it's hard for me to make a guess about the performance implications of just using math.range, I would have to try it and do some simple benchmarks. If there isn't a very large performance penalty it could be interesting to do this change in the API.

As for range.forEach(): nothing stops you from creating an inlined version of the range.forEach() method, you can just read the properties range.start, range.end, range.step from a Range object and execute a for or while loop :)

[rjbaucells]

Yes, doing the inline version will work but at the end I will be creating the same permutation vector (range.toArray()) since I cannot use an index seek in a SparseMatrix.

See the actual code in a proof of concept I am working on based on the proposed changes (not completed yet)