Get rid of ComponentMatrix and higher-order ComponentArrays?

[gdalle]

This is a very disruptive proposition, I know, but hear me out. The reason I'm suggesting it is because higher-order ComponentArrays are

1. still a bit buggy
2. not widely used
3. keeping Base from doing its job

Still a bit buggy

For as long as I can remember, the doc has stated:

Higher dimensional ComponentArrays can be created too, but it's a little messy at the moment.

Indeed, ComponentMatrix gives rise to a number of issues like

Indexing

• FR: Preserve CA-ness when indexing component matrix with shaped components #257
• Indexing ComponentMatrix with FlatAxis components #248
• Mat index bug #249

Concatenation

• ComponentVectors are not stack-able #254
• hcat: type unstable based on the number of arguments #168
• Need reduce(hcat, ... method #113
• Fix *cat inconsistencies #169

Linear algebra operations

• Reexport inv and diag etc from LinearAlgebra? #187 (comment)
• Component Arrays with DAE's #128 (comment)
• [wip] Fix LinearAlgebra.lmul! on GPU #233
• Change Adjoints to be ComponentArrays #170

Basically, as soon as we touch matrices, the number of necessary overloads grows out of control, and no one has time to play whack-a-mole.

Not widely used

If we take a look at the search

https://github.com/search?q=language%3AJulia+ComponentMatrix&type=code

and exclude forks, we find exactly 4 repos where ComponentMatrix is used:

• MuseInference.jl by @marius311: one linalg overload
• MacroModeling.jl by @thorek1: one commented line
• Flight.jl by @e271828e: seems like a legitimate use
• MTKHelpers.jl by @bgctw: not sure what this does exactly but it has type piracy

In addition, @vpuri3 and @dingraha have recently submitted PRs related to ComponentMatrix, so perhaps they use it too.

Keeping Base from doing its job

In several of the bugs outlined above, @mcabbott and I think that just removing the higher-order ComponentArrays would allow functions from Base (like stack or hcat) to take over and do the right thing. But that is still an untested hypothesis.

---

Upsides of ComponentMatrix

On the other hand it is true that this format has good things going for it. Most notably, the ability to index block by block, as demonstrated in the ODE example

[e271828e]

For what it's worth, ComponentMatrix is extremely useful to me. When working with state-space control methods, it allows me to easily extract decoupled reduced models from a complete linear model. Afterwards, it makes the design process itself much less tedious and error-prone.

I understand that higher-dimensional ComponentArrays might not generalize as gracefully as ComponentVectors, but removing them would disrupt non-trivial parts of my workflow.

By the way, @jonniedie, huge thanks for this package. I absolutely love it. If I hadn't found it, I would have tried to develop something similar myself (likely with much worse results).

[vpuri3]

I'm not using component matrix at the moment, although I do like the ability to do block index matrixes.

is it possible to classify ComponentMatrix as an experimental feature with limited support?

[jonniedie]

It's an interesting proposal and one I've considered regularly due to how many problems it's caused over the years. I don't think I could remove it at this point, though. Especially because I personally use it on my own work (like @e271828e, I use it controls work as the state/control jacobians from the linearization of our ComponentArray-based nonlinear sim can be accessed by their blockwise components). But maybe designating it as an experimental feature is a good idea at this point as I've haven't had much time to respond to issues lately.

[gdalle]

Thanks for your answers! As I suspected, there are few people using this functionality but it is important to them. Marking it as experimental would probably be wise.

The main reason I am asking is because I'm developing DifferentiationInterface.jl, which has bindings for every single autodiff package... as long as you use arrays. For anything more complicated, I recommend putting everything inside a ComponentVector. The trouble is that I sometimes need to stack several of those to create Jacobians or Hessians, and that's why I opened #254. Basically, even when you don't work with a ComponentMatrix, it can creep up on you unexpectedly, or cause errors when you combine ComponentVectors. If you give me some pointers, I can try to cook up a PR myself

[marius311]

Thanks for the effort to track down the MuseInference usage and tagingg me so I see this. Fwiw if you could see private repos youd find a bunch more usage and that one overload probably undersells how much ComponentMatrix gets used even in that package, lots of the arrays flowing through there are ComponentMatrices. Mainly in the context of covariance matrices with parameters in each dimension labeled. So Id love it this could be kept around, I use it alot, even if ocassionaly theres papercuts.

(Echoing that this package is amazing, easily one of my top 5 Julia usability packages)

[gdalle]

Yeah I forgot to emphasize it but I obviously am a huge fan of ComponentArrays too!

[vpuri3]

@gdalle if there is somebody who wants to make a conceited effort to fix/reimagine ComponentMatrix, that could resolve this problem

[gdalle]

Of course but I don't think it will be me unfortunately. The reason I suggested a clean slate was because I thought Base might handle the matrix manipulations better than the not-quite-complete ComponentMatrix implementation, but if it is useful to people we're definitely keeping it

[dingraha]

Thanks for pinging me on this discussion. Yes, I'm a big fan of ComponentMatrix and ComponentArrays.jl in general. I would be sad if we lost ComponentMatrix, but it does everything I want at the moment, so I wouldn't mind if it was marked experimental or split off into a separate package.

@gdalle I was able to fix (I think) the issue with stack and ComponentVectors with a few small tweaks to this PR: #249. I've only tried stacking ComponentVectors, though: https://github.com/dingraha/ComponentArrays.jl/blob/ab68e9a0ae9190524127206e1b7ab3bd4c6c9f1a/test/runtests.jl#L694

[dingraha]

A counter-proposal: instead of getting rid of ComponentMatrix and higher-dimensional ComponentArrays, what about restricting the functionality of ComponentArray? The idea would be to treat

• indexing with symbols
• getproperty with symbols
• in-place operations

"specially" as they are now, but have everything else return plain arrays. So, for example, we would expect stack, similar, all the various *cat functions, broadcasting, zero, *, \, / would just return plain arrays (or whatever is appropriate for the object that getdata returns).

The motivation: my use case for ComponentArrays is to keep track of inputs and outputs of functions that will be passed to ForwardDiff.jl or other AD libraries, which require one array input and one array output. The functions I pass will, at the beginning, unpack the inputs based on ComponentArray names and then write the outputs into the output ComponentArray at the end. In between, all the ComponentArray functionality isn't necessary. For example:

function f!(Y, X)
    # unpack inputs
    a = @view X[:a]
    b = @view X[:b]

     # Do something with a and b

    # write outputs to Y
    @view(Y[:c]) .= c
    @view(Y[:d]) .=  d

    return nothing
end

One of the downsides of this proposal is that we couldn't use the adjoint trick to create ComponentMatrix things as described in the docs https://jonniedie.github.io/ComponentArrays.jl/dev/quickstart/. But that wouldn't be too much harder:

julia> x
ComponentVector{Int64}(a = 1, b = [2, 3], c = [4, 5])

julia> y
ComponentVector{Int64}(d = 6, e = [7, 8, 9])

julia> J = ComponentMatrix(getdata(y) .* getdata(x)', getaxes(y)..., getaxes(x)...)
4×5 ComponentMatrix{Int64} with axes Axis(d = 1, e = ViewAxis(2:4, Shaped1DAxis((3,)))) × Axis(a = 1, b = ViewAxis(2:3, Shaped1DAxis((2
,))), c = ViewAxis(4:5, Shaped1DAxis((2,))))
 6  12  18  24  30
 7  14  21  28  35
 8  16  24  32  40
 9  18  27  36  45

julia> 

Obviously this would be a very breaking change. But I personally would much prefer it to removing ComponentMatrix etc..

[dingraha]

Ah, another downside would be that code like the ODE example in the docs https://jonniedie.github.io/ComponentArrays.jl/dev/examples/ODE_jac/ wouldn't work anymore. I'm not sure how the ComponentMatrix Jacobian is getting created there (must be using the x .* x' trick internally in DifferentialEquations?). I wonder if libraries like DifferentialEquations allow the user to communicate how the Jacobian should be constructed?

[thorek1]

I used it to get AD working but once I needed sparse matrices I abandoned it again and wrote a bespoke version instead.
here is an example to pack multiple (sparse) matrices into once vector and pass them to a function (then unpacked inside the function again)

r1,c1,v1 = findnz(ŝ_to_ŝ₃)

coordinates = Tuple{Vector{Int}, Vector{Int}}[]
push!(coordinates,(r1,c1))
        
dimensions = Tuple{Int, Int}[]
push!(dimensions,size(ŝ_to_ŝ₃))
push!(dimensions,size(C))
        
values = vcat(v1, vec(collect(-C)))