Vector adjoint and indexing don't commute?

General Usage
1

I was surprised that vector adjoint doesn’t commute with indexing.

x = [1:5;]
x[1:3]' # is a 1x3 Adjoint
x'[1:3] # is a 3 element Vector

For example if I set y = x' my mental model of y is “row vector”, and I was thinking that y[1:3] * x[1:3] would be a dot product. Instead I must write y[1:3]' * x[1:3] which feels like I’m transposing a second time. What’s the reason for indexing not maintaining the adjointness?

In a similar vein perhaps, reverse:

reverse(x)' # is a 1x5 Adjoint
reverse(x') # is an error for not specifying which dimension

If I’m reversing a vector adjoint, would it be sensible for reverse to default to operating along the vector?

1 Like
2

instead of y[1:3] write y[:, 1:3], the reason is that y[1:3] == y[1, 1:3] is just one row which will give you a vector, this is true for any-shape matrix/array

4 Likes
3

This seems to be a common pitfall for newcomers to Julia. Keeping in mind that your mental model of a row vector is actually a matrix will probably help a lot!

4 Likes
4

Right, but I thought the motivation for the Adjoint type is that it’s not a matrix. It’s a map that operates on vectors to produce scalars. So [1,2,3]' * [1,2,3] is a dot product. But [1 2 3] * [1,2,3] is a matrix-vector product, and the result is a 1-element vector. Maybe this is part of the confusion. As far as indexing is concerned, we do still need to think of Adjoints as if they were matrices, although they don’t operate as matrices do.

5

The conflict here is that a single index also serves as the linear index into a container (which, by convention, always returns a Vector). Thus for an adjoint, this conflicts with the “normal” index and apparently the linear one wins. I’m not sure there is a way to do a “normal” index into a adjoint vector.

3 Likes
6

Ah yes good point, it’s the same syntax for both, so for a “morally one-dimensional” array you can’t distinguish the intention.

7

This will let A[1,:] return a row vector. In case it is useful.

julia> import Base: getindex

julia> getindex(A::AbstractArray,i::Int,j) = adjoint(A'[j,i])

julia> A = randn(4,4)

julia> A[1,:]
1×4 LinearAlgebra.Adjoint{Float64,Array{Float64,1}}:
 -2.12388  1.72479  -0.704399  0.183572

julia> A[1,2:3]
1×2 LinearAlgebra.Adjoint{Float64,Array{Float64,1}}:
 1.72479  -0.704399

julia> A[1,:]*A[:,1]
6.601714378026813
8

This is type piracy, since you don’t own any of the types involved. This will break all kinds of routines in Base, because A[1, :] is always assumed to return an AbstractVector, so use at your own risk.

3 Likes
9

Yes agree it is dangerous. Adding this method breaks the VecOrMat of Vectors testset in matmul.jl.

10

Having A[1,:] return a row vector sounds appealing at first, but it doesn’t generalise to indexing in higher dimensions. The issue was more that an actual vector adjoint was losing its adjointness upon indexing that puzzled me. But now I see that’s an inevitable consequence of using the same square bracket syntax for two distinct forms of indexing: linear indexing and Cartesian indexing. For most array types you can tell the two apart by the number of indices provided, but not for vector adjoints, hence the conflict.

3 Likes
11

Actually digging a little further it doesn’t seem quite as simple as linear indexing always wins. Consider OffsetArrays:

using OffsetArrays
offmat = OffsetArray([1 2 3 ; 4 5 6], 0:1, -1:1)
offvec = OffsetArray([1,2,3], -1:1)

offmat[1:3] # is fine, this is linear indexing
offvec[1:3] # is not fine, this is the array's custom indexing; no syntax for ordinary 1-based linear indexing just because I'm 1-dimensional?

This seems a bit at odds with the documentation on this point:
“When exactly one index i is provided, that index no longer represents a location in a particular dimension of the array. Instead, it selects the ith element using the column-major iteration order that linearly spans the entire array.” i.e. suggesting linear indexing always wins. Or is this strictly describing the base array types, and not applicable to OffsetArray?

It’s water under the bridge now, but it would have been nice to have a separate syntax for 1-based linear indexing whatever the array type. And that would have freed the existing syntax to do the right thing by adjoints. Also it would catch inadvertent under-indexing of two- or higher-dimensional arrays; if you meant linear indexing you’d have to specify it.

2 Likes
12

If there is no reply here, you should open an issue with OffsetArrays.

1 Like
13

I think the thinking is that these really are the linear indices of the offset vector. They just don’t start at 1:

julia> offmat
2×3 OffsetArray(::Array{Int64,2}, 0:1, -1:1) with eltype Int64 with indices 0:1×-1:1:
 1  2  3
 4  5  6

julia> LinearIndices(offmat)
2×3 LinearIndices{2,Tuple{OffsetArrays.IdOffsetRange{Int64,Base.OneTo{Int64}},OffsetArrays.IdOffsetRange{Int64,Base.OneTo{Int64}}}} with indices 0:1×-1:1:
 1  3  5
 2  4  6

julia> LinearIndices(offvec)
3-element LinearIndices{1,Tuple{OffsetArrays.IdOffsetRange{Int64,Base.OneTo{Int64}}}} with indices -1:1:
 -1
  0
  1

julia> parent(offvec)[1:3]
3-element Array{Int64,1}:
 1
 2
 3
2 Likes
14

I always see these posts, but I’m a bit confused why. I almost never carry around 1 x N object in my code. I have no need. What are the use cases?

Honestly you don’t hardly even need to keep track of dimensions anymore. Things like Tullio.jl or TensorOperations.jl just let you write out what you mean. This is the way most people do linear algebra… “Z = A*B*C, or whatever this should be so that you end with the right dimensions for Z”. Or is that just me?

1 Like
15

Oh and just like that, the reverse issue is resolved!
https://github.com/JuliaLang/julia/pull/37367

16

IMO OffsetArrays should not support linear indexing at all, because of the possibility for confusion that you demonstrate. One can always get the parent and index into that if necessary.

17

PRs to straighten out docs confusion would be very welcome. FWIW we also have this:
https://docs.julialang.org/en/v1/devdocs/offset-arrays/#Linear-indexing-(LinearIndices)

Until we have a dedicated “I’m a linear index!” type, there will likely be confusion; but we still need to be able to index OffsetVectors so of course we have to support ::Int index arguments.

3 Likes
18

Also called a covector. See LinearAlgebra and TensorCore future design: Be more greedy · Issue #35763 · JuliaLang/julia · GitHub

3 Likes
19

That explains it very clearly indeed, thanks very much.

20

A linear index can always be converted into a multi-index by using Base._ind2sub to convert.

function Base.getindex(A::Adjoint, I::Int)
    Base.@_inline_meta
    @inbounds getindex(A, Base._to_subscript_indices(A, I)...)
end
Base._to_subscript_indices(A::Adjoint, i::Integer) = (Base.@_inline_meta; Base._unsafe_ind2sub(A, i))
function Base._ind2sub(A::Adjoint, ind)
    Base.@_inline_meta
    Base._ind2sub(axes(A), ind)
end

Thus, the Adjoint should be able to have a linear index.

julia> [1,2,3]'[1]
1

And in fact, it does already have linear indexing built in.