UpperTriangular / LowerTriangular broken on 1.0.0

[penelopeysm]

MWE:

using ForwardDiff, LinearAlgebra
ForwardDiff.gradient(sum, UpperTriangular(randn(3, 3)))

On v0.10.38

julia> ForwardDiff.gradient(sum, UpperTriangular(randn(3, 3)))
3×3 UpperTriangular{Float64, Matrix{Float64}}:
 1.0  1.0  1.0
  ⋅   1.0  1.0
  ⋅    ⋅   1.0

On v1.0.0

julia> ForwardDiff.gradient(sum, UpperTriangular(randn(3, 3)))
ERROR: ArgumentError: cannot set index in the lower triangular part (2, 1) of an UpperTriangular matrix to a nonzero value (Dual{ForwardDiff.Tag{typeof(sum), Float64}}(0.0,0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0))
Stacktrace:
  [1] throw_nonzeroerror(T::Type, x::Any, i::Int64, j::Int64)
    @ LinearAlgebra ~/.julia/juliaup/julia-1.11.4+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/triangular.jl:295
  [2] setindex!
    @ ~/.julia/juliaup/julia-1.11.4+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/triangular.jl:306 [inlined]
  [3] _unsafe_setindex!
    @ ./reshapedarray.jl:297 [inlined]
  [4] setindex!
    @ ./reshapedarray.jl:286 [inlined]
  [5] setindex!
    @ ./subarray.jl:372 [inlined]
  [6] macro expansion
    @ ./broadcast.jl:973 [inlined]
  [7] macro expansion
    @ ./simdloop.jl:77 [inlined]
  [8] copyto!
    @ ./broadcast.jl:972 [inlined]
  [9] copyto!
    @ ./broadcast.jl:925 [inlined]
 [10] materialize!
    @ ./broadcast.jl:883 [inlined]
 [11] materialize!
    @ ./broadcast.jl:880 [inlined]
 [12] seed!(duals::UpperTriangular{…}, x::UpperTriangular{…}, seeds::NTuple{…})
    @ ForwardDiff ~/ppl/ForwardDiff.jl/src/apiutils.jl:52
 [13] vector_mode_dual_eval!(f::typeof(sum), cfg::ForwardDiff.GradientConfig{…}, x::UpperTriangular{…})
    @ ForwardDiff ~/ppl/ForwardDiff.jl/src/apiutils.jl:23
 [14] vector_mode_gradient(f::typeof(sum), x::UpperTriangular{…}, cfg::ForwardDiff.GradientConfig{…})
    @ ForwardDiff ~/ppl/ForwardDiff.jl/src/gradient.jl:91
 [15] gradient
    @ ~/ppl/ForwardDiff.jl/src/gradient.jl:20 [inlined]
 [16] gradient(f::typeof(sum), x::UpperTriangular{…}, cfg::ForwardDiff.GradientConfig{…})
    @ ForwardDiff ~/ppl/ForwardDiff.jl/src/gradient.jl:17
 [17] gradient(f::typeof(sum), x::UpperTriangular{Float64, Matrix{Float64}})
    @ ForwardDiff ~/ppl/ForwardDiff.jl/src/gradient.jl:17
 [18] top-level scope
    @ REPL[2]:1
Some type information was truncated. Use `show(err)` to see complete types.

I think this is caused by #695 and the fact that iszero(ForwardDiff.Dual(0.0, 1.0)) no longer returns true, so the checks here fail https://github.com/JuliaLang/LinearAlgebra.jl/blob/aecb71426201bad8447be446589d3af986886a87/src/triangular.jl#L295-L304

[mcabbott]

It's a bit of an accident that this ever worked. Note that it encodes an assumption that the gradient is to be projected to the subspace defined by UpperTriangular, as the result disagrees with that for collect(x)):

julia> x = UpperTriangular(randn(3, 3));

julia> ForwardDiff.gradient(sum, x)
3×3 UpperTriangular{Float64, Matrix{Float64}}:
 1.0  1.0  1.0
  ⋅   1.0  1.0
  ⋅    ⋅   1.0

julia> ForwardDiff.gradient(sum, collect(x))  # has fewer zeros
3×3 Matrix{Float64}:
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0

julia> ForwardDiff.gradient(x -> sum(UpperTriangular(x)), rand(3,3))  # intermediate UpperTriangular is easier
3×3 Matrix{Float64}:
 1.0  1.0  1.0
 0.0  1.0  1.0
 0.0  0.0  1.0

(@v1.11) pkg> st ForwardDiff
Status `~/.julia/environments/v1.11/Project.toml`
⌅ [f6369f11] ForwardDiff v0.10.38

To preserve this, there might have to be some explicit method like seed!(duals:: UpperTriangular, x:: UpperTriangular, seeds::NTuple)? And the same for Diagonal, Bidiagonal? For these it would be ideal not to generate all N^2 duals.

Note aside that Symmetric didn't work on 0.10, and is trickier:

julia> ForwardDiff.gradient(sum, Symmetric(randn(3, 3)))
ERROR: ArgumentError: Cannot set a non-diagonal index in a symmetric matrix
...
 [11] seed!(duals::Symmetric{…}, x::Symmetric{…}, seeds::NTuple{…})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/UBbGT/src/apiutils.jl:52

julia> ForwardDiff.gradient(x -> sum(Symmetric(x)), rand(3,3))
3×3 Matrix{Float64}:
 1.0  2.0  2.0
 0.0  1.0  2.0
 0.0  0.0  1.0

julia> Zygote.gradient(sum, Symmetric(rand(3,3)))[1]  # this has an explicit projection step
3×3 Symmetric{Float64, Matrix{Float64}}:
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0

julia> Zygote.gradient(x -> sum(Symmetric(x)), rand(3,3))[1]  # much the same as ForwardDiff
3×3 UpperTriangular{Float64, Matrix{Float64}}:
 1.0  2.0  2.0
  ⋅   1.0  2.0
  ⋅    ⋅   1.0

[penelopeysm]

It's a bit of an accident that this ever worked.

Indeed, I also discovered that UnitUpperTriangular and UnitLowerTriangular never worked (even on 0.10) due to similar reasons - just like Symmetric. (I totally appreciate the new equality behaviour is more principled too - it's just that for better or worse we have some code that implicitly relies on that old behaviour.)

In our case we can't do the intermediate UpperTriangular call since we get the triangular matrix directly from rand(LKJCholesky(3, 1.0, :U)). I agree, it seems a custom seed! might be the best thing here. I could take a stab at it...

[mcabbott]

Maybe worth appending that jacobian obeys the same projection, and works the same way on 1.0:

julia> x = UpperTriangular(reshape(1:9,3,3) ./ 1)
3×3 UpperTriangular{Float64, Matrix{Float64}}:
 1.0  4.0  7.0
  ⋅   5.0  8.0
  ⋅    ⋅   9.0

julia> ForwardDiff.jacobian(x -> x[:,2].^3, x)
3×9 Matrix{Float64}:
 0.0  0.0  0.0  48.0   0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  75.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0   0.0  0.0  0.0  0.0  0.0

And derivative can return an UpperTriangular, again behaves the same:

julia> ForwardDiff.derivative(x -> UpperTriangular(fill(x,3,3)) .^ 2, 3.0)
3×3 UpperTriangular{Float64, Matrix{Float64}}:
 6.0  6.0  6.0
  ⋅   6.0  6.0
  ⋅    ⋅   6.0

Would have to dig to understand the code paths these take. These two don't seem to behave differently with chunked mode. I didn't try out pre-prepared GradientConfig etc. But it might be worth adding tests for all these things.