JuliaLang / julia

The Julia Programming Language
https://julialang.org/
MIT License
45.68k stars 5.48k forks source link

Unexpected transposition when combining convert() and copy() for adjoints #26294

Closed rdeits closed 6 years ago

rdeits commented 6 years ago

I just noticed this while trying to fix WoodburyMatrices.jl for v0.7. Here's a reduced example of the behavior:

julia> struct Woodbury{Vt} <: AbstractMatrix{Float64}
         V::Vt
       end

julia> Woodbury(V::AbstractMatrix) = Woodbury{typeof(V)}(copy(V))

Previously, this worked fine. However, something interesting happens if V happens to be an Adjoint:

julia> V = zeros(3, 2)
3×2 Array{Float64,2}:
 0.0  0.0
 0.0  0.0
 0.0  0.0

julia> Vprime = V'
2×3 LinearAlgebra.Adjoint{Float64,Array{Float64,2}}:
 0.0  0.0  0.0
 0.0  0.0  0.0

julia> Woodbury(Vprime).V
3×2 LinearAlgebra.Adjoint{Float64,Array{Float64,2}}:
 0.0  0.0
 0.0  0.0
 0.0  0.0

Note that even though we constructed Woodbury with the 2x3 Vprime, we got a 3x2 matrix in the V field. What's happening is:

  1. We call Woodbury(Vprime::Adjoint)
  2. typeof(Vprime) is Adjoint
  3. copy(Vprime) returns a plain Matrix
  4. we call Woodbury{Adjoint{...}}(vprime_copy::Matrix)
  5. that implicitly calls convert(Adjoint, vprime_copy::Matrix)...
  6. ...which in turn transposes vprime_copy, yielding a matrix which is transposed from the original Vprime

In essence, this is a case where convert(typeof(V), copy(V)) returns the adjoint of V, not something of the same shape as V.

This is all perfectly correct given the currently defined behaviors, but it seems like an edge case that might be worth fixing. For example, if copy(Vprime) returned another Adjoint, then I believe this particular problem would go away (but I imagine it's not as easy as that).

Of course, the answer may just be "don't do that". It's easy to fix this particular case now that I understand what's going on, but it may be a trap that others fall into.

andreasnoack commented 6 years ago

I think this is a "don't do that" and the copying constructor should be

function Woodbury(V::AbstractMatrix)
    cV = copy(V)
    Woodbury{typeof(cV)}(cV)
end
rdeits commented 6 years ago

That is indeed the fix I went with. Thanks!

JeffBezanson commented 6 years ago

I think this convert method is wrong, as in #26178. The result of convert should be "basically the same thing", just converted to a different type. It shouldn't change values.

rdeits commented 6 years ago

Thank you!