*(a, b, c, xs...) should have methods for matrices and vectors

[andreasnoack]

There are two problems right now. First of all, it would be great to have A*B*x=A*(B*x) instead of (A*B)*x because the former is Θ(n^2)-ish whereas the latter is Θ(n^3)-ish. The second problem is that the optimized transposed multiplication methods not called for e.g. A'*B*C. Right now, A is transposed before the multiplication so typically (A'*B)*C is much faster than A'*B*C.

[StefanKarpinski]

As a generalization of this, we could use the classic dynamic programming
algorithm to decide what order to multiply matrices in. This could be
memoized, or it might be fast enough to compute every time. I'd guess that
it is since the sizes of arrays needed to compute this would generally
dwarf those required for the new matrix.

On Wednesday, July 8, 2015, Andreas Noack notifications@github.com wrote:

There are two problems right now. First of all, it would be great to have
A_B_x=A_(B_x) instead of (A_B)_x because the former is Θ(n^2)-ish whereas
the latter is Θ(n^3)-ish. The second problem is that the optimized
transposed multiplication methods not called for e.g. A'_B_C. Right now, A
is transposed before the multiplication so typically (A'_B)_C is much
faster than A'_B_C.

—
Reply to this email directly or view it on GitHub
#227.

[tpapp]

Matrix multiplications with floats are not associative, so technically A*(B*x) is not equal to (A*B)*x. My understanding is that the proposal would make Julia select the grouping with the best performance, but this would make A op B op C either (A op B) op C or A or (B op C) depending on the values of A,B, and C for op=*. I don't know if there is a name for what this violates... not type stability, but something similar. That said, I can only think of very contrived examples for when this matters in practice.

[jiahao]

Everyone knows that floating point multiplication is not associative, but the reason why it rarely matters in practice is because the error in nonassociativity is bounded essentially by two ulps. See JuliaLang/julia#12375 (comment) for a reference to TAoCP.

[LMescheder]

I'm new to Julia and took this as an exercise to get used to it. I created a simple package that uses the dynamic programming algorithm mentioned by Stefan. It think it works quite well and gives huge speedups, when the default ordering is suboptimal. Maybe you find it useful!

[kshyatt]

Hi @LMescheder, this looks cool. I don't see your package in https://github.com/JuliaLang/METADATA.jl - have you registered it yet?

[StefanKarpinski]

Yes, very cool! I would note that the left-to-right multiplication result is not necessarily more correct, in fact, I would guess that the method that uses fewer operations would tend to be more accurate, but I'm sure that people like @alanedelman, @stevengj, @andreasnoack and @jiahao could say better than I.

[jiahao]

Well, oestensibly the reason for our special parsing of *(varargs...) was precisely to enable tricks of the form @LMescheder has played with, so I'm glad that someone has been able to make use of it! In theory, reordering is unsafe because of the nonassociativity of floating point roundoff, but in practice it doesn't seem to be a big problem in the end. And @alanedelman has yet to see a case where reordering is problematic in real problems.

[StefanKarpinski]

You can always force a particular association with parentheses.

[jiahao]

Yes, that seems to be what is done in practice, rather than being left ambiguous for the reader to deduce from parsing rules.

[LMescheder]

What I really like about the idea is, that it is very reasonable to interpret the absence of parentheses as "I don't care about association". In fact, I think this interpretation is far more natural than the usual "act as though I wanted multiplication from the left" interpretation (even though most programming languages and frameworks interpret it in the second way).
@kshyatt I've just created a pull request to METADATA.jl

[LMescheder]

However, it should be noted, that the dynamic programming algorithm's running time is cubic and the required space quadratic in the number of factors. So there should be some heuristic preventing the algorithm from running if someone e.g. wants to multiply 10000 2x2 matrices...

[Jutho]

I guess nobody will explicitly write the multiplication of 10000 matrices as
A1 * A2 * ... * A10000, i.e. the heuristic might be that the algorithm is worth running for as many multiplications as people can actually write after each other. Unless of course some kind of code generation is used, but that certainly seems bad practice.

[andreasnoack]

Excluding sqaure matrices seems like good conservative heuristic. Have you ever looked at the Hu and Shing (1981) algorithm under the More Efficient Algorithms section on Wikipedia? Apparantly it's n log n instead of n^3.