Nice way to iterate through a triangular matrix

[RobertGregg]

I often need to loop through the upper/lower triangular section of a matrix and usually end up doing something like:

A = reshape(1:25,5,5)

for j in axes(A,2)
    for i in axes(A,1)
        if i>j
            @show A[i,j]
        end
    end
end

This is fine but can lead to code with way too many indents. You could get rid of the if statement, but its still very clunky and takes a second to figure out whats going on:

for j in axes(A,2)
    for i in j+1:size(A,1)
        @show A[i,j]
    end
end

I came up with simpler solution and wanted to share it. Basically Iterators does all the heavy lifting.

lowTriIter(ax1, ax2) = Iterators.filter(i -> first(i)>last(i),
                                        Iterators.product(ax1,ax2) )

lowTriIter(A::AbstractMatrix) = lowTriIter(axes(A)...)

Which leads to code that looks like this:

for (i,j) in lowTriIter(A)
    @show A[i,j]
end

It’s pretty versatile like you could change the comparison to include the diagonal and there’s no fussing around with adding one to indices. You could probably even get it to work with CartesianIndices if you wanted. Feel free to suggest improvements

[Dan]

Another way could be:

lowTriIter2(A::AbstractMatrix) = 
  Iterators.flatten(Iterators.map(j->(j,i), i+1:size(A,1)) for i=1:size(A,2))

julia> collect(lowTriIter2(A))
10-element Vector{Tuple{Int64, Int64}}:
 (2, 1)
 (3, 1)
 (4, 1)
 (5, 1)
 (3, 2)
 (4, 2)
 (5, 2)
 (4, 3)
 (5, 3)
 (5, 4)

This avoids looping over the skipped indices, and therefore might save a little time.

[RobertGregg]

Oh I don’t think I would have ever thought of making a generator of maps, interesting. I guess you could skip right to the generator as well:

lowTriIter3(A::AbstractMatrix) =
    ((i,j) for j in axes(A,2) for i in j+1:size(A,1))

I do like that these methods skip over the additional iterations but I do wish they avoided the i+1:size(A,1) bit. Also I made my version more complicated than it need to be. I think I was worried about axes(A) having more than two iterators, but it won’t because A is already a matrix:

lowTriIter(A::AbstractMatrix) = 
    Iterators.filter(i -> first(i)>last(i), Iterators.product(axes(A)...))

[baggepinnen]

@RobertGregg

This is fine but can lead to code with way too many indents.

You know you can write both iterations on one line?

for j in axes(A,2), i in j+1:size(A,1)
    @show A[i,j]
end

[rocco_sprmnt21]

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

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], Carte…

[RobertGregg]

I did, but that’s a fair point! It’s wild how many ways you can accomplish the same task.

I think I liked my version because it abstracted away the need to think about indices. Like if I wanted to change this nested loop to iterate over the upper triangle it would take me a minute to figure out versus just flipping a less-than sign. I get that’s just a preference though

[RobertGregg]

Ah i didn’t see this. Sorry for the repeated topics. That formula for the ut2ci(idx) function is both amazing and headache inducing

[bertschi]

Just some additional thoughts …

First, you could also use iterator comprehensions and the definition of your iterator looks just like your for loop:

lowertriag(A::AbstractMatrix) = ((i,j) for j in axes(A, 2) for i in axes(A, 1) if i > j)

for (i, j) in lowertriag(A)
    @show A[i, j]
end

Further, flatmap, i.e., flatten ∘ map, is a very powerful combinator and most other operations on iterables can be defined in terms of it:

mymap(f, xs) = Iterators.flatmap(x -> (f(x),), xs)
myfilter(pred, xs) = Iterators.flatmap(x -> if pred(x); (x,) else empty((x,)) end, xs)

In the functional programming jargon a la Haskell, flatmap is the monadic bind whereas map is the way more general functor instance for lists.

[mkitti]

julia> (A[ci] for ci in CartesianIndices(A) if first(ci.I) ≤ last(ci.I)) |> collect
15-element Vector{Int64}:
  1
  6
  7
 11
 12
 13
 16
 17
 18
 19
 21
 22
 23
 24
 25

julia> (A[ci] for ci in CartesianIndices(A) if first(ci.I) ≥ last(ci.I)) |> collect
15-element Vector{Int64}:
  1
  2
  3
  4
  5
  7
  8
  9
 10
 13
 14
 15
 19
 20
 25

[RobertGregg]

Oh I did see that Iterators.flatmap was new to v1.9. It definitely looks useful but I’m not sure how I would apply it to this situation.