minor fixes in multiplication with Diagonals

[dkarrasch]

• there is no *mul! with vectors
• order of multiplication
• removed todo note, since that has been resolved already

[StefanKarpinski]

These seem like significant bug fixes. Would you mind adding some tests that would have caught them?

[dkarrasch]

Sure. BTW, I'd appreciate if anyone could explain to me when/where Adjoint/Transpose{<:Any, <:Diagonal} are generated. The type Diagonal is closed under adjoint and transpose, and the Adjoint wrapper only arises when calling the constructor Adjoint(::Diagonal) explicitly. Do we often call the constructor (instead of adjoint or ' or transpose) in other places in LinearAlgebra? I wonder why we have so many methods for this "strange" combination.

[dkarrasch]

I realized a couple of things. In my first commit, I overlooked the recursive action of Adjoint and Transpose, so those versions were actually correct, but not tested. Second, in rmul!(A::AbstractMatrix, D::Diagonal) the transpose was superfluous in the definition of broadcast, since it transposed also the elements in the case of diagonal matrix of matrices. What we really want there is to permutedim. I don't think this was tested with matrix elements, but now is. Third, I continue to be confused about the possibility to both adjoint(D) and Adjoint(D). The first creates a new vector with adjoint elements wraps the diag vector into an Adjoint and constructs a Diagonal out of that. The second simply creates an Adjoint wrapper around the Diagonal. IIUC, we don't recommend to users to construct Adjoints directly, but to use adjoint. Now, adjoint is subject to multiple dispatch, and we seem to have chosen deliberately to not wrap the Diagonal by an Adjoint, but to create another Diagonal. (Similarly, we don't wrap Hermitian matrices when we adjoint them, but return the Hermitian itself.)

IMO, we should discourage to wrap a Diagonal by an Adjoint, either by throwing an error, or to define Adjoint(D::Diagonal) = adjoint(D). We would then no longer need methods for Adjoint{<:Any, <:Diagonal}, because this was unconstructable. Now that we can still wrap, I defined

function LinearAlgebra.lmul!(adjA::Adjoint{<:Any,<:Diagonal}, B::AbstractMatrix)
    A = adjA.parent
    return lmul!(adjoint(A), B)
end

which is quite weird, because what we do here is to first "correct" the wrong wrapper, and then multiply as usual, with the caveat that this is not allocation-free, which is expected for adjoint, but not for lmul!. Before, we had

function LinearAlgebra.lmul!(adjA::Adjoint{<:Any,<:Diagonal}, B::AbstractMatrix)
    A = adjA.parent
    return lmul!(conj(A.diag), B)
end

which would give unexpected/"wrong" results in the matrix of matrix case (and it threw because conj(A.diag) is a vector).

[StefanKarpinski]

@andreasnoack, @fredrikekre — can one of you comment on whether this is "just a bug fix" and therefore backportable?

[StefanKarpinski]

Does this change the behavior? Was this previously non-recursive?

[dkarrasch]

For Number eltypes this does not change behaviour, since transpose doesn't do anything to them. For block matrices, this is a bugfix. It makes this Transpose line consistent with the Adjoint line 552 above. Both these are now also covered by tests, testing Complex numbers (where the adjoint boils down to conjugation on the element level) and 2x2 matrices (where both adjoint and transpose have an effect).

[mbauman]

Here's why I call this a bugfix and think it's worthy of a backport:

julia> A = reshape([[1 2; 3 4], zeros(Int,2,2), zeros(Int, 2, 2), [5 6; 7 8]], 2, 2)
2×2 Array{Array{Int64,2},2}:
 [1 2; 3 4]  [0 0; 0 0]
 [0 0; 0 0]  [5 6; 7 8]

julia> adjoint(1:2) * A
1×2 Adjoint{Adjoint{Int64,Array{Int64,2}},Array{Array{Int64,2},1}}:
 [1 2; 3 4]  [10 12; 14 16]

julia> transpose(1:2) * A
1×2 Transpose{Transpose{Int64,Array{Int64,2}},Array{Array{Int64,2},1}}:
 [1 2; 3 4]  [10 12; 14 16]

julia> adjoint(1:2) * Diagonal(A)
1×2 Adjoint{Adjoint{Int64,Array{Int64,2}},Array{Array{Int64,2},1}}:
 [1 2; 3 4]  [10 12; 14 16]

julia> transpose(1:2) * Diagonal(A)
1×2 Transpose{Transpose{Int64,Array{Int64,2}},Array{Array{Int64,2},1}}:
 [1 3; 2 4]  [10 14; 12 16]

[fredrikekre]

LGTM, just one comment.

[andreasnoack]

Thanks for fixing this. I also think this should be considered a bug fix.

I'd appreciate if anyone could explain to me when/where Adjoint/Transpose{<:Any, <:Diagonal} are generated.

I don't recall if we generate these anywhere but I have an idea why they might be there. The problem is that the purpose of Adjoint/Transpose is laziness but applying adjoint/transposeto the diagonal vector aDiagonalwill create a new vector unless we can complete avoid the operation such as when it'sAbstractVector{<:Real}`. Hence, the extra allocation in

julia> Dr = Diagonal(ones(3));

julia> Dc = Diagonal(complex.(ones(3)));

julia> @btime Dr'*Dr.diag;
  93.711 ns (2 allocations: 128 bytes)

julia> @btime Dc'*Dc.diag;
  136.575 ns (3 allocations: 272 bytes)

I'm not sure what the best solution is since it is generally best to push the wrapping to the "simples" array type. Maybe we should consider reintroducing the lazy array conjugation.

[andreasnoack]

@dkarrasch There is a test failure in the diagonal tests on 32bit. Could you take a look and identify the issue? Except for that, I think this one is good to go.

[dkarrasch]

It says that there is an issue with this

Expression: x' * D == x' * Array(D) == copy(x') * D == copy(x') * Array(D)
 Evaluated: 
Complex{Float64}[0.28264113873776175 - 0.08940514975097541im 0.1389278716409064 + 0.23843302149008275im … 0.25966490344645743 - 0.23196007257190188im 0.9204553072599387 - 0.06755870535805081im] == 
Complex{Float64}[0.28264113873776175 - 0.08940514975097541im 0.1389278716409064 + 0.23843302149008275im … 0.25966490344645743 - 0.23196007257190188im 0.9204553072599387 - 0.06755870535805084im] ==
Complex{Float64}[0.28264113873776175 - 0.08940514975097541im 0.1389278716409064 + 0.23843302149008275im … 0.25966490344645743 - 0.23196007257190188im 0.9204553072599387 - 0.06755870535805081im] ==
Complex{Float64}[0.28264113873776175 - 0.08940514975097541im 0.1389278716409064 + 0.23843302149008275im … 0.25966490344645743 - 0.23196007257190188im 0.9204553072599387 - 0.06755870535805084im]

The displayed numbers look good, except for the last digit in the imaginary part of the last entries (scroll to the very right). Can/should I simply replace all == tests by approximate equality tests? I assume there are slightly different computational routes (sometimes BLAS, sometimes explicit loops), which may yield slightly different results. So it's probably even dangerous to have exact tests which may then fail on different machines/OSs etc.

[dkarrasch]

Ok, I was under the wrong impression that we, when taking the adjoint of a Diagonal, simply wrap D.diag by an adjoint. But we can't do that, because adjoint(D.diag) (which would be lazy), would yield an AbtractMatrix, not AbstractVector, which would then give the wrong Diagonal matrix. A bit involved to figure that out, but now I agree that we currently need both Adjoint and adjoint of D.

[dkarrasch]

Travis looks even worse now, but the failure does not seem to be related to the Diagonal tests.

[mbauman]

Boy, that CI run looks way worse than I've been seeing recently, but the ones I sampled do seem unrelated. I still don't feel great about merging this with so many failures, so I'm going to re-trigger CI — I think things have stabilized over the past two days.

[dkarrasch]

I agree we should be careful, and I don't mind if we wait for another few days if other things are still in the flow.

[nalimilan]

I see failures similar to the ones you posted above when building Julia on AArch64 Fedora. I've bisected them to this commit. Any ideas?

Error in testset LinearAlgebra/diagonal:
Test Failed at /builddir/build/BUILD/julia/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/test/diagonal.jl:453
  Expression: lmul!(transform(D), copy(M)) == transform(Matrix(D)) * M
   Evaluated: Complex{Float64}[0.495913+0.570551im 1.9883+0.829459im … 0.605464+0.0268929im -0.798705+0.141228im; -0.268808-0.935347im 0.150569+0.127275im … 0.9115+0.381105im -1.03646-0.663045im; … ; 0.591381+0.739548im -1.90991+0.323593im … -0.144823-1.2759im 0.426071-0.12725im; 0.0475895+0.0491165im -0.378787-0.801472im … 0.17056+0.285193im 0.0111192-0.100995im] == Complex{Float64}[0.495913+0.570551im 1.9883+0.829459im … 0.605464+0.0268929im -0.798705+0.141228im; -0.268808-0.935347im 0.150569+0.127275im … 0.9115+0.381105im -1.03646-0.663045im; … ; 0.591381+0.739548im -1.90991+0.323593im … -0.144823-1.2759im 0.426071-0.12725im; 0.0475895+0.0491165im -0.378787-0.801472im … 0.17056+0.285193im 0.0111192-0.100995im]

(and so on)
https://koji.fedoraproject.org/koji/getfile?taskID=34930452&volume=DEFAULT&name=build.log&offset=-4000

[dkarrasch]

I can't spot the slightest difference in the printed numbers, so I assume the difference is in the last digit(s). Maybe we should have been slightly more generous with approx. tests? It seems, however, that tests passed on most platforms, and in those which failed, I can't seem to find error messages as you report here. What's the recommendation? Should I make a "relaxation" PR?

[nalimilan]

OK, thanks. That would make sense to me since AFAIK we allow BLAS to return slightly different results depending on many details (e.g. because of the use of SIMD and of variations in the number of threads).

Have you found an explanation about the failures you've seen above? Why have they disappeared?

[dkarrasch]

I can't find any failures that would be related to the tests. So, unfortunately, no, I have no explanation of why there were so many failures at the time.