Re: The "star" problem, was Re: [eigen] ISO C++ working groups on Linear

Le sam. 16 févr. 2019 à 13:41, Mark Borgerding <mark@xxxxxxxxxxxxxx> a écrit :

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

Le jeu. 14 févr. 2019 à 18:29, Gael Guennebaud <gael.guennebaud@xxxxxxxxx> a écrit :

Hi Patrik,

On Wed, Feb 13, 2019 at 6:19 PM Patrik Huber <patrikhuber@xxxxxxxxx> wrote:

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:

SG14 mailing list: http://lists.isocpp.org/sg14/

SG19 mailing list: https://groups.google.com/a/isocpp.org/forum/?fromgroups#!forum/sg19

There are two repos where they started mapping out ideas/code/paper:

https://github.com/BobSteagall/wg21/tree/master/linear_algebra

https://github.com/cpp-linear-algebra

Best wishes,

Patrik

--

Dr. Patrik Huber

Founder & CEO 4dface Ltd

3D face models & personal 3D face avatars for professional applications

United Kingdom

Web: www.4dface.io

--

Quantitative analyst, Ph.D.
Blog: http://blog.audio-tk.com/
LinkedIn: http://www.linkedin.com/in/matthieubrucher