Unexpected? result of multiplication with adjoint()

#1

On Julia Version 1.0.2 the output of the following code

using LinearAlgebra
unblock(A) = mapreduce(identity, hcat, [mapreduce(identity, vcat, A[:,i]) for i = 1:size(A,2)])
x=Vector{Any}(undef,2)
x[1]=[1+2im 1-im;-3 4+3im]
x[2]=[1+im 0;1 1+im]
y=unblock(x)
println(x)
println(y)
println(adjoint(y)*y)
println(transpose(conj(x))*x)
println(adjoint(x)==transpose(conj(x)))
println(adjoint(x)*x)

is

Any[Complex{Int64}[1+2im 1-1im; -3+0im 4+3im], Complex{Int64}[1+1im 0+0im; 1+0im 1+1im]]
Complex{Int64}[1+2im 1-1im; -3+0im 4+3im; 1+1im 0+0im; 1+0im 1+1im]
Complex{Int64}[17+0im -12-11im; -12+11im 29+0im]
Complex{Int64}[17+0im -12-11im; -12+11im 29+0im]
true
46 + 0im

The correct result is the one obtained by line 9 using the unblocked array y.
The same result is obtained by line 10 using blocked array x.
Although by definition and in Julia adjoint equals conjugate transpose, the result of the last line is scalar.

Multiplication by adjoint works fine is all elements of x are scalars. Could it work correctly in the blocked case, too?

0 Likes

#2

See https://github.com/JuliaLang/julia/pull/27401 for the change and a discussion of this change. The aim was to make dot return a scalar. However, I agree that this is unfortunate would have preferred that the result was a matrix. An example that really makes the issue with the current solution clear is

julia> x
2-element Array{Any,1}:
 Complex{Int64}[1 + 2im 1 - 1im; -3 + 0im 4 + 3im]
 Complex{Int64}[1 + 1im 0 + 0im; 1 + 0im 1 + 1im]

julia> xm = reshape(x, 2, 1)
2×1 Array{Any,2}:
 Complex{Int64}[1 + 2im 1 - 1im; -3 + 0im 4 + 3im]
 Complex{Int64}[1 + 1im 0 + 0im; 1 + 0im 1 + 1im]

julia> x'x
46 + 0im

julia> (xm'xm)[1]
2×2 Array{Complex{Int64},2}:
  17+0im   -12-11im
 -12+11im   29+0im

so for array elements, we treat vectors and one-column matrices completely different.

A usecase I’m sad isn’t possible is to store a matrix as a vector of rowvectors.

julia> X1 = randn(3,2);

julia> X2 = [X1[i,:]' for i in 1:size(X1, 1)];

julia> X3 = reshape(X2, length(X2), 1);

julia> X1'X1
2×2 Array{Float64,2}:
 2.58747   0.733045
 0.733045  0.420716

julia> X2'X2
3.0081843634412744

julia> (X3'X3)[1]
2×2 Array{Float64,2}:
 2.58747   0.733045
 0.733045  0.420716

(this would be really useful for SVector elememts.)

0 Likes