Applied

[dlfivefifty]

This adds Applied, which will fix #11 and #14.

[codecov-io]

Codecov Report

Merging #17 into master will increase coverage by 5.5%.
The diff coverage is 53.89%.

@@            Coverage Diff            @@
##           master      #17     +/-   ##
=========================================
+ Coverage   50.57%   56.07%   +5.5%     
=========================================
  Files          11       13      +2     
  Lines         785      815     +30     
=========================================
+ Hits          397      457     +60     
+ Misses        388      358     -30

Impacted Files	Coverage Δ
src/lazybroadcasting.jl	43.47% <ø> (-2.36%)	⬇️
src/linalg/lazymul.jl	3.22% <0%> (+3.22%)	⬆️
src/linalg/blasbroadcasting.jl	57.14% <16.66%> (+3.8%)	⬆️
src/linalg/blasmul.jl	50.32% <40%> (-0.63%)	⬇️
src/linalg/inv.jl	50% <44.23%> (ø)	⬆️
src/linalg/add.jl	55% <55%> (ø)
src/linalg/mul.jl	58.62% <68.18%> (+12.96%)	⬆️
src/lazyapplying.jl	76.31% <76.31%> (ø)
... and 3 more

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 9d0fb6a...5c630b0. Read the comment docs.

[dlfivefifty]

@tkf Could you have a look and say if you have any suggestions for the changes made to Add. I've made it an Applied which makes it consistent with Mul.

[tkf]

Is Add still the primary API? Or users should call AddArray instead of Add? I'm guessing that Add is the primary one as it looks like so for Mul.

If that's the case, isn't it better to test this with Add (or maybe just test it with both)? When I tested it locally with Add, there are some errors.

[dlfivefifty]

Yes good point. This needs a redesign of the code, though it's touching against a broader point: given a general applied "tree" of operations, how's best to simplify it? For example, in this case what we actually want is

applied(*, applied(+, A, B,C) , D)

to be "simplified" to

applied(+, applied(*, A, D), applied(*, B, D), applied(*,C,D))

Another such case is choosing multiplication orders to minimize cost.

[tkf]

I'd suggest to store the tree of operations as-is and do the simplification/transformation inside materialize. We can then define "optimization passes" to look at the tree and turn it into more efficient one (maybe in the instantiate phase?). If there is no transformation at "tree build time" then it's straight-forward to implement your own case-specific materialize function for ApplyArray or Applied.

Or, maybe you are wondering how ApplyArray and Applied interact? I think it's better to always use ApplyArray in the outer most tree. That is to say,

ApplyArray(applied(*, ApplyArray(applied(+, A, B, C)) , D))

is "canonicalized" to

ApplyArray(applied(*, applied(+, A, B, C) , D))

at the tree build time. And, repeating my point in the first paragraph, it then is turned into

applied(+, applied(*, A, D), applied(*, B, D), applied(*, C, D))

inside materialize later.

[dlfivefifty]

At this point in time I think we should only support 1 version, which should be applied.

Though your fast implementation of (A+B+C)*D is arguably premature to include, without having a well-defined "optimization" infrastructure. Is this feature crucial to you, or could that part be moved to your own code for the time being?

[tkf]

I can do, e.g.,

eq = Mul(Add(A, B, C), D)
my_custom_materialize!(Y, eq)

so it's not super crucial. If my code is getting in a way to re-structure LazyArrays, I'm fine with removing the code until there is a better way to do it.

[dlfivefifty]

It's not that it's getting in the way, so I can update it. It's more that it's an outlier at the moment. But maybe it's better to update it so your code isn't made more complicated, and we can delete it later if we decide it was a mistake.

[dlfivefifty]

PS Arguably the lowering of c .= Mul(A,b) to materialize!(MulAdd(1.0, A, b, 0.0, c)) also fits the pattern of a "simplifying" operation. Probably MulAdd should be type-aliased to applied(+, applied(*, α, A, b), applied(*, β, c))

[tkf]

it's an outlier at the moment

I was kind of hoping that LazyArrays.jl would accumulate optimized/fused operations. For example, MatrixChainMultiply.jl or similar can be hooked into materialize! mechanism. I also wrote a few fused operations for sparse matrices which can in principle be used on the computation graph generated by LazyArrays. But adding such optimizations in an organized manner would be a hard project by itself. So I thought solidifying "operation tree" data structure and API would be much more important so that code outside LazyArrays can use it.

Probably MulAdd should be type-aliased to applied(+, applied(*, α, A, b), applied(*, β, c))

I think it makes sense to treat MulAdd to be the "lowered representation" of c = αAb + βc when c in the lhs and rhs are detected to be the same object in materialize!. That is to say, tree of operations representing the rhs αAb + βc always is applied(+, applied(*, α, A, b), applied(*, β, c)) and it's lowered to MullAdd or copy!+MullAdd depending on if the destination array is c or not. But this part is more like an implementation detail so I'm sure you know better.

[dlfivefifty]

I was kind of hoping that LazyArrays.jl would accumulate optimized/fused operations

Yes, I think that's a great idea. But not in this PR 🤣

For example, MatrixChainMultiply.jl or similar can be hooked into materialize! mechanism.

The ApplyStyle setup gives a lot of flexibility for allowing different behaviour, so it's possible to do this without it (yet) being supported in LazyArrays.jl.

So I thought solidifying "operation tree" data structure and API would be much more important so that code outside LazyArrays can use it.

I think we're a long ways from "solidifying": we'll probably have to iterate the design based on usage patterns.

think it makes sense to treat MulAdd to be the "lowered representation"

OK I think I'll leave it as is for now.

[tkf]

we'll probably have to iterate the design based on usage patterns.

Ah, yes, I wasn't trying to argue LazyArrays.jl should become 1.0 right after this PR.

[zsoerenm]

Great work!
Are there plans to optimize Mul(A, A') to BLAS.syrk('U', 'N', 1, A)?

I often have weighting matrix W too, which is either diagonal or symmetric. Unfortunately there is no BLAS equivalent for Mul(A, W, A'). But could I force the output to be symmetric something like B = Symmetric(Matrix{Float64}(undef, 10, 10)); B .= Mul(A, W, A')?

[dlfivefifty]

Yes, supporting custom matrix types is a big motivation (especially when they live on the GPU or parallel memory). The syntax you asked for should eventually work, though also possible is the shorthand SymmetricMatrix(Mul(A,W,A’)).

An example of this already in action is converting MulArray of banded matrices to a BandedMatrix.