Hermitian, ishermitian now check diagonals

[jiahao]

Fix #10671

I also had to fix some tests which had incorrectly neglected to check the diagonals.

[andreasnoack]

I think we should do this without allocation.

[jiahao]

Can anyone shed light on the 32-bit build failure?

julia-debug: ast.c:686: jl_lam_vinfo: Assertion (((jl_value_t*)(((jl_taggedvalue_t *) ((char )((l)) - __builtin_offsetof (jl_taggedvalue_t, value)))->type_bits&~(size_t)3))==(jl_value_t)(jl_expr_type))' failed.`

[jiahao]

Build failure seems unrelated.

[tkelman]

This is causing intermittent failures, the tolerances look a little too tight

[jiahao]

Could you point me to a failure?

[tkelman]

Here's 4, haven't even checked appveyor.
https://travis-ci.org/JuliaLang/julia/jobs/57350258
https://travis-ci.org/JuliaLang/julia/jobs/57350849
https://travis-ci.org/JuliaLang/julia/jobs/57354524
https://travis-ci.org/JuliaLang/julia/jobs/57354525

[tkelman]

Bump. Any suggestions for mitigating these?

[jiahao]

I wrote a simple PR but I got an inscrutable error.

[andreasnoack]

I think we should revert because it gives some annoying exceptions this are often unwanted. We don't check that a matrix is symmetric in cholfact so I don't think we should check the input of Hermitian. We should simply ignore the imaginary part of the diagonal and document the behavior.

cc: @bnels

[tkelman]

Ignoring the diagonal sounds like a pretty bad idea - the constructors here should verify the invariants that their types represent.

[andreasnoack]

No. The point is that the user already knows that the matrix has the property or that the matrix should be interpreted in that way. Similarly to allowing

julia> LowerTriangular(randn(3,3))
3x3 LowerTriangular{Float64,Array{Float64,2}}:
 -0.193406   0.0       0.0    
 -0.225881  -1.29177   0.0    
  1.83933    1.77965  -1.24135

[tkelman]

Unlike LowerTriangular, you aren't proposing ignoring uninitialized data - you're proposing ignoring actual imaginary parts of values of elements that matter. You'd have to zero them and make the Hermitian constructor mutating, which is surprising and inconsistent with the rest of the wrapper types.

[andreasnoack]

It's not necessary to mutate anything. The methods for Hermitian should just ignore any imaginary components in the diagonal. That is also what LAPACK is doing and what we are doing for Cholesky

julia> cholfact(fill(complex(1.0,1.0), 1, 1))
Base.LinAlg.Cholesky{Complex{Float64},Array{Complex{Float64},2}} with factor:
1x1 UpperTriangular{Complex{Float64},Array{Complex{Float64},2}}:
 1.0+0.0im

[tkelman]

Every single method then. getindex and everything else you might ever call on a Hermitian will need extra special casing to ignore the imaginary component any time it touches a diagonal. Sounds very very easy to get wrong for not much benefit, easier to explicitly verify up front.

[andreasnoack]

No. It's not easy to verify without a lot of unnecesary exceptions. One of the main reasons to use Hermitian is when floating point errors results in a computation like A = QΛQ' begin slightly non-hermitian. Then you want to use Hermaitan(A), but it defeats the purpose to throw an exception in that case.

I dont' think ingoring imaginary parts in the diagonal will be that demanding. Making the very small change to getindex will probably take care of most of it.

[tkelman]

tril and friends as well, copytri!, probably others. And when you send the data to a C library?

[andreasnoack]

There might be a few other places in symmetric.jl, but the file is only 239 lines so it's very doable. As far as I can tell, we are calling out to three libraries with Hermitian: BLAS, LAPACK and CHOLMOD. The first two use the assumption I advocate for here and CHOLMOD checks the imaginary part like we are doing now.

[tkelman]

Methods that use in any way other than getindex, or leave alone without checking, diagonal elements for any abstract supertype of Hermitian would need to be made aware of this special case, or have a more specific method written for Hermitian. I'm anticipating there will be quite a few places in base that would get this wrong. This is the kind of thing where we start needing type-based coverage stats.

I'm more referring to users sending their data to C libraries outside of base. The proposal is to store undefined data that we intentionally want to be meaningless inline within the same diagonal scalar values. This could be documented and tested thoroughly for corner cases, and left up to the user to manually zero the imaginary components for this ccall case, but this isn't something we've done any other places to my knowledge.