Vector of upper triangle

General Usage
1

Is there an elegant way of extracting the upper triangle of a matrix as a vector using Base? Besides, of course, writing something like

function vec_triu(M::AbstractMatrix{T}) where T
    m, n = size(M)
    m == n || throw(error("not square"))
    l = n*(n+1) ÷ 2
    v = Vector{T}(l)
    k = 0
    for i in 1:n
        v[k + (1:i)] = M[1:i, i]
        k += i
    end
    v
end

so that

julia> vec_triu(reshape(1:9, 3, :))
6-element Array{Int64,1}:
 1
 4
 5
 7
 8
 9
2 Likes
2

M[triu(trues(M))]

See also Half Vectorization

3 Likes
3

Not that I know, but just FYI, you need to write out the loop; your function allocates (on 0.7 from this week). This does not even make your code longer or harder to read.

@inbounds for ell=1:i v[k + ell] = M[ell, i] end
1 Like
4

[A[i, j] for j in indices(A, 1) for i in 1:j]

It’s allocating a little more than expected, though.

2 Likes
5

Thanks everyone for the replies. Results below, related code in gist.

function time (ns) allocation
vec_triu 232 528
vec_triu_collect 278 848
vec_triu_inbounds 201 528
vec_triu_loop 37 128
2 Likes
6

Hi everyone, I needed this functionality a lot recently and extracted my code for this into a separate package (I hope this is not overkill) but I removed the allocations in my implementation.
I hope it helps someone else as well and feedback is much appreciated.

7

It seems risky to use @inbounds together with one-based indexing, while accepting AbstracArrays.

I think it would be wise to either restrict inputs to one-based, or use general indexing, or remove @inbounds/@turbo.

Simplest solution may be to sprinkle some require_one_based_indexing calls in there.

8

Thank you very much for your feedback. I have added calls to require_one_based_indexing in (hopefully) appropriate places.

1 Like