How best to fix similar for triangular matrices?

[artkuo]

In response to #258, the bug was traced to mean which relies on similar (through reducedim) to prepare an array for its output, e.g. mean of a 4x3 matrix along dim 1 should yield a 1x3. It ends up calling something like similar(triangularmatrix, LowerTriangular, (1,3)) (the template matrix, the desired type, and the dimensions), which fails because a 1x3 triangular matrix doesn't make sense. In this context, similar means "give me a matrix with the same eltype" not "give me another triangular just like this one."

1. An easy fix (@tkelman) is to adjust similar so that when optional dimensions are used, it automatically drops the triangular, which results in AbstractMatrix. When no dimension is given, it could yield a lower triangular. This actually works, and even preserves sparsity automatically. I find this a bit yucky, because the return type depends implicitly on how the call is made.
2. Another fix is to change similar to always drop the triangular, again preserving sparsity where applicable. I've tried this fix and it doesn't break any existing tests (other than needing a new test for similar). Downside is it is a potentially breaking change for user code where similar is intended to mean "give me another of same structure." Yucky as well.
3. Yet another possibility is to delve into the details of how similar is called and root out the problem at the level of reducedim. Where appropriate, it could be replaced by sameeltype, meaning a more explicit version of what is meant already. In other cases, the replacement might be samecontainer or samestructure or some such. Idea is to make explicit what you really want. This seems like a longer term solution, because similar occurs 600+ times in julia project (between method calls and definitions) and it may not be straightforward to know what is intended when.

Previous discussions have complained about similar JuliaLang/julia#11574 (comment) and JuliaLang/julia#10064. The problem is that similar is vague. Present docs say that behavior for special types should be decided case by case. Right now similar for triangulars is undocumented, so the user must guess or look at the source code. Might be better to clarify and perhaps split off a sameeltype or samestructure in the long run.

There are safe, conservative ways to do option 2, such as when triangular is input, to return a matrix that's not triangular but whose off-triangle has been zeroed. This makes it safe for the user expecting a triangular back and for the most part won't blow up on them. (Triangulars store a full matrix in their .data field, so there is potential danger if things are not zeroed.)

What is the right solution?

[kshyatt]

No one responded to this issue - @andreasnoack, @jiahao, @tkelman I'm pinging the LinAlg crew. Is this still a problem? What do we want to do?

[tkelman]

I prefer options 1 or 3 here.

[kshyatt]

BTW I would also like to say that this is a very nice issue with a well stated definition of the problem & various solutions. Nice work, @artkuo.

[andreasnoack]

I would like to postpone the introduction of a new method until we have better evidence for the need so I'll vote for 1 or 2. To make progress, I'll propose that we try out 1 and see how it goes. @artkuo Do you want to try prepare a pr?

[artkuo]

Thanks for the replies, but the more I delve into it, the less I like option 1. It is a temporary bandaid, and doesn't address a fairly deep problem, which is that sometimes similar is used to allocate memory, say to hold the results of a BLAS call, and for some structured matrices the similar is not fully mutable. Meanwhile, few people seem be reporting bugs now, perhaps because the structured matrices are only used in particularly safe situations or infrequently. So I'd rather work on my solution, which is closer to option 3. I am slowly working towards a PR, which will probably be quite unpopular but at least give an idea of how wide-ranging the issue might be.

So for now I'd like to put this on hold. I'd like to keep the issue open just as a benchmark for when the real problem is fixed, hopefully in a more meaningful way than option 1.

[andreasnoack]

I think the problem is quite limited and therefore that it's best to refrain from adding complexity by introducing a new generic function. Could you provide more examples for when similar causes problems right now? I don't think anything in base relies on similar(XTriangular) being an XTriangular. Do you know of such examples?

[tkelman]

I agree with the opinion that similar is underspecified and has conflicting uses, and choosing the type it should return is all too often a subjective judgement call or dependent on how many arguments it's called with. For example the case of what to return for similar(::SharedArray) has gone back and forth, ref JuliaLang/julia#5893 and JuliaLang/julia#12964, and I think either choice has left multiple people unsettled and disagreeing. That's a pretty good sign that something is off with the abstraction.

[nalimilan]

I agree with @tkelman. Even if in practice it more or less works, we'd better clarify what similar means. (In most cases, I suspect it's used to create an Array of the same eltype and (sometimes) size. Everything would be clearer if a separate function was used for that.)

[artkuo]

@andreasnoack As an example where a mutable similar would be great, take A*B for both A and B AbstractMatrix. Right now the structured matrices try to code for every single combination, and there are quite a few missing ones. Ideally there would be a fall-back function (*){T,S}(A::AbstractMatrix{T}, B::AbstractMatrix{S}). I do have one waiting in the wings. It's 4 lines of code that will safely produce the desired result if no specific method exists, but it only if there is a proper similar (not option 1).

There is no case I could find in base where you need similar(XTriangular) to be XTriangular. The only fear is that user code expects it, which @tkelman pointed out in the bug report that led to this issue. However, I should point out that the current documentation says that similar should produce a mutable array. It also says that custom types can choose what they want to do, with default to be mutable. The current XTriangular behavior is not mutable, and is not documented as so.

If you believe that nobody should be relying on undocumented behavior, then the safer thing to do is to make it mutable across the board (option 2). I also think that the most basic and necessary closed algebra for triangulars is already implemented in triangular.jl, and does so without itself using similar. But elsewhere in Base, there is code that expects similar to return mutable, and will not work correctly with some structured matrices.

My draconian self prefers option 2, which instantly facilitates the matrix multiplication fix, is consistent with documentation, and is reliable behavior that most people should expect. I have a couple dozen examples of these and other operations between AbstractMatrix and/or AbstractVector that have no fall-back, or a fall-back that breaks because of similar (or copy) not being mutable. If you like, I can file a PR for my multiplication fix, which is extremely simple other than including option 1. I wanted to start this issue first, because I felt like this is an executive decision above my pay grade. (I also don't want to file a PR unless there is >10% chance of being accepted.)

[tkelman]

I think the documentation for similar used a poor choice of words, since mutable in the sense of Julia types is not what matters here. Triangular is an immutable-type wrapper around a mutable array, but it is writable and modifiable according to the abstraction of mutable array storage - but with additional structural information. The type system exists to provide additional information about objects, and throwing away that structural information in similar strikes me as absolutely the wrong thing to do. If you always want a dense Array without taking into account any structural properties, use the Array constructor. similar should preserve type-level structural information wherever it can.

[nalimilan]

If you always want a dense Array without taking into account any structural properties, use the Array constructor. similar should preserve type-level structural information wherever it can.

We could add a method like Array(a::AbstractArray, T=eltype(a), dims=size(a)) = Array{T}(dims), deprecate similar in favor of it, and add a new stricter version of similar always returning an AbstractArray of the same type and structure.

[andreasnoack]

@artkuo I think it makes sense to always return a (completely) mutable array determined by the wrapped array, i.e. option 2, but I also wanted to close #258 and suggested option 1 as a compromise with little harm done. #258 is really a tiny issue and we'll probably never be able to completely avoid issues like that. I don't really want to defend option 1, but wouldn't your matrix multiply fix use a dimension argument such that option 1 would actually work?

@tkelman @nalimilan Array wouldn't be generic enough for e.g. XTriangular{DArray}. I'm not sure that throwing away the wrapper type is such a disaster. The place where I expected similar(XTriangular)=XTriangular to be useful was in linalg/triangular.jl and I realized that I didn't even use the method there. When do you think it would be useful? As things are now, you can't use a generic setindex! strategy because the structurally zeroed part of the matrix will give you errors.