Matrix*Diagonal does not properly check the dimensions

[StephenVavasis]

The following unexpected behavior from Diagonal in Julia caused wrong output from my grad student's code for a period of months. When she finally determined that there was a bug, it took her hours to find it. Please fix this! EIther make the Diagonal work as expected, or else make it generate an error.

julia> [1 0;0 1] * Diagonal([2 3])
2×2 Array{Int64,2}:
 2  0
 0  2

[oscardssmith]

The problem here is that Diagonal([2 3]) doesn't do what you think it does. [2 3] is a matrix, not a Vector, so Diagonal is getting the diagonal of the matrix (which is just [2]). You want Diagonal([2, 3])

[oscardssmith]

Actually there might be a bug here. This means that you are multiplying a 1x1 matrix with a 2x2 matrix, which should be throwing an error.

[cafaxo]

The result of Diagonal([2 3]) is a bit surprising to me as well. I would have expected it to not change the dimension of the matrix and return a 1x2 diagonal matrix. If there has been a discussion about this in the past, I would be interested in the reasoning for this.

[StephenVavasis]

I am saying that, regardless of what the author of Diagonal intended, this behavior of silently giving a wrong answer for Diagonal applied to a row vector is not acceptable. Julia is supposed to 'just work' in the same way that Matlab does.

[oscardssmith]

Row vectors aren't a thing in Julia. There are only Vectors and Matrix. Diagonal is clearly documented to return the diagonal of the matrix.

[andreasnoack]

I agree with @oscardssmith that the behavior of Diagonal(::Matrix) makes sense. The alternative would be to error. However, the multiplication should have thrown an error the same way that

julia> [1 0;0 1]*fill(2, 1, 1)
ERROR: DimensionMismatch("matrix A has dimensions (2,2), matrix B has dimensions (1,1)")

does.

[cafaxo]

Diagonal is clearly documented to return the diagonal of the matrix.

This is not correct. The diagonal of a matrix is not a matrix: It is a list of the diagonal elements of the matrix. Diagonal(::Matrix) is supposed to return a view of a Matrix (not the diagonal of a matrix). In my opinion, a view should not silently change the dimensions of the original matrix.

But in this case, the Diagonal view of the 1x2 matrix gives a 1x1 matrix. Was it a conscious decision to do it this way?

[StephenVavasis]

I will point out that in math papers, it is common to use Diag as the constructor for a diagonal matrix and diag as the operator that returns a vector that is the diagonal entries of a matrix. See e.g. p. 94 of L. Tuncel, "Polyhedral and Semidefinite Programming Methods in Combinatorial Optimization," 2010.

[andreasnoack]

The current behavior is just a short version of

julia> Diagonal(diag([2 3]))
1×1 Diagonal{Int64, Vector{Int64}}:
 2

which isn't ambiguous since [2 3] is a Matrix.

Diagonal(::Matrix) is supposed to return a view of a Matrix

We have chosen a different behavior.

In my opinion, a view should not silently change the dimensions of the original matrix.

The Diagonal matrix type is square.

[StephenVavasis]

It seems that a sequence of design decisions has left a trap in the language. It's difficult for me to discern which individual design decision should be altered, but somehow you need to figure out how to remove this trap. People are going to design nuclear power plants and vaccines with Julia.

[StefanKarpinski]

Hopefully those people do a little bit of testing of their programs before applying them. I agree that this isn't ideal, but please try not to catastrophize this in order to make your case. Since changing this particular behavior would be breaking, it seems unlikely that we can change this before Julia 2.0 (or at least LinearAlgebra 2.0).

[StephenVavasis]

I disagree. A big matrix somewhere in a code come out wrong. How would you test for that? The only way is to re-implement the same algorithm in a different language and then compare. Even when my grad student figured out that the matrix was wrong because of weird inconsistencies between runs, it was still not obvious what was the mistake.

[oscardssmith]

If this correctly errored on multiplication, I think it would have been pretty clear what went wrong. It would give an error message that would inform you that the diagonal matrix wasn't what you expected.

[StephenVavasis]

As near as I can tell, the reason that there is no error message is because right-multiplication by a diagonal matrix is implemented with a . operation, and . is overloaded to do something special if one of the operands has a single entry. So possibly to make this case error out, instead of using .*, right-multiplication by a diagonal matrix could use one of the more specialized broadcast routines like _bcs1 in broadcast.jl? I am putting a question mark because broadcast.jl is really difficult to understand.

Note that multiplying a 3-by-3 matrix times an ordinary 1-by-1 matrix (using *) gives the expected error message, so the lack of the error message in the diagonal case is indeed an inconstency in LinearAlgebra, and fixing it would address my concern.