Add tril and triu for special matrices

[kshyatt]

As discussed in JuliaLang/LinearAlgebra.jl#161, Diagonal, Bidiagonal, Tridiagonal, SymTridiagonal, UpperTriangular, LowerTriangular, UnitUpperTriangular, UnitLowerTriangular, Symmetric, and Hermitian types didn't have tril or triu defined. Diagonal had tril! and triu! but the ! was incorrect since the underlying matrix wasn't modified.

I also added tests.

[jakebolewski]

👍 My only complaint would be that it might be too readable.

[andreasnoack]

👍

[StefanKarpinski]

My only complaint would be that it might be too readable.

heh

[ScottPJones]

Maybe this could be:

    elseif k > 0
          return M
    elseif M.isupper
          if k < 0
               return Bidiagonal(zeros(M.dv),zeros(M.ev),true)
          else
               return Bidiagonal(M.dv,zeros(M.ev),true)
          end
    elseif k < -1
       return Bidiagonal(zeros(M.dv),zeros(M.ev),false)
    elseif k == -1
       return Bidiagonal(zeros(M.dv),M.ev,false)
    else
       return M
    end

This avoids the rechecking of M.isupper and passing it to Bidiagonal.

[tkelman]

Should these generally be returning copies of or references to the input data?

[tkelman]

Answering my own question, judging by

julia/base/linalg/generic.jl

Lines 41 to 44 in db787d8

	triu(M::AbstractMatrix) = triu!(copy(M))
	tril(M::AbstractMatrix) = tril!(copy(M))
	triu!(M::AbstractMatrix) = triu!(M,0)
	tril!(M::AbstractMatrix) = tril!(M,0)

I think this really needs to be consistent about copying rather than returning a reference. Sometimes-copy, sometimes-reference is not good.

edit: Or these can consistently implement tril!(A, k) for various subtypes of AbstractMatrix and let the general version take care of the copying for tril.

[kshyatt]

@tkelman so I should have all of these return copy(M) instead of M? I'm happy to submit a new PR to do that.

[tkelman]

I think the right solution here is to define tril!(A, k) for the various specialized subtypes, which is free to do in-place modification, return references, etc. The general-purpose fallback definitions for tril!(A::AbstractMatrix) and tril(A::AbstractMatrix, k) tril(A::AbstractMatrix) should handle the copying for all of the subtypes.

edit: apparently there isn't a fallback for tril(A::AbstractMatrix, k) right now, but there probably could be tril(A::AbstractMatrix, k) = tril!(copy(A), k) ?

[kshyatt]

I guess the reason I kept many of these as tril and not tril! is for cases like:

elseif !M.isupper && k == 0
        return M
elseif M.isupper && k == 0
        return Bidiagonal(M.dv,zeros(M.ev),M.isupper)

where we either return M or a new matrix. It seems like it might make sense to just copy M in the return? Or should we just do M = Bidiagonal(M.dv,zeros(M.ev),M.isupper); return M in the lower one and put the ! in?

[tkelman]

I think if one of these methods ever returns the original matrix itself or a reference to some part of it in the output, then that method should be named tril! (resp. triu!) to denote the resulting potential for aliasing. If it only ever returns copies then it can be a method of tril (resp. triu).

There are cases where tril or triu could get away without copying the entire input matrix, and only copy some parts of it like in the case you've shown. Do those deserve extra versions of the same basic method implementations? Maybe it could be simple and short with macro loops where tril! is defined in terms of running identity and tril is defined in terms of running copy on any piece that references the original matrix.

[StefanKarpinski]

It looks like mutating versions of all of these could modify the argument matrix in-place and always return the argument matrix. In which case, the non-mutating version could always just be written as calling the mutating version on a copy of the input matrix. Might be slightly slower since the input has to be copied and you can't just allocate vectors of zeros, but I'm not sure that it would be a major difference.

[tkelman]

I think I like that option. Instead of using zeros, the in-place versions could do fill!(0,.

The only tricky case is the non-type-stable branches for SymTridiagonal.

[StefanKarpinski]

Perhaps tril and triu for SymTridiagonal should always just return Tridiagonal (which implies copying). Making it type-stable seems likely to be better for performance than the more precise type.

[andreasnoack]

I thank I didn't read carefully enough yesterday. There are some problems here, e.g.

julia> A = LinAlg.UnitLowerTriangular(randn(4,4))
4x4 Base.LinAlg.UnitLowerTriangular{Float64,Array{Float64,2}}:
 1.0        0.0        0.0       0.0
 1.1291     1.0        0.0       0.0
 0.195707  -1.72338    1.0       0.0
 3.31616    0.252167  -0.598446  1.0

julia> triu(A, 1)
4x4 Base.LinAlg.UnitLowerTriangular{Float64,Array{Float64,2}}:
 1.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0
 0.0  0.0  1.0  0.0
 0.0  0.0  0.0  1.0

julia> triu(full(A), 1)
4x4 Array{Float64,2}:
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0

and I'm in doubt if we can, in general, support the trix(AbstractMatrix, Integer) methods for the special matrix types. A solution could be to split up the methods such that triu(A,0) returns a Matrix and triu(A)) potentially returns a special matrix type. That should be type safe.

@tkelman I think you might be extending the meaning of !. Usually the meaning is that the input is modified, not that a potential aliasing is happening. I can try to find some earlier discussions of this, but I think the view is that your solution would cause too many copies and that you should explicitly copy when necessary.

[tkelman]

That should be type safe.

I think you mean type stable.

I think you might be extending the meaning of !. Usually the meaning is that the input is modified, not that a potential aliasing is happening.

True. So in cases where in-place modification is possible for these methods then tril! should do so, but I think tril should always copy. We really shouldn't be leaving aliasing up to the user to figure out and doing a sometimes-copy, sometimes-reference. Shouldn't have to read the source to determine whether or not modifying the output of tril is going to change the input.

[kshyatt]

I'm working on a PR to address your feedback. For the Unit{Lower/Upper}Triangular types we also cannot maintain type stability - should I have trix! return a {Upper/Lower}Triangular as appropriate?

[tkelman]

I think so, that sounds appropriate.

[kshyatt]

Ok, almost done!

[tkelman]

I don't think these are right. Aren't half the elements in A.data potentially garbage?

julia> A = rand(4,4)
4x4 Array{Float64,2}:
 0.992492  0.88037   0.40931   0.311544
 0.623913  0.93559   0.310044  0.598238
 0.393736  0.051534  0.103072  0.74386
 0.747731  0.684462  0.813315  0.0329609

julia> Symmetric(A)
tril(4x4 Base.LinAlg.Symmetric{Float64,Array{Float64,2}}:
 0.992492  0.88037   0.40931   0.311544
 0.88037   0.93559   0.310044  0.598238
 0.40931   0.310044  0.103072  0.74386
 0.311544  0.598238  0.74386   0.0329609

julia> tril(Symmetric(A))
4x4 Array{Float64,2}:
 0.992492  0.0       0.0       0.0
 0.623913  0.93559   0.0       0.0
 0.393736  0.051534  0.103072  0.0
 0.747731  0.684462  0.813315  0.0329609

julia> full(Symmetric(A))
4x4 Array{Float64,2}:
 0.992492  0.88037   0.40931   0.311544
 0.88037   0.93559   0.310044  0.598238
 0.40931   0.310044  0.103072  0.74386
 0.311544  0.598238  0.74386   0.0329609

julia> tril(full(Symmetric(A)))
4x4 Array{Float64,2}:
 0.992492  0.0       0.0       0.0
 0.88037   0.93559   0.0       0.0
 0.40931   0.310044  0.103072  0.0
 0.311544  0.598238  0.74386   0.0329609

[kshyatt]

Wow. K, on it.

The lower triangle is garbage.