Taking matrix transposes seriously

[StefanKarpinski]

Currently, transpose is recursive. This is pretty unintuitive and leads to this unfortunateness:

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, 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:

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.]

[jiahao]

The logical successor to a previous issue... 👍

[stevengj]

The behavior of conj operating elementwise on arrays is kind of an unfortunate holdover from Matlab and isn't really a mathematically justifiable operations

This is not true at all, and not at all analogous to exp. Complex vector spaces, and conjugations thereof, are a perfectly well-established mathematical concept. See also JuliaLang/julia#19996 (comment)

[StefanKarpinski]

Complex vector spaces, and conjugation thereof, are a perfectly well-established mathematical concept.

Unless I'm mistaken, the correct mathematical conjugation operation in that context is ctranspose rather than conj (which was precisely my point):

julia> v = rand(3) + rand(3)*im
3-element Array{Complex{Float64},1}:
 0.0647959+0.289528im
  0.420534+0.338313im
  0.690841+0.150667im

julia> v'v
0.879291582684847 + 0.0im

julia> conj(v)*v
ERROR: DimensionMismatch("Cannot multiply two vectors")
Stacktrace:
 [1] *(::Array{Complex{Float64},1}, ::Array{Complex{Float64},1}) at ./linalg/rowvector.jl:180

[stevengj]

The problem with using a keyword for recur is that, last I checked, there was a big performance penalty for using keywords, which is a problem if you end up recursively calling these functions with a base case that ends up on scalars.

[StefanKarpinski]

We need to fix the keyword performance problem in 1.0 anyway.

[stevengj]

@StefanKarpinski, you're mistaken. You can have complex conjugation in a vector space without having adjoints — adjoints are a concept that requires a Hilbert space etc., not just a complexified vector space.

Furthermore, even when you do have a Hilbert space, the complex conjugate is distinct from the adjoint. e.g. the conjugate of a complex column vector in ℂⁿ is another complex vector, but the adjoint is a linear operator (a "row vector").

(Complex conjugation does not imply that you can multiply conj(v)*v!)

[StefanKarpinski]

Deprecating vectorized conj is independent of the rest of the proposal. Can you provide some references for the definition of complex conjugation in a vector space?

[stevengj]

https://en.wikipedia.org/wiki/Complexification#Complex_conjugation

(This treatment is rather formal; but if you google "complex conjugate matrix" or "complex conjugate vector" you will find zillions of usages.)

[StefanKarpinski]

If mapping a complex vector to the conjugate vector (what conj does now) is an important operation distinct from mapping it to is conjugate covector (what ' does), then we can certainly keep conj to express that operation on vectors (and matrices, and higher arrays, I guess). That provides a distinction between conj and adjoint since they would agree on scalars but behave differently on arrays.

In that case, should conj(A) call conj on each element or call adjoint? The representation of complex numbers as 2x2 matrices example would suggest that conj(A) should actually call adjoint on each element, rather than calling conj. That would make adjoint, conj and conj. all different operations:

1. adjoint: swap i and j indices and recursively map adjoint over elements.
2. conj: map adjoint over elements.
3. conj.: map conj over elements.

[stevengj]

conj(A) should call conj on each element. If you are representing complex numbers by 2x2 matrices, you have a different complexified vector space.

For example, a commonplace usage of conjugation of vectors is in the analysis of eigenvalues of real matrices: the eigenvalues and eigenvectors come in complex-conjugate pairs. Now, suppose you have a block matrix represented by a 2d array A of real 2x2 matrices, acting on blocked vectors v represented by 1d arrays of 2-component vectors. For any eigenvalue λ of A with eigenvector v, we expect a second eigenvalue conj(λ) with eigenvector conj(v). This won't work if conj calls adjoint recursively.

[stevengj]

(Note that the adjoint, even for a matrix, might be different from a conjugate-transpose, because the adjoint of a linear operator defined in the most general way depends also on the choice of inner product. There are lots of real applications where some kind of weighted inner product is appropriate, in which case the appropriate way to take the adjoint of a matrix changes! But I agree that, given a Matrix, we should take the adjoint corresponding to the standard inner product given by dot(::Vector,::Vector). However, it's entirely possible that some AbstractMatrix (or other linear-operator) types would want to override adjoint to do something different.)

[stevengj]

Correspondingly, there is also some difficulty in defining an algebraically reasonable adjoint(A) for e.g. 3d arrays, because we haven't defined 3d arrays as linear operators (e.g. there is no built-in array3d * array2d operation). Maybe adjoint should only be defined by default for scalars, 1d, and 2d arrays?

Update: Oh, good: we don't define ctranspose now for 3d arrays either. Carry on.

[timholy]

(OT: I'm already looking forward to "Taking 7-tensors seriously," the next installment in the highly-successful 6-part miniseries...)