# Arrays proposal: array broadcasting, and adding new axes

**URL:** <https://racket.discourse.group/t/arrays-proposal-array-broadcasting-and-adding-new-axes/3948>\
**Category:** General\
**Tags:** arrays-proposal\
**Created:** [September 21, 2025, 12:50am UTC](https://racket.discourse.group/t/arrays-proposal-array-broadcasting-and-adding-new-axes/3948 "2025-09-21T00:50:17Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![Gambiteer](https://avatars.discourse-cdn.com/v4/letter/g/779978/32.png) [@Gambiteer](https://racket.discourse.group/u/Gambiteer)\
**Post date:** [September 21, 2025, 12:50am UTC](https://racket.discourse.group/t/arrays-proposal-array-broadcasting-and-adding-new-axes/3948/1 "2025-09-21T00:50:17Z")

</div>

I didn't know whether to start a new thread for each topic, or to stick with one thread. I started a new thread, but can do the other if more appropriate.

When composing this message, I didn't know how much background to give, so I probably erred on the long side.

First the background, then the proposal.

Brad

Background:

SRFI 231 does not have array broadcasting, either explicit or explicit, or a way to add new axes to an existing array.

Array broadcasting repeats an existing array or scalar as array "slices" along one or more new axes (or existing, length one, axes). Broadcasting is very common in array libraries---it's in Racket's math/array, Python's NumPy, Stan, Julia, Javascript stdlib, Matlab, etc.; links can be found here:

> **[Re: Array broadcasting | Simplelists](https://srfi-email.schemers.org/srfi-231/msg/32788225/)**
>
> Re: Array broadcasting

But I didn't see a need for broadcasting in SRFI 231 given that SRFI's other features.

For example, in Python broadcasting is used to implement both APL-style inner and outer products. (References are in the message linked above, and there are a few other comments in that thread.) SRFI 231 has APL-style array-inner-product and array-outer-product directly.

Or, broadcasting is used to simplify notation, to write `A x s` for an array `A` and a scalar `s` and have the value of `s` "broadcast" to an implicit array the same size and shape as `A` for component-wise multiplication.

Or, to use another common example, one could have a two-dimensional array `A` and a one-dimensional vector `v`, interpreted as a row vector, and write `A-v` to subtract `v` from each row of `A`.

But in SRFI 231 these examples are not so hard to do, because the SRFI has `array-map` (which deals directly with the first example) and `array-curry`, which decomposes an array by dimension; in the second example one would interpret the two-dimensional array `A` as a one-dimensional array of one-dimensional arrays, and then use array-map to subtract `v` from each of those one-dimensional rows (as arrays themselves). The array broadcasting approach would extend the vector `v` to a two-dimensional array of the right size and shape and the subtract that new array from `A` component-wise.

A further complication is that SRFI 231 arrays can have nonzero lower bounds in indices, while the other systems have (implicit) zero lower bounds in each dimension.

SRFI 231 does have a very general

```scheme
(array-share old-array new-domain mapping-from-new-domain-to-old-domain)

```

but the SRFI requires that the affine `mapping-from-new-domain-to-old-domain` from new multi-indices to existing multi-indices be one-to-one, so that mutating one array entry does not change a different array entry. To write broadcasting in this way necessarily requires that the mapping be many-to-one. (By the way, it appears that checking whether the mapping is 1-1 is equivalent to the small lattice vector problem, which I hear is NP-hard.)

Proposal:

In the mail thread starting here:

> **[First stab at API for adding new array axes and broadcasting arrays |...](https://srfi-email.schemers.org/srfi-231/msg/32833489/)**
>
> First stab at API for adding new array axes and broadcasting arrays

I laid out some of the issues I found problematic when considering integrating array broadcasting into SRFI 231 and then made the following initial proposal for new functionality.

I invite comments about this, or any other, aspect of the array SRFI.

==============================

Procedure: `(array<-scalar obj)`

Returns `(make-array (make-interval '#()) (lambda () obj))`. (People seem  
to like treating scalars like arrays, so let's make it easier to do so.)

Procedure: `(array-add-axes array vector)`

Sample calls:

```scheme
(array-add-axes A '#(* 0 1 * 2 * 3))
(array-add-axes A (vector '* 0 1 '* 2 '* 3))
(array-add-axes A '#(* * *)) ;; if A is zero-dimensional

```

Assumes that the first argument `A` is an array with dimension `d`.

Assumes the second argument is a vector `v` whose elements are the exact  
integers `0, 1, ..., d-1`, in order, interspersed with zero or more  
instances of the symbol `*`. [In a later message I note that there is no real reason to require that the integers be in order, in which case this routine would generalize `array-permute`.]

Returns a new array `A*` that shares its elements with `A` in the same  
lexicographical order.

Any entry in `v` with a nonnegative integer `j` entry denotes an axis in the  
domain of `A*` with the same lower and upper bounds as the `j`th axis of  
domain of `A`.

Any entry in `v` with the symbol `*` denotes a new axis in the domain of `A*`  
with lower bound 0 and upper bound 1.

[In a later message, I ask whether `array-squeeze`, which is a partial inverse to `array-add-axes` and is given in the SRFI document as an example, should be added to the SRFI if `array-add-axes` is added to the SRFI.]

Procedure: `(array-broadcast list-of-As)`

Assumes that `list-of-As` is a nonempty list of immutable arrays.

Let `D` be the maximum of the dimensions of the arrays in `list-of-As`.

A new list of arrays, `list-of-A*s`, is constructed by adding axes to the  
left of existing axes of the arrays in `list-of-As` with `array-add-axes` so  
that all arrays in `list-of-A*s` have dimension `D`.

The arrays in list-of-A\*s are now required to have the following  
properties, checked in order:

1. If the domain of any `A*` in `list-of-A*s` has an axis `j` with nonzero  
lower bound, then the domains of all arrays `A*` in `list-of-A*s` must have  
the same lower and upper bounds for axis `j`.

2. If the domain of any `A*` in `list-of-A*s` has an axis `j` with an  
upper-bound `u_j` different from 1, then axis `j` in the domain of all  
arrays `A*` in `list-of-A*s` must have either upper bound `u_j` or upper bound 1.

We then use the usual broadcast rules (which we have not yet specified  
here) to form a new `list-of-A **s` containing arrays `A** ` all with the same  
domain.

---

<div class="post-metadata">

**Author:** ![soegaard](https://yyz2.discourse-cdn.com/free1/user_avatar/racket.discourse.group/soegaard/32/19_2.png) [@soegaard](https://racket.discourse.group/u/soegaard)\
**Post date:** [September 21, 2025, 3:39pm UTC](https://racket.discourse.group/t/arrays-proposal-array-broadcasting-and-adding-new-axes/3948/2 "2025-09-21T15:39:06Z")

</div>

This is a detail, but I am curious about the choice of `array<-scalar` over `scalar->array`.  
Am I missing some convention?

---

<div class="post-metadata">

**Author:** ![Gambiteer](https://avatars.discourse-cdn.com/v4/letter/g/779978/32.png) [@Gambiteer](https://racket.discourse.group/u/Gambiteer)\
**Post date:** [September 21, 2025, 4:55pm UTC](https://racket.discourse.group/t/arrays-proposal-array-broadcasting-and-adding-new-axes/3948/3 "2025-09-21T16:55:52Z")

</div>

I can't recreate my thinking, it should really be `scalar->array`.

---

<div class="post-metadata">

**Author:** ![Gambiteer](https://avatars.discourse-cdn.com/v4/letter/g/779978/32.png) [@Gambiteer](https://racket.discourse.group/u/Gambiteer)\
**Post date:** [September 21, 2025, 9:55pm UTC](https://racket.discourse.group/t/arrays-proposal-array-broadcasting-and-adding-new-axes/3948/4 "2025-09-21T21:55:45Z")

</div>

Maybe `scalar->array` isn't needed, it's a one-line definition:

```scheme
> (define (scalar->array obj) (list*->array 0 obj))
> (array->list* (scalar->array 42))
42

```

I guess the question is how would people want to use it.

---

<div class="post-metadata">

**Author:** ![soegaard](https://yyz2.discourse-cdn.com/free1/user_avatar/racket.discourse.group/soegaard/32/19_2.png) [@soegaard](https://racket.discourse.group/u/soegaard)\
**Post date:** [September 22, 2025, 5:04pm UTC](https://racket.discourse.group/t/arrays-proposal-array-broadcasting-and-adding-new-axes/3948/5 "2025-09-22T17:04:18Z")

</div>

I would include `scalar->array` since the operations matches a conceptual operation.  
The other expression `(list*->array 0 obj)` requires switching to a more low-level way of thinking.

---

<div class="post-metadata">

**Author:** ![Gambiteer](https://avatars.discourse-cdn.com/v4/letter/g/779978/32.png) [@Gambiteer](https://racket.discourse.group/u/Gambiteer)\
**Post date:** [September 22, 2025, 5:22pm UTC](https://racket.discourse.group/t/arrays-proposal-array-broadcasting-and-adding-new-axes/3948/6 "2025-09-22T17:22:55Z")

</div>

Yes, good.

But the operation wraps any value, not just a scalar and possibly another array, in an array of dimension zero, as in

```scheme
heine:~/text/courses> /usr/local/racket/bin/racket
Welcome to Racket v8.16.0.2 [cs].
> (require math/array)
> (array-shape (array (array #[0 1 2])))
'#()

```

So maybe `object->array`? Just `->array`? (But the latter might be confusing, because people might think that `(->array '#(0 1 2))` would process the vector into an array with three elements.)

---

<div class="post-metadata">

**Author:** ![soegaard](https://yyz2.discourse-cdn.com/free1/user_avatar/racket.discourse.group/soegaard/32/19_2.png) [@soegaard](https://racket.discourse.group/u/soegaard)\
**Post date:** [September 26, 2025, 9:13am UTC](https://racket.discourse.group/t/arrays-proposal-array-broadcasting-and-adding-new-axes/3948/7 "2025-09-26T09:13:45Z")

</div>

> But the operation wraps any value, not just a scalar ...  
> You are right, I didn't further than scalars.

I think `object->array` is a fine choice.

---

<div class="post-metadata">

**Author:** ![Gambiteer](https://avatars.discourse-cdn.com/v4/letter/g/779978/32.png) [@Gambiteer](https://racket.discourse.group/u/Gambiteer)\
**Post date:** [October 1, 2025, 10:16pm UTC](https://racket.discourse.group/t/arrays-proposal-array-broadcasting-and-adding-new-axes/3948/8 "2025-10-01T22:16:17Z")

</div>

Should `object->array` have the calling sequence

```scheme
(object->array object [storage-class [ mutable?] ] )

```

like the other conversion routines, or just the simpler

```scheme
(object->array object)

```

which returns an immutable array with `generic-storage-class`?

The advantage of the second is that there is no error checking needed, no possible way to raise an exception, and I can't see the benefit of the first except for a certain uniformity of calling sequence.

Maybe I'm just tired, but I've been going back and forth on this while implementing it.

---

<div class="post-metadata">

**Author:** ![Gambiteer](https://avatars.discourse-cdn.com/v4/letter/g/779978/32.png) [@Gambiteer](https://racket.discourse.group/u/Gambiteer)\
**Post date:** [October 7, 2025, 12:55am UTC](https://racket.discourse.group/t/arrays-proposal-array-broadcasting-and-adding-new-axes/3948/9 "2025-10-07T00:55:44Z")

</div>

The most recent github commit [Add interval-insert-axis, compute-broadcast-array, object-\>array, arr… · gambiteer/srfi-231@90c70e0 · GitHub](https://github.com/gambiteer/srfi-231/commit/90c70e02dd975562c1c4af02e48c4af1486d9e97) adds `interval-insert-axis`, `array-insert-axis`, `compute-broadcast-array`, and `object->array`, all routines that are needed to implement general array broadcasting.

In the next commit, I will add a stripped-down version of the array library, with all the functionality but _no_ error checking and _no_ "optimizations". (I put "optimizations in quotes because it seems that sometimes the stripped-down library is faster than the general library.) The usual `generic-arrays.scm` has much code specialization for special cases and tries to catch and raise errors as soon as possible.

The idea behind `mini-arrays.scm` is to offer a second implementation, so that independent implementers can choose where on the axis between `mini-arrays.scm` and `generic-arrays.scm` they would like their implementations to sit.

There will also be a file `test-mini-arrays.scm` that is instrumented to include or elide tests to see that appropriate exceptions are raised, so it can be used to test both `generic-arrays.scm` and `mini-arrays.scm`.

---

<div class="post-metadata">

**Author:** ![Gambiteer](https://avatars.discourse-cdn.com/v4/letter/g/779978/32.png) [@Gambiteer](https://racket.discourse.group/u/Gambiteer)\
**Post date:** [October 11, 2025, 6:27pm UTC](https://racket.discourse.group/t/arrays-proposal-array-broadcasting-and-adding-new-axes/3948/10 "2025-10-11T18:27:54Z")

</div>

TL;DR: Should `array-broadcast` and `array-insert-axis` return only immutable array results?

NumPy has deprecated allowing a broadcasted array be mutable ("writable"), probably because, in general, changing one element of a broadcasted array changes other elements. See [numpy.broadcast\_arrays — NumPy v2.3 Manual](https://numpy.org/doc/stable/reference/generated/numpy.broadcast_arrays.html) . Strictly speaking, this restriction is not needed when the broadcast result has no more elements than the argument.

`math/array` has a similar routine, `array-broadcast`, [6.10&nbsp;Pointwise Operations](https://docs.racket-lang.org/math/array_pointwise.html#%28def._%28%28lib._math%2Farray..rkt%29._array-broadcast%29%29) , on which I modeled my own version.

I think NumPy's restriction is a good idea, so I thought I would have my `array-broadcast` return only immutable arrays.

NumPy also has the function `numpy.expand_dims`, which adds an axis of width 1 to its array argument. The result of `numpy.expand_dims` has the same elements as the argument, in the same order, so there is no need to restrict the result to be immutable, and indeed NumPy does not restrict the result to be immutable.

`math/array` has the similar `(array-axis-insert arr k [dk])`, which adds a new axis with an optional parameter that specifies an axis width greater than 1.

I modeled the my new `array-axis-insert` on `math/array`'s procedure.

The question is whether arrays returned by my `array-axis-insert` should be immutable, or should be immutable only when `dk>1`, to use `math/array`'s notation.

One could similarly ask whether the result of `array-broadcast` should be immutable only when there are more elements in the result than in the argument. Perhaps a similar question led to the following design decision in `math/array`: [6.10&nbsp;Pointwise Operations](https://docs.racket-lang.org/math/array_pointwise.html#%28def._%28%28lib._math%2Farray..rkt%29._array-broadcast%29%29)

> When [array-strictness](https://docs.racket-lang.org/math/array_strict.html#%28def._%28%28lib._math%2Farray..rkt%29._array-strictness%29%29) is #t, [array-broadcast](https://docs.racket-lang.org/math/array_pointwise.html#%28def._%28%28lib._math%2Farray..rkt%29._array-broadcast%29%29) returns a strict array when arr is nonstrict and the result has more elements than arr.

---

<div class="post-metadata">

**Author:** ![Gambiteer](https://avatars.discourse-cdn.com/v4/letter/g/779978/32.png) [@Gambiteer](https://racket.discourse.group/u/Gambiteer)\
**Post date:** [November 1, 2025, 2:45am UTC](https://racket.discourse.group/t/arrays-proposal-array-broadcasting-and-adding-new-axes/3948/11 "2025-11-01T02:45:46Z")

</div>

Here is an example of broadcasting that does not involve mutation and illustrates a situation that has been troubling me.

TL;DR: I think broadcasting, either explicit or implicit, should be applied only to specialized (strict) arrays, and the resulting arrays should be immutable.

You start with an `N x N` array of 3-vectors, a vector-field.

The task is to subtract the average of the 3-vectors from the `N x N` array of 3-vectors.

With broadcasting, one might write:

```scheme
(define (recenter A)
  
  ;; A is an array of dimension N x N x 3 
  ;; representing a two-dimensional N x N field of three-vectors.
  
  ;; We want to subtract the mean of the three-vectors from all the
  ;; entries of A, which we could write as follows.

  ;; The "problem" is that the array-map marked by
  ;; <<<<<<<<<<<<<<<
  ;; returns a generalized array (and I don't want to change that)
  ;; so every time an element of the average of the three-vectors is
  ;; accessed, the internal array-reduce is recomputed.

  (array-copy
   (array-map fl- A (array-broadcast
                     (array-map ;; <<<<<<<<<<<<<<<
                      ;; calculate the average of each coordinate
                      ;; of the three-vectors.
                      (lambda (slice)
                        (/ (array-reduce fl+ slice)
                           (interval-volume (array-domain slice))))
                      (array-curry
                       ;; write A as a three-dimensional array of NxN arrays
                       (array-permute
                        ;; put the coordinates of the three-vectors first
                        A '#(2 0 1))
                       2))
                     (array-domain A)))))

```

Here I've explicitly called `array-broadcast`, but with implicit broadcasting that explicit call would be elided.

I now have an efficient implementation of array broadcasting for specialized (strict) arrays, but not for generalized (nonstrict) arrays (which is not so easy, in general). So here are two routines with as efficient implementations of broadcasting generalized and specialized arrays as I can imagine:

```scheme
(define (recenter-generalized A)
  (array-copy
   (array-map - A (let* ((average
                          ;; the generalized array result is not copied to a specialized array
                          (array-map
                           (lambda (slice)
                             (/ (array-reduce fl+ slice)
                                (interval-volume (array-domain slice))))
                           (array-curry
                            (array-permute
                             A '#(2 0 1))
                            2)))
                         (average-getter ;; use internal unsafe getter
                          (%%array-unsafe-getter average)))
                    (make-array (array-domain A)
                                (lambda (i j k)
                                  (average-getter k)))))))

(define (recenter-specialized A)
  (array-copy
   (array-map - A (array-broadcast
                   ;; copy the array to a specialized array first
                   (array-copy
                    (array-map
                     (lambda (slice)
                       (/ (array-reduce fl+ slice)
                          (interval-volume (array-domain slice))))
                     (array-curry
                      (array-permute
                       A '#(2 0 1))
                      2)))
                   (array-domain A)))))

```

And we can test it on very small and not quite so small arrays:

```scheme
(define A (list->array (make-interval '#(10 10 3)) (map inexact (iota 300))))
(define B (time (recenter-specialized A)))
(define C (time (recenter-generalized A)))

(define A* (list->array (make-interval '#(100 100 3)) (map inexact (iota 30000))))
(define B* (time (recenter-specialized A*)))
(define C* (time (recenter-generalized A*)))

```

With a `10 x 10` array of 3-vectors, the timings are not sooooo different:

```scheme
(time (recenter-specialized A))
    0.000084 secs real time
    0.000085 secs cpu time (0.000047 user, 0.000038 system)
    no collections
    37376 bytes allocated
    12 minor faults
    no major faults
    292467 cpu cycles
(time (recenter-generalized A))
    0.000894 secs real time
    0.000894 secs cpu time (0.000666 user, 0.000228 system)
    no collections
    574640 bytes allocated
    71 minor faults
    no major faults
    3202257 cpu cycles

```

But with a `100 x 100` array of 3-vectors, the timings differ dramatically:

```scheme
(time (recenter-specialized A*))
    0.001860 secs real time
    0.001860 secs cpu time (0.001860 user, 0.000000 system)
    no collections
    783088 bytes allocated
    84 minor faults
    no major faults
    6671997 cpu cycles
(time (recenter-generalized A*))
    6.449153 secs real time
    6.433722 secs cpu time (6.428914 user, 0.004808 system)
    no collections
    55809952 bytes allocated
    6907 minor faults
    no major faults
    23162995368 cpu cycles

```

A library should not allow the programmer to fall inadvertently into such a hole,, so I conclude that generalized (nonstrict) arrays should not be broadcast.

But if you have a call

```scheme
(array-map f A B)

```

and `B` is supposed to be broadcast to the domain (shape) of `A`, then if `B` is generalized one has two choices: raise an error and say that `B` is not specialized; or silently copy `B` to a specialized `B*`, and then broadcast `B*`.

The thing is, in all other circumstances, the array elements of the arguments to `array-map` are not evaluated _by_ `array-map`---`array-map` is supposed to set up work to do be done _later_ by routines like `array-fold` or `array-every`. Copying `B` to `B*` would require `array-map` to evaluate `B`'s getter on its domain, which would break this convention of the library. I conclude that an error should be raised, that `B` should not be implicitly copied to `B*` by `array-map`.

Brad

---

<div class="post-metadata">

**Author:** ![Gambiteer](https://avatars.discourse-cdn.com/v4/letter/g/779978/32.png) [@Gambiteer](https://racket.discourse.group/u/Gambiteer)\
**Post date:** [December 28, 2025, 10:30pm UTC](https://racket.discourse.group/t/arrays-proposal-array-broadcasting-and-adding-new-axes/3948/12 "2025-12-28T22:30:11Z")

</div>

I've now added an `array-broadcast` procedure to the follow-on to SRFI 231: [Add array-broadcast · gambiteer/srfi-231@5ec0ad1 · GitHub](https://github.com/gambiteer/srfi-231/commit/5ec0ad1ac34378e99f060a68e6fc3c01604f7823)

The commit message details some of the decisions that were made after studying Python's NumPy library and Racket's `math/array`:

> [@](#):
>
> We implement array broadcasting in the style of Python's NumPy or Racket's math/array.
> 
> Arrays in both those systems have only zero lower bounds on indices, and while Racket has what we call generalized arrays, NumPy does not. Additionally, Racket allows R-style "permissive" broadcasting.
> 
> So we had to make some decisions about how to broadcast arrays with nonzero lower bounds and generalized arrays.
> 
> We decided:
> 
> 1. Axes initially added at the left of an array's domain have lower bounds 0 and upper bounds 1.
> 
> 2. Axes with nonzero lower bounds cannot be expanded in broadcasting.
> 
> 3. Generalized arrays can be broadcast only to new intervals with the same volume (so, effectively, can only have length-1 axes added on the left), returning a generalized array. I want only 1-1 mappings of indices for generalized arrays to avoid surprises when those broadcast arrays are "consumed" by some later operation.
> 
> 4. We don't allow R-style "permissive" broadcasting, which requires a non-affine mapping of indices that is not composable with other array transformations.
> 
> 5. If an array is broadcast to an interval with more elements (of greater volume), then the resulting array is set to be immutable.
> 
> There is not yet automatic broadcasting of array arguments, only the new procedure array-broadcast.
> 
> There is no new documentation.

Any comments about these choices are welcome.

---

<div class="post-metadata">

**Author:** ![Gambiteer](https://avatars.discourse-cdn.com/v4/letter/g/779978/32.png) [@Gambiteer](https://racket.discourse.group/u/Gambiteer)\
**Post date:** [January 4, 2026, 7:25pm UTC](https://racket.discourse.group/t/arrays-proposal-array-broadcasting-and-adding-new-axes/3948/13 "2026-01-04T19:25:23Z")

</div>

With this commit the followup to SRFI 231,

> <https://github.com/gambiteer/srfi-231/commit/36c7522136a06f6b58b664366084b0cc49f18186>
>
> We add implicit broadcasting of array arguments, in the style of Python's NumPy …and Racket's math/array. We add the parameter array-broadcasting? that controls whether broadcasting of arguments is allowed.
> 
> At the moment we have decided not to broadcast array arguments of array-stack and array-append.
> 
> 231-bis.sld:
> 
> 1. Add array-broadcasting?
> 
> generic-arrays:
> 
> 1 Add parameter array-broadcasting?, initial value #t.
> 
> 2. Add %%broadcast-array-arguments, use it in array-map, array-for-each, array-every, array-any, array-fold-left, and array-fold-right.
> 
> mini-arrays:
> 
> 1 Add parameter array-broadcasting?, initial value #t.
> 
> 2. Add %%broadcast-array-arguments, use it in array-map, array-for-each, array-every, array-any, array-fold-left, and array-fold-right.
> 
> test-arrays.scm:
> 
> 1. Add minimal testing for array-broadcasting in array-{map|for-each|every|any|fold-left|fold-right}. Remove tests that are no longer appropriate.
> 
> test-mini-arrays.scm:
> 
> 1. Add minimal testing for array-broadcasting in array-{map|for-each|every|any|fold-left|fold-right}. Remove tests that are no longer appropriate.

I've now added implicit array broadcasting to array arguments of

1. `array-map`
2. `array-for-each`
3. `array-every`
4. `array-any`
5. `array-fold-left`
6. `array-fold-right`

I've followed Python's NumPY library, and not Racket's `math/array`, in not implicitly broadcasting array arguments to `array-stack` and `array-append`.

Implicit broadcasting of array arguments is controlled by a parameter, `array-broadcasting?`, as in Racket's `math/array`.

Comments are welcome.