Operations on jagged arrays

[dvd101x]

Currently mathjs has support for mapping jagged arrays (non rectangular arrays).

config({matrix:"Array"})
a = [[1,2], 3]

map(a, _(x) = x+1)
# [[2, 3], 4]

But operations are not possible,

a + 1
# DimensionError: Dimension mismatch (1 < 2)

The broadcasting algorithm could be tweaked to allow jagged arrays, something like this.

config({matrix:"Array"})
a = [[1,2], 3]

a + 1
# [[2, 3], 4]

a + [[1],[2]]
# [[2, 3], 5]

Would this functionality be of interest?

Reference:
https://awkward-array.org/doc/main/index.html

[josdejong]

Interesting idea. Some first thoughts:

• Are there actual use case for this?
• I can understand [[2, 3], 4] + 1, but I think it can become hard to understand what will happen in cases like [[1,2], 3] + [[1],[2]]. I'm not sure if we should want to facilitate such cases.
• How complex would it be to implement something like this?

[dvd101x]

Sorry for the multiple posts, there is no rush.

Here is a more elaborated example mimicking some current functions.

const arr1 = [[1, 2], 3];
const arr2 = [1, 2];

console.log("add", add(arr1, arr2));
console.log("subtract", subtract(arr1, arr2));
console.log("dotMultiply", dotMultiply(arr1, arr2));
console.log("dotDivide", dotDivide(arr1, arr2));
console.log("sum", sum(arr1));
forEach((a, b) => console.log(a, b), arr1, arr2);

function add(...args) {
  return map((...X) => X.reduce((acc, val) => acc + val, 0), ...args);
}

function subtract(...args) {
  return map((...X) => X.reduce((acc, val) => acc - val), ...args);
}

function dotMultiply(...args) {
  return map((...X) => X.reduce((acc, val) => acc * val, 1), ...args);
}

function dotDivide(...args) {
  return map((...X) => X.reduce((acc, val) => acc / val), ...args);
}

function sum(array) {
  let total = 0;
  forEach((x) => total += x, array);
  return total;
}

function map(callback, ...args) {
  return broadcast(callback, args);
}

function forEach(callback, ...args) {
  broadcast(callback, args);
}

function broadcast(callback, args) {
  if (args.length === 0) return;
  if (args.length === 1) return deepMap(callback, args);
  if (args.length === 2) return broadcastTwo(callback, args[0], args[1]);
  if (args.length > 2) return broadcastMultiple(callback, args);
}

function deepMap(callback, array) {
  if (Array.isArray(array)) {
    return array.map((item) => deepMap(callback, item));
  }
  return callback(array);
}

function broadcastTwo(callback, arg1, arg2) {
  const arg1IsArray = Array.isArray(arg1);
  const arg2IsArray = Array.isArray(arg2);
  const hasArrays = arg1IsArray || arg2IsArray;
  if (hasArrays) {
    const length1 = arg1IsArray ? arg1.length : null;
    const length2 = arg2IsArray ? arg2.length : null;
    const maxLength = Math.max(length1, length2);
    if (length1 > 1 && length1 < maxLength)
      throw `Not possible to broadcast item 1 with length ${length1} to length ${maxLength}`;
    if (length2 > 1 && length2 < maxLength)
      throw `Not possible to broadcast item 2 with length ${length2} to length ${maxLength}`;
    const result = Array(maxLength);
    for (let i = 0; i < maxLength; ++i) {
      result[i] = broadcastTwo(
        callback,
        length1 > 1 ? arg1[i] : length1 === 1 ? arg1[0] : arg1,
        length2 > 1 ? arg2[i] : length2 === 1 ? arg2[0] : arg2
      );
    }
    return result;
  }
  return callback(arg1, arg2);
}

function broadcastMultiple(callback, args) {
  const argIsArray = args.map(Array.isArray);
  const hasArrays = argIsArray.some((x) => x);
  if (hasArrays) {
    const lengths = args.map((arg, argIndex) =>
      argIsArray[argIndex] ? arg.length : null
    );
    const maxLength = Math.max(...lengths);
    lengths.forEach((length, i) => {
      if (length > 1 && length < maxLength)
        throw `Not possible to broadcast item ${i} with length ${length} to length ${maxLength}`;
    });
    const result = Array(maxLength);
    for (let i = 0; i < maxLength; ++i) {
      result[i] = broadcastMultiple(callback,
        args.map((arg, argIndex) =>
          lengths[argIndex] > 1
            ? arg[i]
            : lengths[argIndex] === 1
              ? arg[0]
              : arg
        )
      );
    }
    return result;
  }
  return callback(...args);
}

Performance wise the current algorithms are faster when comparing them with the ragged version working on rectangular arrays of the same size. But this proposal is faster when broadcasting is needed, thus maybe a good strategy is to assume arrays are rectangular and if they fail then try the jagged version. Also there might be a benefit of implementing a broadcasting algorithm for rectangular arrays (that doesn't copy in memory as currently happens in matrixAlgorithmSuite).

jagged arrays performance jsbench

[josdejong]

It's interesting to me that in mathjs is possible to map, flatten, sum, max, min, std etc. jagged arrays (not possible in numpy) but it's not possible to do operations like add.

Yeah, good point.

You're explanation of this [[1, 2], 3] + [1, 2] example makes sense indeed.

Looking at your code examples it looks like it's not very complex, so it can be interesting to look into this. How would this be integrated in the existing functionality? I think this would (partially?) replace the matrixAlgorithmSuite methods that implement element-wise operations?

[dvd101x]

As you mentioned, I think adding the feature would be a modification of matrixAlgorithmSuite when it has Arrays don't convert them to Matrices. But attempt something like the proposed algorithm.

A challenge is that depending on the information of the inputs, it's possible to use a more optimized algorithm with less uses of Array.isArray and array.length.

type	example	check if array per dimension	check array length per dimension
jagged heterogeneous	[[1, 2], 3]	all elements	all elements
jagged homogeneous	[[1, 2], [3]]	only the first	all elements
rectangular	[[1, 2], [3, 3]]	only the first	only the first
matrix	[[1, 2], [3, 3]] size [2, 2]	known	known

Something interesting is that broadcasting arrays together will create a new array that has the most information of all the inputs. For example jagged heterogeneous + jagged homogeneous + rectangular = rectangular.

Maybe this is thinking ahead, and the first step could be only to broadcast two arrays assuming jagged heterogeneous and return a new one, implement that into matrixAlgorithmSuite and see if it works.

[dvd101x]

I'm interested in working on a proof of concept that for matrixAlgorithmSuite does the following

• doesn't apply broadcastTo to both arrays or matrices, but applies broadcast directly (for rectangular arrays of known sizes)
• doesn't apply broadcastTo to arrays and then convert array to matrices, but applies broadcast directly (for non homogeneous jagged arrays).

Even if the jagged functionality is not needed in the end, my guess is that the change for rectangular arrays of known sizes will be faster.

[josdejong]

Sounds good to create a minimal proof of concept to get a concreate idea on what this involves and now complex it will be 👍

[gwhitney]

To chime in along the lines of Jos's first comment, I am concerned about an example like [[1, 2], [3, 4]] + [[1], [3, 4]] You reference a python awkward arrays package, and skimming its docs it seems the second summand would be regarded as missing its 0,1 entry and the sum would be [[2,2], [6,8]] (not 100% certain, didn't run the code, but just saying this would be a completely reasonable interpretation of how to add jagged arrays). On the other hand, I think you are suggesting it be [[2, 3], [6, 8]] from the perspective of broadcasting.

My takeaways:

• I think trying to define operations on jagged arrays is going to be very prone to ambiguous cases with a high likelihood of surprising some users, so we need to tread lightly and perhaps think about throwing exceptions in the more ambiguous cases, to avoid user surprises.
• there may be an inherent conflict between the perspectives of "wanting to broadcast" and "treating jagged arrays as missing data" that might mean, for example, that these operations should only broadcast scalars to corresponding array parts, not portions of extent 1 to corresponding array parts.

Thanks for considering.

[dvd101x]

Hi Glen, thanks for your consideration.

In your example awkward throws an error

import awkward as ak
ak.Array([[1, 2], [3, 4]]) + ak.Array([[1], [3, 4]])
# ValueError: cannot broadcast nested list

This is my misunderstanding. I think the package is more inline with your description, as replacing the portion of extend 1 to a scalar does run a result, broadcasting the scalar (not filling the blanks)

import awkward as ak
ak.Array([[1, 2], [3, 4]]) + ak.Array([1, [3, 4]])
# [[2, 3],
#  [6, 8]]

So I think Jos's first comment, your latest comment and the awkward array package have a shared idea that these operations would recursively apply broadcasting rules but only to broadcast scalars.

A downside that I've just found out is that if the arrays are rectangular it doesn't broadcast due to the same limitation.

import awkward as ak
ak.Array([[1, 2], [3, 4]]) + ak.Array([[1], [2]])
# ValueError: cannot broadcast nested list

But it would simply work for numpy arrays (rectangular)

import numpy as np
np.array([[1, 2], [3, 4]]) + np.array([[1], [2]])
# array([[2, 3],
#        [5, 6]])

It's a limitation that to be honest is not clear to me, maybe there is a mathematical basis for this.

The proposed algorithm would be a superset as it would work the same for rectangular arrays and jagged heterogeneous arrays without knowing beforehand. At the moment I haven't started with the proof of concept in the codebase but I think it could work either way:

• recursively applying broadcasting rules
• recursively applying broadcasting rules but:
  • if there is a mix of scalars and arrays in the N dimension just broadcast the scalars
  • if there are only arrays only broadcast if all arrays have the same length

Something along those lines.

Appendix

Another mistake I made is that I was saying that

a = [[1,2], 3]
a + [[1],[2]]
# [[2, 3], 5]

But I should have said either:

a = [[1,2], 3]
a + [[1],[2]]
# [[2, 3], [5]]

or

a = [[1,2], 3]
a + [[1], 2]
# [[2, 3], 5]

[gwhitney]

Something along those lines.

Yes, well, my recommendation is just to try to keep it as simple and intuitive as possible and in doubtful/ambiguous cases err on the side of throwing. Otherwise i feel the facility could become too complicated to really be useful. That's all.

[josdejong]

Makes sense. It should indeed be easy to understand what the output will be. And implementing it should not result in a lot of added complexity.

[dvd101x]

Thanks,

While working on it I found a fundamental flaw in the previous algorithms on this thread. I will fix in the example.

[dvd101x]

Hi, just to let you know I was having a hard time with the algorithm of broadcasting jagged heterogeneous arrays correctly. The main difficulty is that my initial examples where not aligning the sizes correctly as it needs to know at least the maximum depth beforehand. Even though I think it's possible, it might not be as easy.

A simpler version with just operations between a jagged array and a scalar is provided at #3607

I'm still working on the algorithm to see if it's possible to do it in a simple way but working on it sporadically.

[dvd101x]

My proposed algorithm was broadcasting from the left (which makes sense when you don't know the number of dimensions) but rectangular arrays are broadcasted from the right (the number of dimensions is known). That's the fundamental flaw I was referring to.

It happens that awkward arrays are broadcasted from the left. Now it's starting to make sense.

So one possible algorithm could be.

• if A and B are rectangular then use the current broadcasting algorithm that depends on the size.
• else use the jagged broadcasting algorithm (because at least one size is not known)

[gwhitney]

One other important point: if there is a PR along these lines, it should establish the behavior for jagged arrays with empty sections, even for functions like map and filter and foreach that currently have some behavior on jagged arrays. The submitter of PR #3567 found for example that [[], [1]] and [[1],[]] behave differently in such contexts, which should be corrected.