Feature: add/subtract every row of matrix by vector

[JasonShin]

I am getting a dimension error while using math.subtract()

const X = [[1, 2], [3, 4], [5, 6]]
const A = math.transpose(X);
// A = [ [ 1, 3, 5 ], [ 2, 4, 6 ] ] 
const M = math.mean(A, 1);
// M = [ 3, 4 ]
const C = math.subtract(X, M);

Running the code raises the below error

....mathjs/lib/function/arithmetic/subtract.js:172
    throw new DimensionError(xsize.length, ysize.length);

I believe math.subtract(X, M) should return

[[-2., -2.],
[ 0.,  0.],
[ 2.,  2.]]

same goes for math.sum math.multiply

[JasonShin]

In the meantime, I will write an extension method to achieve my goal

Temporary extension method that I am using:

/**
 * Subtracts two matrices
 * @param X
 * @param y
 */
const subtract = (X, y) => {
  const _X = _.clone(X);
  for (let rowIndex = 0; rowIndex < _X.length; rowIndex++) {
    const row = X[rowIndex];
    for (let colIndex = 0; colIndex < row.length; colIndex++) {
      const column = row[colIndex];
      // Supports y.length === 1 or y.length === row.length
      if (y.length === 1) {
        const subs = y[0];
        _X[rowIndex][colIndex] = column - subs;
      } else if (y.length === row.length) {
        const subs = y[colIndex];
        _X[rowIndex][colIndex] = column - subs;
      } else {
        throw Error(`Dimension of y ${y.length} and row ${row.length} are not compatible`);
      }
    }
  }
  return _X;
}

[josdejong]

Yes, that can be a useful extension. Right now math.subtract, math.add, etc simply do element-wise subtraction/addition. Allowing to subtract a vector row-wise from a matrix is something you can do in Matlab and Octave.

I will change the title of this topic a bit.

Anyone interested in implementing this feature?

[dvd101x]

An explanation of how this could work inside math.evaluate( )

A = [1, 2; 3, 4; 5, 6];
M = mean(A, 1);
M = [M, M, M];      # broadcast would do something like this
C = A - M           # yields [[-2, -2], [0, 0], [2, 2]]

And once broadcasting is implemented the original expression could be reduced like this.

A = [1, 2; 3, 4; 5, 6];
C = A - mean(A, 1)   # could yield [[-2, -2], [0, 0], [2, 2]]

[dvd101x]

I have an example of the algorithm working like so:

A = [1, 2; 3, 4; 5, 6];
C = A - mean(A, 1)   # yields [[-2, -2], [0, 0], [2, 2]]

To validate the behavior, please copy paste the the above example into the left pane of this project
https://dvd101x.github.io/web-thermodynamics/

I included broadcasting between scalars, vectors and N dimensional matrixes to math.add( ) , math.subtract( ), math.dotMultiply( ) and math.dotDivide( ).

I will be looking into ways of including math.typed( ) instead of some functions I made to check the types and try to make it easier to read.

Note: I edited some errors in my comments regarding divide, when I should have used dotDivide

Off topic: I will be looking into other functions that could benefit from elementwise operations with two or more inputs like log(N, base) that might be nice to accept broadcasting.

[josdejong]

1. broadcast_to(A, dimensions)

So it is a sort of a repeat function. I guess that can be useful to create a filled matrix with a certain repeated pattern.

2. broadcast_matrices(A, B)

Ahh, now I get how you can use broadcasting to make this work. It is a very neat solution. Technically, I think it is not the fastest, temporarily creating a large broadcasted matrix. But we can start with this, and if needed, we can always switch to a faster solution similar to the reduce that we have for mean and other functions.

This functionality will be very useful, I'm looking forward to have this working!

We can go through the list with functions so see which functions can benefit from this: search for all text occurrences of (x, y) and (a, b) on the page 😁

3. add(a, b, dimension)

One question: suppose we have a square matrix, and we're going to broadcast a vector. Which dimension should we take, and should we make this customizable, like with mean you can provide the dimension as second argument? I.e. should we implement add(a, b, dimension)?

[dvd101x]

Thank you for reviewing this, I hope you don't mind but here is an overview of the algorithm.

broadcast sizes

Given a set of matrix sizes, find the broadcasted size that can fit all dimensions for all matrices. All matrices would need to check if they need to grow, for that it's need to increase size empty or 1 to the maximum of all sizes in that dimension.

matrix	$dim_1$	$dim_2$ (rows)	$dim_3$ (columns)
A	2	3	4
B		3	1
C		1	4
broadcasted size	2	3	4

Thus is only possible to broadcast a dimension that is empty or 1 to a higher size.

matrix	$dim_1$	$dim_2$ (rows)	$dim_3$ (columns)
A	2	3	4
B		3	1
C	1	2	4
broadcasted size	2	fail	4

error: Failed because it's not possible to broadcast shape 2 to shape 3. For this to work that dimension should be empty or 1

The algorithm must work for N dimensions and M matrices

1. broadcast_to(A, dimensions)

Yes, it first checks if it's "broadcastable" and then repeats the numbers, the issue here is that it doesn't repeat a patter as it's not supposed to broadcast an existing dimension that is not of size one.

Personally I think it might be useful for something like [1, 2] + [3, 4, 5, 6] as [1, 2, 1, 2] + [3, 4, 5, 6] and which could yield [4, 6, 6, 8] but the typical convention is that it doesn't work like that and should throw an error of not matching dimensions so I kept it that way.

For something like patterns we can use orthogonal dimensions and then reshaping back

A = [1, 2] + [[3,4],[5,6]]   # [[4, 6], [6, 8]]
reshape(A, [-1])             # [4, 6, 6, 8]

2. broadcast_matrices(A, B)

I agree it's not very efficient in that sense, it was easier for me to lean on the capabilities to do operations on matrices of the same size. Maybe a better way would be to modify the code that does operations on matrices of the same size to first check for sizes, calculate only the broadcasted size for all matrices and if it doesn't throw an error then use the modulus operator to find the correct indexes of the matrices, like math.mod(indexPositionInThisDimension, broadcastedSizeInThisDimension) in the case of 0-indexing in a way to optimize for memory, maybe affecting a bit of performance. Maybe just checking if the broadcasted index in the dimension is greater than 1 on the matrix then use it, otherwise use index 0 depending on what is more computationally expensive, a mod or an if. I haven't thought about this thoroughly.

I can check on the list of functions for candidates.

3. add(a, b, dimension)

I think there is no need to specify a dimension as it takes the already available dimensions for all matrices and the matrices can be reshaped beforehand.

In the example of a square matrix and a vector

The dimension that gets repeated is the first one (counting on the final broadcasted size) and one that does not exist yet on the vector.

matrix	$dim_1$	$dim_2$
A	2	2
B		2
broadcasted size	2	2

If we want to repeat on the second one we could first reshape(B, [-1, 1])

matrix	$dim_1$	$dim_2$
A	2	2
B	2	1
broadcasted size	2	2

[dvd101x]

Hi, I found these candidates for broadcasting:

• add
• dotDivide
• dotMultiply
• dotPow
• gcd
• lcm
• log: maybe
• floor
• round
• subtract
• bitAnd
• bitOr
• bitXor
• and
• or
• xor
• compare
• compareText
• deepEqual
• equal
• equalText
• larger
• largerEq
• smaller
• smallerEq
• unequal
• atan2: it has not lost it's elemen twise operations
• to

Method: for every function with multiple inputs, check if the second argument could be a matrix and operations are element wise

[JasonShin]

Once I finish the getMatrixDataType pull request, I am happy to implement it =). Sorry I was busy for the last few days =)...

[josdejong]

👍 thanks for the offer. and for keeping focus :)

[josdejong]

We can consider two cases in this regard, as noted by @alvinhui in #1228:

Case 1:

const a = [[1, 2, 3], [4, 5, 6], [7, 8, 9]];
const b= [[9], [8], [7]]; 

math.subtract(a, b); // should be [[-8, -7, -6],  [-4, -3, -2], [ 0,  1,  2]]

Case 2

const a = [[1, 2, 3], [4, 5, 6], [7, 8, 9]];
const b= [9, 8, 7];

math.subtract(a, b); // should be [[-8, -6, -4], [-5, -3, -1], [-2,  0,  2]]

This behavior is similar to subtraction between Python arrays:

a = np.array([[1,2,3], [4,5,6], [7,8,9]])
b = np.array([9,8,7])
c = [[9], [8], [7]]; 

a - b
a - c

[dvd101x]

Thanks. I have it in my repository is it ok if I make a pull request?

[josdejong]

Yes please!

[dvd101x]

I included pull request #2895 regarding this topic.

[josdejong]

Published now in v11.6.0 🎉

[dvd101x]

Thanks!

[ChristopherChudzicki]

This could definitely be useful in some cases. But I suggest that there should be a way to turn it off.

Example: I mainly use MathJS for math webapps to be used by high-school and undergraduate students. I generally want students to think that adding matrix + vector is not valid. Great, MathJS throws a useful error!

I could be one of a small minority, but I would guess that MathJS's current behavior is useful to other people.

(Admittedly, if people like me want this error, it's not too hard to do our own validation. But that makes upgrading harder.)

Incidentally, numpy calls this broadcasting. (Although, Numpy's broadcasting is more general.)

Edit: To be clear, doing my own validation really is not a big deal. Just something to consider :)

[JasonShin]

Thanks for the inputs guys =)

[josdejong]

Thanks for sharing your thoughts @ChristopherChudzicki. Very good point to weight convenience against more explicit notations preventing you from accidental mistakes. We could also think about a solution where we have separate, explicit functions (like subtractVector(matrix, vector)) to do operations on matrices and vectors, or pass an extra option to these functions
(like subtract(matrix, vector, { broadcast: true })
or subtract(matrix, vector, { rowWise: true })).

Just thinking aloud here :)

[ChristopherChudzicki]

That said, it occurs to me that allowing matrix + row or matrix + vector consistent with how MathJS handles matrix + scalar. What I mean is, in the matrix+scalar case MathJS already does not demand strict shape-matching.

So maybe this is the "right" behavior.

[josdejong]

Yes that's true.

I'm just exploring different options, not sure yet what solution would best serve the audience of mathjs.

[johnathanalyst]

Any progress on this? Am getting into mathjs as a numpy replacement in the browser for implementing ML/AI stuffs. Broadcasting is a verrry handy feature that makes a lot of learning algorithms a bit simpler to implement. Happy to pitch in.

[johnathanalyst]

Also might add the word "Broadcasting" to the title for the numpy folks looking for help on this. I'd imagine they'll become a larger and larger part of your target audience, if they aren't currently.

[josdejong]

@chivington thanks for your input. There is no progress on this feature so far. I've marked the issue with a label "help wanted", now.

Anyone interested in implementing this feature?

[joeGermany]

Hello!! My case is similar to @chivington. I am getting into machine learning on the web and I need an alternative for NumPy. The above feature (broadcasting) is a real help for easier manipulation of matrices. Is there any way I can implement it alone in the meantime until it is a mainstream MathJS function (if it will ever be). 😃

[josdejong]

Is there any way I can implement it alone in the meantime until it is a mainstream MathJS function (if it will ever be).

The solution would be to implement your own version of all operators/functions like subtract, add, etc, like Jason suggests in #1165 (comment).