Matrix versus matrix transpose

PetrKryslUCSD · January 18, 2018, 6:05pm

I’m a bit confused: how is the matrix transpose expected to work?

Let’s say I form a transpose:

julia> a = rand(3,3)
3×3 Array{Float64,2}:
 0.138434  0.0714774  0.570142
 0.169383  0.0834496  0.1411
 0.704381  0.0703364  0.0379085

julia> transpose(a)
3×3 Transpose{Float64,Array{Float64,2}}:
 0.138434   0.169383   0.704381
 0.0714774  0.0834496  0.0703364
 0.570142   0.1411     0.0379085

But in my mind that transpose IS an ARRAY, so I should be able to pass it to a function that expects an Array{Float64,2}. However that function is not going to accept that as an argument, because now I’m going to pass a Transpose{Float64,Array{Float64,2}} instead of a Array{Float64,2}.

Am I missing something?

(Sorry, I forgot to say that I’m talking about 0.7.)

rdeits · January 18, 2018, 6:23pm

A Transpose is still an AbstractArray, so you can pass it to any function that expects an AbstractArray. If a function expects exactly an Array{Float64, 2} then either (1) that function is doing something that relies on the exact internal representation of an Array{Float64, 2}, in which case you can copy() the transpose to get a “plain” array, or (2) that function’s input type signatures are too restrictive, in which case the function’s input types should be updated.

PetrKryslUCSD · January 18, 2018, 6:49pm

OK, that makes sense. Thanks.

ChrisRackauckas · January 18, 2018, 7:38pm

Don’t type your functions to ::Array. Instead type them to ::AbstractArray. And make your types have type parameters instead of ::Array. It makes this all a non-issue without changing performance, and gives you a lot of new features for free.

PetrKryslUCSD · January 18, 2018, 7:50pm

That is good advice.

turtle · January 18, 2018, 8:23pm

As I understand this is a version 0.7 behaviour,
because for version 0.6 transpose(a) returns a 3×3 Array{Float64,2}

I do not find any trace of this change in the release notes for the master.
Could you point me where the change happened.

ChrisRackauckas · January 18, 2018, 9:37pm

github.com/JuliaLang/julia

Taking matrix transposes seriously

opened 08:54PM - 10 Mar 17 UTC

closed 02:43AM - 15 Jan 18 UTC

StefanKarpinski

decision breaking linear algebra

Currently, `transpose` is recursive. This is pretty unintuitive and leads to thi…s unfortunateness: ```jl julia> A = [randstring(3) for i=1:3, j=1:4] 3×4 Array{String,2}: "J00" "oaT" "JGS" "Gjs" "Ad9" "vkM" "QAF" "UBF" "RSa" "znD" "WxF" "0kV" julia> A.' ERROR: MethodError: no method matching transpose(::String) Closest candidates are: transpose(::BitArray{2}) at linalg/bitarray.jl:265 transpose(::Number) at number.jl:100 transpose(::RowVector{T,CV} where CV<:(ConjArray{T,1,V} where V<:(AbstractArray{T,1} where T) where T) where T) at linalg/rowvector.jl:80 ... Stacktrace: [1] transpose_f!(::Base.#transpose, ::Array{String,2}, ::Array{String,2}) at ./linalg/transpose.jl:54 [2] transpose(::Array{String,2}) at ./linalg/transpose.jl:121 ``` For some time now, we've been telling people to do `permutedims(A, (2,1))` instead. But I think we all know, deep down, that's terrible. How did we get here? Well, it's pretty well-understood that one wants the `ctranspose` or "adjoint" of a matrix of matrices to be recursive. A motivating example is that you can [represent complex numbers using 2x2 matrices](https://en.wikipedia.org/wiki/Conjugate_transpose#Motivation), in which case the "conjugate" of each "element" (actually a matrix), is its adjoint as a matrix – in other words, if ctranspose is recursive, then everything works. This is just an example but it generalizes. The reasoning seems to have been the following: 1. `ctranspose` should be recursive 2. `ctranspose == conj ∘ transpose == conj ∘ transpose` 3. `transpose` therefore should also be recursive I think there are a few problems here: * There's no reason that `ctranspose == conj ∘ transpose == conj ∘ transpose` has to hold, although the name makes this seem almost unavoidable. * The behavior of `conj` operating elementwise on arrays is kind of an unfortunate holdover from Matlab and isn't really a mathematically justifiable operations, much as `exp` operating elementwise isn't really mathematically sound and what `expm` does would be a better definition. * It is actually taking the conjugate (i.e. adjoint) of each element that implies that `ctranspose` should be recursive; in the absence of conjugation, there's no good reason for transpose to be recursive. Accordingly, I would propose the following changes to remedy the situation: 1. Rename `ctranspose` (aka `'`) to `adjoint` – that is really what this operation does, and it frees us from the implication that it must be equivalent to `conj ∘ transpose`. 2. Deprecate vectorized `conj(A)` on arrays in favor of `conj.(A)`. 3. Add a `recur::Bool=true` keyword argument to `adjoint` (née `ctranspose`) indicating whether it should call itself recursively. By default, it does. 4. Add a `recur::Bool=false` keyword argument to `transpose` indicating whether it should call itself recursively. By default, it does not. At the very minimum, this would let us write the following: ```jl julia> A.' 4×3 Array{String,2}: "J00" "Ad9" "RSa" "oaT" "vkM" "znD" "JGS" "QAF" "WxF" "Gjs" "UBF" "0kV" ``` Whether or not we could shorten that further to `A'` depends on what we want to do with `conj` and `adjoint` of non-numbers (or more specifically, non-real, non-complex values). [This issue is the second of an ω₁-part series.]

Transposes are lazy now.

mbauman · January 18, 2018, 10:05pm

The PR is here: https://github.com/JuliaLang/julia/pull/25364

I just flagged it as needing a NEWS update.

glemieux · February 3, 2019, 6:11am

I just found this conversation when trying to understand why I couldn’t pass a transposed array to lyap. Why shouldn’t it be able to accept a transpose or adjoint of an array?

julia> A = rand(2,2)
2×2 Array{Float64,2}:
 0.14733  0.885718
 0.79308  0.426695

julia> Q = I+zeros(2,2)
2×2 Array{Float64,2}:
 1.0  0.0
 0.0  1.0

julia> lyap(A,Q)
2×2 Array{Float64,2}:
  0.44532   -0.638588 
 -0.638588   0.0151193

julia> lyap(transpose(A),Q)
ERROR: MethodError: no method matching lyap(::Transpose{Float64,Array{Float64,2}}, ::Array{Float64,2})
Closest candidates are:
  lyap(::Union{DenseArray{T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64},2}, ReinterpretArray{T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64},2,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, ReshapedArray{T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64},2,A,MI} where MI<:Tuple{Vararg{SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, SubArray{T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64},2,A,I,L} where L where I<:Tuple{Vararg{Union{Int64, AbstractRange{Int64}, AbstractCartesianIndex},N} where N} where A<:Union{ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, DenseArray}}, ::Union{DenseArray{T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64},2}, ReinterpretArray{T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64},2,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, ReshapedArray{T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64},2,A,MI} where MI<:Tuple{Vararg{SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, SubArray{T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64},2,A,I,L} where L where I<:Tuple{Vararg{Union{Int64, AbstractRange{Int64}, AbstractCartesianIndex},N} where N} where A<:Union{ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, DenseArray}}) where T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64} at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/dense.jl:1476
Stacktrace:
 [1] top-level scope at none:0

rdeits · February 3, 2019, 6:57am

From glancing at the source code, the implementation of lyap seems to rely on some LAPACK routines which probably expect strided matrices: julia/dense.jl at 69ac379c9d38f758b32e769711c8b1ac90afa25d · JuliaLang/julia · GitHub

You should be able to do lyap(copy(transpose(A)), Q) to make this work in your case.