TL;DR: I'm considering turning the procedures
array-{inner|outer}-product into procedures that return a specialized
array.  I'm also tempted to change the calling sequence for
array-inner-product to something that can be implemented more easily.

There is a discussion thread comparing some aspects of Racket's
math/array and math/matrix libraries and SRFI 231 at

https://racket.discourse.group/t/requesting-help-to-translate-small-array-program-into-typed-racket/4079

The part of the conversation perhaps most relevant to this list is
==============================
Studying this code leads me to reconsider the specifications of
array-outer-product and array-inner-product in SRFI 231, which now seem
a bit sloppy. Some things perhaps to copy from how matrix* is coded:

1. If necessary, copy argument arrays to specialized (strict) arrays,
after permuting the second array argument to place its first axis last.

2. Return a specialized (strict) array.

3. Right now, the calling sequence is

(array-inner-product A f g B)

and the inner loop is computed by (a and b are one-dimensional subarrays
of A and B):

(lambda (a b)
   (array-reduce f (array-map g a b)))

where, for the matrix* loop, g is *, and f is +.

Perhaps the calling sequence should be

(array-inner-product A f-of-g g B)

where for matrix*, for example, g would be * and f-of-g would be

(lambda (sum a b) (+ sum (* a b)))

so it would be not two but one function call (other than the + and *,
which can be inlined by the compiler, as in Gambit) per iteration.

I'll have to think about things more.
==============================