How to update a symmetric matrix respecting API boundaries?

General Usage
, ,
1

Suppose I have a symmetric matrix H, and I would like to update it in place. For the purposes of this MWE, the update is adding f(i,j) to elements, based on index.

How can I do this while respecting API boundaries? setindex! only allows setting diagonals for symmetric matrices (I am not sure why). I could get the parent, but .uplo is not exposed by LinearAlgebra, if it was, I could do this:

Z = parent(H)
n = size(Z, 1)
if H.uplo == 'U'
    for j in 1:n
        for i in 1:j
            Z[i, j] += f(i, j)
        end
    end
else
    for j in 1:n
        for i in j:n
            Z[i, j] += f(i, j)
        end
    end
end
3 Likes
2 
Z = parent(H)
n = size(Z, 1)
for j in 1:n, i in 1:n
    Z[i, j] += f(i, j)
end

It’s stupid but it’s the only way. The real answer is of course to not respect API boundaries.

EDIT: On a second thought, you can do this

function find_uplo(M)
    temp12 = M[1,2]
    temp21 = M[2,1]
    P = parent(M)
    P[1,2] = zero(P[1,2])
    P[2,1] = one(P[2,1])
    if M[1,2] == P[1,2]
        P[1,2] = temp12
        return 'U'
    else
        P[2,1] = temp21
        return 'L'
    end
end
3

I see the following general solutions to fix this:

  1. allow setindex! to work off-diagnal for Symmetric. The idea is that it automatically does the right thing, eg if the inner representation is upper, and you are trying to set the lower triangle, it just exchanges indices. Cheap (just checking a conditional), convenient, hassle free, does not expose the implementation. (It is understood that if the parent is immutable, this would error). Disadvantage: access pattern may not be optimal (rows vs columns).
  2. Exposing .uplo as part of the API.
  3. Combination of both: setindex! always does the right thing, but you can make it optimal by conditioning on which triangle stores stuff.
1 Like
4

something I have been missing occasionally is a way to iterate the structural indices in an efficient way. For

  • Symmetric this would be the indices of the upper/lower triangles (in the most performant order)
  • Diagonal, it would be equivalent to diagind
  • SparseMatrixCSC it would be equivalent to the nested loops I always need to look up again in the docs of nzrange

Such a mechanism in combination with extending setindex! to offdiagonal indices for Symmetric would be a nice generic interface for this usecase.

2 Likes
5

I think that this would be great to experiment with in a package (that would support relevant types in Base and the standard libraries, and allow other packages to hook into the API).

First I thought that nonzero (sparse, diagonal) and redundant (eg Symmetric) indices would need different handling, but I think that an iterator that would go through all indices that

  1. are not redundant, and
  2. you are allowed to change

would be enough.

6

setindex! for off-diagonals isn’t supported for a Symmetric since A[i,j] and A[j,i] are linked, and setindex! is expected to modify only one index and not both. I think exposing uplo is a better idea, since that allows one to wrap the parent in a triangular matrix and modify the off-diagonal elements. This is a well-defined setindex! operation.

julia> S = Symmetric([1 2; 3 4])
2×2 Symmetric{Int64, Matrix{Int64}}:
 1  2
 2  4

julia> U = UpperTriangular(parent(S))
2×2 UpperTriangular{Int64, Matrix{Int64}}:
 1  2
 ⋅  4

julia> U[1,2] = 10
10

julia> S
2×2 Symmetric{Int64, Matrix{Int64}}:
  1  10
 10   4
7

I think

function Base.setindex!(A::Symmetric, val, i::Integer, j::Integer)
           if (i < j && A.uplo == 'L') || (i > j && A.uplo == 'U')
               setindex!(A.data, val, j, i)
           else
               setindex!(A.data, val, i, j)
           end
           return A
       end

is a perfectly sensible interface, if you didn’t want the matrix to stay symmetric you shouldn’t have wrapped it with Symmetric.

The problem with it is that it will silently send performance to hell, so it is indeed better to expose uplo and let the user handle indexing.

8

The issue with reindexing [i,j] to [j,i] when accessing the wrong triangle is that there are a lot of sensible operations that become easy to mess up. For example, X::Symmetric .+= Y would (without a new specialized implementation) end up adding Y twice to each nondiagonal.

We could try to make specializations that handle these, but I’m not certain we could find something much more sensible or ergonomic than determining the uplo and then modifying just the relevant triangle of the parent. While there are a lot of simple cases where the correct thing would be obvious and trivial, such features usually die under the scrutiny of complex or edge cases.

Without a new interface, I don’t think I’d go much further than to allow setindex! to the “live” triangle and throw a bounds error if accessing the “ghost” half. A new interface might include functions like setindexsym!(with the semantics of setting both [i,j] and [j,i], regardless of matrix type), but at this point it’s still not much better than a manual implementation (assuming a public uplo-getter).

Up to now, if I will do mutations I will tend not to wrap my arrays in Symmetric except at the call sites of functions where the wrapper is relevant.

2 Likes
9
2 Likes