Summary: The right choice for what `operator*` should do depends
on how many dimensions the language wants to elegantly handle.
Matlab and Eigen have chosen the path of notation elegance.
Numpy has chosen the path of higher dimensions.
Using operator* for matrix multiplication is concise and aligns
with how linear algebra is taught and thought. i.e. using vectors
and matrices.
Eigen, which deals exclusively with rows and columns[1], feels
quite natural for linear algebra and optimization problems. It
matches how most of us (except Einstein) think about the
operations involved.
Using higher dimensions elegantly comes at a cost of making some
linear algebra operations (specifically matrix multiplication)
more awkward to describe and implement.
If a library wants to handle the regimes {dims<=2, dims>2},
it make little sense to use operator* for matrix multiplication.
The meaning is unclear in higher dimensions.
There are, after all, only two sides of a '*' character: left and
right.
In numpy and scipy, there is a conceptual divide between the data
(numpy arrays) and the operations (e.g. scipy.linalg, numpy
matrices[3]).
This is easy to forget since the linear algebra functions [4]
just work as expected, assuming you send in 1d and 2d data.
Unfortunately, this conceptual divide makes it difficult to define
`operator*` as matrix multiplication. For example, if `a` and `b`
are 1d arrays of the same length, should `a*b` be an inner product
or an outer product?
Numpy arrays allow one to store data on meaningful hypercube axes
and perform inner products and summations over those axes.
In a cell tower simulation[2], this helped me consider channel
responses for different users, antennas, base stations, paths, and
multiple simulation realizations as axes on a hypercube.
This would also be possible in Matlab, but it would be clunky:
lots of `transpose`, `reshape`, and `bsxfun`.
Some of you might be thinking that Matlab handles higher
dimensions just fine -- try adding a singleton axis at the end:
`size( reshape( 1:10 ,5,2,1 ) )`.
I tried to be fairly objective and unbiased in my discussion
above, but I do have a preference, of course.
That preference is based on how *I* use my tools in the years
I've spent with Matlab, Numpy, and Eigen.
I tend to prefer the numpy model with its bias toward elegance in
higher dimensions vs convenient notation for matrix multiplication.
Explicit is better than implicit.
-- Mark Borgerding
P.S. I do not belong to any of the working groups below. If you
do and think this worth sharing, please do so.
P.P.S. Anyone learning numpy after Matlab should repeat the
mantra "numpy arrays are not matrices" a thousand times a day
until it sticks.
[1] There is the notable exception of the unsupported tensor
module. I've only used it indirectly from TensorFlow, so I cannot
comment further.
[2]
https://github.com/mborgerding/onsager_deep_learning/blob/c94ea234c7e1ac0fe9b7b3f9bc6847f406d20c4e/tools/raputil.py
[3] FWIW, the numpy matrix class overloads `operator*`, but few
people use the matrix class, instead opting for arrays and
`matmul`, `dot`, or `einsum`
[4] https://docs.scipy.org/doc/scipy/reference/linalg.html
On 2/15/19 5:51 PM, Matthieu Brucher
wrote:
Haha, it's not that we can't stand the matrix
multiplication, it's that objectively a huge chunk of the
scientists (probably the majority) is not using '*' as the
matrix multiplication (from Fortran to Numpy...).
One of the other things to notice from data scientists is
that we don't consider that everything is 2D, or 2D+ like
Matlab, but we use 1D arrays as much as nD, which makes a big
difference.
I use Eigen a lot for my personal projects that are doing
some linear algebra, and even there, it can be annoying to
jump from arrays to matrices to arrays. That's one of the
biggest advantages of lumpy, just one type (the Numpy matrix
type is not relevant nowadays).
Cheers,
Matthieu
I recently found out that the C++
standardisation committee now created a
Special Interest Group (SIG) for Linear
Algebra within SG14 (the study group for
game-dev & low-latency), and that
there's also a new study group on machine
learning, SG19.
Both groups have been formed within the
last few months, and I don't think there
was too much noise around it, so I thought
it might be quite interesting for some
people on the Eigen list. I also just
joined their forums/list,
Thanks a lot for these informations! I joined both
SG14 and SG19 lists too.
and I didn't recognise any familiar
name for the Eigen list so far.
There
is Matthieu Brucher who is a member of this list and
posted here a few times.
On a first glance, I saw that they seem
to make a few design decisions that are
different from Eigen (e.g. operator* is
only for scalar multiplications; or there
are separate row/col_vector classes
currently).
Regarding operator*, from their discussions we can
already see clear disagreements between "linear algebra"
people and more general "data scientists"... Some cannot
stand operator* as being a matrix-matrix product and are
basically seeking for a Numpy on steroid. Personally, as
I mostly do linear algebra I almost never use the
component-wise product and I'd have a hard time giving
up operator* for matrix-matrix products. On the other
hand, I found myself using .array() frequently for
scalar addition, abs, min/max and comparisons... and
I've to admit that our .array()/.matrix() approach is
not ideal in this respect.
Nevertheless, following the idea of a "numpy on
steroid", if that's really what developers want, we
might thought about making our "array world" world more
friendly with the linear-algebra world by:
- adding a prod(,) (or dot(,)) function
- moving more MatrixBase functions to DenseBase (most
of them could except inverse())
- allowing Array<> as input of decompositions
- enabling "safe" binary operations between Array and
Matrix (and returning an Array)
This way people who don't want operator as a matrix
product, or with a strong experience of numpy, could
simply forget about Matrix<>/.matrix()/.array()
and exclusively use Array<>. Then time will tell
us if, as with numpy::matrix vs numpy::array, everybody
will give up about Matrix<>... (I strongly doubt).
Gaël.
I think it would be really great to get
some people from Eigen (us!) aboard that
process.
Here are some links:
There are two repos where they started
mapping out ideas/code/paper:
Best wishes,
Patrik
--
Dr. Patrik Huber
Founder & CEO 4dface Ltd
3D face models & personal 3D face
avatars for professional applications
United Kingdom
--