Iterating over elements of upper triangular matrix, but cartesian indices are needed

lmiq · July 29, 2021, 3:00pm

Basically what I (think) I need is this function:

https://stackoverflow.com/questions/27086195/linear-index-upper-triangular-matrix

Of course I can just implement it, two lines of code. But that seems something that probably is needed quite frequently in matrix computations, so I am surprised I didn’t find it as a Base function.

First I thought that most people would prefer constructing an iterator like (for a 3x3 matrix, for example):

julia> for c in Iterators.filter(c -> c[1] < c[2], CartesianIndices((1:3,1:3)))
         @show c
       end
c = CartesianIndex(1, 2)
c = CartesianIndex(1, 3)
c = CartesianIndex(2, 3)

That would be fine for my current purposes, except that it does not play well with @threads:

julia> Threads.@threads for c in Iterators.filter(c -> c[1] < c[2], CartesianIndices((1:3,1:3)))
         @show c
       end
ERROR: TaskFailedException
Stacktrace:
 [1] wait
   @ ./task.jl:322 [inlined]
 [2] threading_run(func::Function)
   @ Base.Threads ./threadingconstructs.jl:34
 [3] top-level scope
   @ ./threadingconstructs.jl:93

    nested task error: MethodError: no method matching length(::Base.Iterators.Filter{var"#19#20", CartesianIndices{2, Tuple{UnitRange{Int64}, UnitRange{Int64}}}})

What I am missing here? Should I just use a custom implementation of the indexing, is there a Base function already for that, or some alternative way to iterate that is accepted by @threads?

mcabbott · July 29, 2021, 4:37pm

I don’t think there’s a built-in iterator for this, although it’s been discussed (including for sparse arrays) e.g. here but probably elsewhere too. Maybe someone knows a better link.

While it won’t directly help with threading, something like (CartesianIndex(x,y) for x in 1:3 for y in x+1:3) (i.e Iterators.flatten) may be more efficient than filtering. You can use Iterators.partition and Threads.@spawn to multi-thread this.

lmiq · July 29, 2021, 6:26pm

Thanks. The custom function resulted to be the most practical alternative for now. For the records, it is this one:

@inline function iuppert(k::Integer,n::Integer)
  i = n - 1 - floor(Int,sqrt(-8*k + 4*n*(n-1) + 1)/2 - 0.5)
  j = k + i + ( (n-i+1)*(n-i) - n*(n-1) )÷2
  return i, j
end

julia> iuppert(1,3)
(1,2)

julia> iuppert(3,3)
(2,2)

julia> iuppert(3,4)
(1,4)

Per · June 23, 2022, 2:38pm

I also needed an iterator over indices of an upper triangular matrix, so after finding this thread, I wrote this small package: GitHub - perrutquist/TriangularIndices.jl

It uses a slightly different formula (compared to the above) which generates the indices for an upper triangle in column-first order.

Larbino1 · April 20, 2023, 1:18pm

Thank you for this! Perfect for what I need.

ilanggear · June 15, 2023, 11:44pm

I don’t need an iterator, I need a direct formula mapping a linear Index of a triangular matrix into cartesian coordinates. After reading various sources on the web like this one I realized that the formula will vary based on quadrant, 0- vs 1-based indexing, row vs column major, and whether the diagonal was included or not. So none of the explicit formulas matched my particular use case (upper right quadrant, 1-based indexing, column major ordering, without the diagonal). I was able to derive a formula to map from cartesian coordinates to a linear index (Index = Int(x*(x-1)/2 -x +y +1) ), but I’m having trouble doing the reverse. This pair of formulas are close, but include the diagonal:

i = Int(ceil(sqrt(2 * idx + 0.25) - 0.5))
j = Int(idx - (i-1) * i / 2)

Any help would be much appreciated!

rocco_sprmnt21 · June 16, 2023, 7:27am

I tried to remove the diagonal from your function, although I don’t understand the characteristic of the upper right quadrant

function urq1cmowd(idx)
r=findfirst(n->n*(n+1)/2>idx, 1:idx)+1
for i=2:r, j=1:i-1
    println((i,j))
end
end
urq1cmowd(11)

There is obviously to fix the output to filter only the first idx elements eventually

function urq1cmowd(idx)
r=findfirst(n->n*(n+1)/2>idx, 1:idx)+1
for i=2:r, j=1:i-1
    (i-1)*(i-2)/2+j<=idx ? println((i,j)) : break
end
end

rocco_sprmnt21 · June 16, 2023, 10:20am

It should not be iuppert(3,3) == (2,3) ?

A slightly different alternative proposal

function urq1cmowd(idx)
    #M=round(Int,sqrt(2*idx))
    #rnk=findfirst(n->n*(n+1)/2>idx, 1:M)
    rnk=findfirst(n->n*(n+1)/2>idx, 1:idx)
    res=Vector{Tuple{Int64, Int64}}(undef,idx)
    k=1
    for i=1:rnk, j=i+1:rnk
        (i-1)*(i-2)/2+j<=idx ? res[k]=(i,j) : break
        k+=1
    end
    res
end
using BenchmarkTools
julia> @btime urq1cmowd(10)
  42.525 ns (1 allocation: 224 bytes)
10-element Vector{Tuple{Int64, Int64}}:
 (1, 2)
 (1, 3)
 (1, 4)
 (1, 5)
 (2, 3)
 (2, 4)
 (2, 5)
 (3, 4)
 (3, 5)
 (4, 5)

julia> @btime iuppert.(1:10,5)
  52.632 ns (1 allocation: 224 bytes)
10-element Vector{Tuple{Int64, Int64}}:

lmiq · June 16, 2023, 11:44am

Yes, actually it does that, I don’t know what I copied there:

julia> @inline function iuppert(k::Integer,n::Integer)
         i = n - 1 - floor(Int,sqrt(-8*k + 4*n*(n-1) + 1)/2 - 0.5)
         j = k + i + ( (n-i+1)*(n-i) - n*(n-1) )÷2
         return i, j
       end
iuppert (generic function with 1 method)

julia> iuppert(1,3)
(1, 2)

julia> iuppert(2,3)
(1, 3)

julia> iuppert(3,3)
(2, 3)

(I’ll fix the post, if possible still)

lmiq · June 16, 2023, 11:45am

I think the formula there is actually for this case.

rocco_sprmnt21 · June 16, 2023, 5:32pm

after having better read the specifics of the problem, I too tried to derive a formula for (i,j).
I ended up with the same result used in the solution.
I leave my notes as an insight on the demonstration process

ilanggear · June 16, 2023, 6:05pm

The first two urq1cmowd functions provide the correct answers, but I thought a purely formulaic approach was possible. I need something efficient because I’m processing a data structure with millions (if not billions) of elements and I don’t have any memory to spare. (Btw, the definition of the upper right quadrant I’m using is simply x ranging from 1:n and j ranging from 1:i-1). Here is a sample of coordinates:

(j, i) <=> k
(2, 1) <=> 1
(3, 1) <=> 2
(3, 2) <=> 3
(4, 1) <=> 4
(4, 2) <=> 5
(4, 3) <=> 6
(5, 1) <=> 7
(5, 2) <=> 8
(5, 3) <=> 9
(5, 4) <=> 10
(6, 1) <=> 11

ilanggear · June 16, 2023, 6:11pm

What I need:

x  1  2  4  7 11 16
x  x  3  5  8 12 17
x  x  x  6  9 13 18
x  x  x  x 10 14 19
x  x  x  x  x 15 20
x  x  x  x  x  x 21

Generated by “isuppert” (row-major ordering, in contrast column major ordering like I need above):

x  1  2  3  4  5  6
x  x  7  8  9 10 11
x  x  x 12 13 14 15
x  x  x  x 16 17 18 
x  x  x  x  x 19 20
x  x  x  x  x  x 21

rocco_sprmnt21 · June 16, 2023, 7:13pm

try this

function ut2ci(idx)
    j = Int(ceil(sqrt(2 * idx + 0.25) - 0.5))+1
    i = Int(idx - (j-2) * (j-1) / 2)
    CartesianIndex(i,j)
end

ci=ut2ci.(1:21);
m=zeros(Int,7,7);
m[ci].=1:21;
m

julia> m
7×7 Matrix{Int64}:
 0  1  2  4   7  11  16
 0  0  3  5   8  12  17
 0  0  0  6   9  13  18
 0  0  0  0  10  14  19
 0  0  0  0   0  15  20
 0  0  0  0   0   0  21
 0  0  0  0   0   0   0

I don’t know if, apart from “aesthetic” issues, for arrays like these, it’s preferable to use SparseArrays, since you already have the coordinates of the non-zero values.

julia> m=zeros(Int,7,7);

julia> using SparseArrays

julia> sm=sparse(m)
7×7 SparseMatrixCSC{Int64, Int64} with 0 stored entries:
 ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅

julia> sm[ci].=1:21;

julia> sm
7×7 SparseMatrixCSC{Int64, Int64} with 21 stored entries:
 ⋅  1  2  4   7  11  16
 ⋅  ⋅  3  5   8  12  17
 ⋅  ⋅  ⋅  6   9  13  18
 ⋅  ⋅  ⋅  ⋅  10  14  19
 ⋅  ⋅  ⋅  ⋅   ⋅  15  20
 ⋅  ⋅  ⋅  ⋅   ⋅   ⋅  21
 ⋅  ⋅  ⋅  ⋅   ⋅   ⋅   ⋅

Per · June 16, 2023, 7:18pm

using TriangularIndices
ix = UpperTriangularIndices(6)
z = zeros(Int, (6,7))
for k in 1:21
      z[(ix[k] .+ (0,1))...] = k
end
z

should do the trick.

ilanggear · June 16, 2023, 7:47pm

Yes!! This does what I need. Thank you so much!!

ilanggear · June 16, 2023, 7:49pm

Thanks for the reply, but I can’t use something that allocates memory like the call to zeros because the sizes that I need to deal with are too large.