Triangular: no implementation of similar, setindex!

[timholy]

I had code that looked like this:

L = chol(A, :L)
X = similar(L)

This broke when chol started returning a Triangular type. I started to implement similar for Triangular (creating another Triangular type), but then noticed there is no setindex!. Since the next thing my code did was stuff some values in X, that would fail.

So before doing anything, I should ask: is the omission of setindex! deliberate?

[andreasnoack]

Initially, Triangular was for dispatch when solving linear systems and therefore setindex! was not necessary. I think we should aim for full coverage of the special matrices so please continue.

By the way, I can see that chol(Matrix) in line 44 returns a Matrix. That must be a mistake. That line could be removed and the method in line 43 could have :U as default value.

[jiahao]

Isn't Triangular immutable now? Is it a good idea to define setindex! on an immutable type?

[StefanKarpinski]

Depends on what you mean by immutable. Do you mean that the type is an immutable Julia type? Or that it is immutable-by-convention like strings?

[andreasnoack]

ArrayViews are also immutable and it would be a significant limitation if setindex! couldn't work for them.

However, I see another problem here. How should we handle the situation when someone tries to set an element in the zero part of the matrix? An error, just ignore it since it doesn't change anything or return a non-Triangular type.

[StefanKarpinski]

Yes, that's the only reasonable thing to do. I think we should allow changing the parts of these special kinds of matrices that can be modified because someone who knows what they're doing may want to do that – and it's often easier to modify a diagonal or off-diagonal element using matrix coordinates than trying to reach in and change the underlying arrays directly.

[jiahao]

Triangular, like many of the other special matrices, are currently defined as an immutable Julia type.

How should we handle the situation when someone tries to set an element in the zero part of the matrix? An error, just ignore it since it doesn't change anything or return a non-Triangular type.

@JeffBezanson didn't like the idea of setindex! mutating the type, and I can see why (it can be quite an expensive type change in general, e.g. from a special matrix like Diagonal to the dense Matrix). Triangular is a little special though, since we don't use a compact memory representation for them (the full matrix is stored) and the type contains a label saying which half (upper or lower triangular) should be considered. So I think I would prefer a silent update of the ignored half for Triangular, and throwing an error in general.

[StefanKarpinski]

I'm not suggesting changing the type of an object – you can't do that. The definition of Triangular is:

immutable Triangular{T<:Number} <: AbstractMatrix{T}
    UL::Matrix{T}
    uplo::Char
    unitdiag::Char
end

The fact that Triangular is immutable means that you can't change its fields – once you construct a Triangular object, its .UL field will always refer to the same Matrix object and it's .uplo and .unitdiag will always be the same Char values. The .UL Matrix object is, however, still mutable. What I'm suggesting is that you can allow mutation of the underlying .UL data via calling setindex! on the Triangular object. If ut is an upper-triangular matrix and you do ut[i,j] = v where i <= j it would do ut.UL[i,j] = v; if, on the other hand, j < i, that would be an error.

[andreasnoack]

@StefanKarpinski I think that would be the best solution for Triangular. However, Hermitian/Symmetric is more tricky. Should setindex! be defined? When setting element i,j are you then automatically setting j,i?

@timholy How far are you in this implementation? I have a rewrite of Triangular on my todo list because I want to parametrise it on the matrix type and probably also split out upper and lower Triangular matrices into two types. If you haven't written your changes yet, I can try to rewrite Triangular today.

[timholy]

I haven't even had a chance to start yet, so do go ahead.

[jiahao]

However, Hermitian/Symmetric is more tricky. Should setindex! be defined? When setting element i,j are you then automatically setting j,i?

Setting A[i,j] and A[j,i] concurrently seems like a recipe for guaranteeing poor performance due to the alternating memory access pattern

[andreasnoack]

Setting A[i,j] and A[j,i] concurrently seems like a recipe for guaranteeing poor performance due to the alternating memory access pattern

I didn't manage to write what I had in mind. A Hermitian only references one the triangles of the storage matrix, so I am not suggesting changing elements of the storage matrix that are not referenced. The problem is that when changing an element of a Hermitian matrix the matrix is (in general) no longer Hermitian. However, it could implicitly be assumed that you change j,i element as well such that the matrix remains Hermitian, but I don't think that is a good idea.

[jiahao]

Ah, that makes more sense. I do think that setindex! on Symmetric/Hermitian, if supported, should be thought of as implicitly changing both elements, rather than mutating the matrix to be no longer of the original symmetry.

[bisraelsen]

Hopefully this is the right place for this.

Was setindex! ever implemented for Symmetric? I haven't been able to figure out how to use it. If it hasn't been implemented, is there another way?

[jiahao]

Not yet

julia> methods(setindex!, (Symmetric,))
0-element Array{Any,1}

PRs to implement setindex! welcome.

[bisraelsen]

Could you give me an idea of where to start. I would be happy to work on it, although I am fairly new to the language.

EDIT: So, I found symmetric.jl I think this is the right place. I'll start there and try to use ffc9778 as a template.

[andreasnoack]

@bisraelsen Could you sketch how you think setindex! should work for Symmetric? If you change an off-diagonal element then it is not symmetric anymore so mimicking setindex! for XTriangular would only allow you to edit the diagonal.

[bisraelsen]

I only meant to use the XTriangular solution as a guide for me to add setindex! support for Symmetric, as I am new to the project. However, setindex! definitely should not have the same implementation.

I had planned only on allowing changes that maintain symmetry. For example if someone changes index (i,j) then (j,i) would also be changed...