Julia: Broadcast with LowerTriangular, UpperTriangular, and Diagonal gives Diagonal

[sethaxen]

Broadcasting addition for UpperTriangular, LowerTriangular, and Diagonal produces a Diagonal. Note that normal addition is fine, and if either of these matrices are excluded, the result is fine.

julia> using Random, LinearAlgebra
julia> Random.seed!(42);
julia> U = UpperTriangular(randn(3,3));
julia> L = LowerTriangular(randn(3,3));
julia> D = Diagonal(randn(3));
julia> U+L+D
3×3 Array{Float64,2}:
 0.857936  -0.299484  -0.468606
 1.08238    2.41289    0.156143
 0.187028   0.367563  -5.28356 
julia> U.+L.+D
3×3 Diagonal{Float64,Array{Float64,1}}:
 0.857936   ⋅         ⋅     
  ⋅        2.41289    ⋅     
  ⋅         ⋅       -5.28356

Version info:

julia> versioninfo()
Julia Version 1.3.0-rc2.0
Commit a04936e3e0 (2019-09-12 19:49 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin18.6.0)
  CPU: Intel(R) Core(TM) i7-4870HQ CPU @ 2.50GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, haswell)

Also confirmed on 1.2.

[andreasnoack]

That is a pretty bad bug. The bug was introduced between 0.6 and 0.7. @mbauman Any idea why the structure detection is wrong here.

[mbauman]

Oh yikes. So the issue is three-fold:

1. Broadcasts that include both upper- and lower-triangular matrices promote to a default matrix style.
2. Broadcasts that combine a default matrix style with some other structured matrix attempt to preserve the structure
3. We assume that if all arguments are structured matrices, then the promoted style describes all the locations where nonzero may fall.

Clearly assumption 3 is not true in the presence of points 1 and 2. Not sure how to fix this yet...

[goggle]

I have spent quite some time now to investigate this bug.

One possible solution that came to mind is to introduce a new StructuredMatrix type StructuredMatrix{Matrix}, so that combinations of lower and upper triangular matrices can promote to this type, which should (hopefully) solve this issue. Any opinions on that?

I'm currently trying to understand, what the function fzeropreserving in stdlib/LinearAlgebra/structuredbroadcast.jl exactly does... It would be great if anyone could explain this to me, because I think there lies another possibility to solve this issue.

[timholy]

Yes, I think StructuredMatrix{Matrix} is the right way to go here.

I'm currently trying to understand, what the function fzeropreserving in stdlib/LinearAlgebra/structuredbroadcast.jl exactly does

I don't think fzeropreserving is the way to go, use your styles approach. (For reference, fzeropreserving tests whether fzero(mtrx) will return 0. An array-of-arrays is an example of something for which that would not be true.)

[mbauman]

Yes, doing a StructuredMatrix{Matrix} seems like a very reasonable plan of attack!

The pair of functions isstructurepreserving and fzeropreserving test to see if structured zeros will remain zero for a given structured broadcast style. Both are overly conservative — they may return false negatives but should never return a false positives. isstructurepreserving is the most conservative and done entirely on the type domain. fzeropreserving will actually evaluate f (and is therefore not type stable) — but only in cases where we know what the values will be for the structured zeros (e.g., scalars and structured matrices).

That's why they don't really help here — both assume that the all values outside of the style type will be structured zeros for structured matrices.

[goggle]

Thanks @timholy and @mbauman, I will try this. If it succeeds, I will set up a PR to solve this issue.