fix crash of eigen with #undef elements

[araujoms]

The functions _structured_value and _structured_partials access the parent of the matrix, which causes a crash when it contains #undef elements, as is the case if we only partially initialize a BigFloat matrix. Here I'm switching map to broadcasting, which can deal with it correctly.

For testing I assume you don't want GenericLinearAlgebra as a test dependency, so I wrote some mock eigvals and eigen functions instead.

The MWE to get the crash is to run test_gradient(BigFloat) below:

using LinearAlgebra, GenericLinearAlgebra

function test_gradient(::Type{T}) where {T}
    d = 2
    
    function barrier(point)
        M = similar(point, d, d)
        counter = 0
        for j in 1:d, i in 1:j
            counter += 1
            M[i, j] = point[counter]
        end
        return real(tr(log(Hermitian(M))))
    end

    p0 = randn(T, div(d*(d+1), 2))

    return ForwardDiff.gradient(barrier, p0)
end

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 91.33%. Comparing base (a3c0f4f) to head (3520531).

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #851      +/-   ##
==========================================
+ Coverage   91.23%   91.33%   +0.10%     
==========================================
  Files          11       11              
  Lines        1072     1073       +1     
==========================================
+ Hits          978      980       +2     
+ Misses         94       93       -1     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[devmotion]

It seems indexing A instead of parent(A) is noticeably slower: it adds a branch and a strided read per element, and on top of that broadcast over a Matrix seems ~1.4x slower than map over it. For _structured_value of Symmetric/Hermitian that's a factor of ~2.2, and for _structured_partials of Hermitian{Complex} ~1.66.

I think we can avoid this regression by specializing on isbitstype(V). Something like

_maptri(f, A::Union{Symmetric,Hermitian}) = _maptri(f, A, parent(A))
_maptri(f, A, P::AbstractArray{V}) where {V} = isbitstype(V) ? map(f, P) : broadcast(f, A)

_structured_value(A::Symmetric{Dual{T,V,N}}) where {T,V,N} = Symmetric(_maptri(value, A), A.uplo === 'U' ? :U : :L)
_structured_value(A::Hermitian{Dual{T,V,N}}) where {T,V,N} = Hermitian(_maptri(value, A), A.uplo === 'U' ? :U : :L)
_structured_value(A::Hermitian{Complex{Dual{T,V,N}}}) where {T,V,N} = Hermitian(_maptri(z -> complex(value(real(z)), value(imag(z))), A), A.uplo === 'U' ? :U : :L)

_structured_partials(A::Symmetric{Dual{T,V,N}}, j::Int) where {T,V,N} = Symmetric(_maptri(a -> partials(a, j), A), A.uplo === 'U' ? :U : :L)
_structured_partials(A::Hermitian{Dual{T,V,N}}, j::Int) where {T,V,N} = Hermitian(_maptri(a -> partials(a, j), A), A.uplo === 'U' ? :U : :L)
function _structured_partials(A::Hermitian{Complex{Dual{T,V,N}}}, j::Int) where {T,V,N}
    return Hermitian(_maptri(z -> complex(partials(real(z), j), partials(imag(z), j)), A), A.uplo === 'U' ? :U : :L)
end

My impression was also that the two complex methods have to be fused into a single closure, otherwise A[i, j] is evaluated twice per element.

Regarding the tests: eigvals(::Hermitian{Complex{BigFloat}}) returns complex eigenvalues. This only works because Calculus.finite_difference_jacobian writes into a Matrix{Float64} and the imaginary parts happen to be exactly zero - with a rounding-level imaginary part the test errors with InexactError instead of failing. I suspect this is also what #workaround for a bug in julia 1.10 is about: the error it avoids is a MethodError in _to_duals caused by these complex eigenvalues, and it occurs on 1.12 as well. If eigvals returns real values the analytic eigen works on both versions, and then Hermitian(real(Matrix(M))) can be dropped as well - it discards imag(M), so currently the imaginary parts in _structured_value/_structured_partials are only exercised with zeros.

It also seems the #undef tests cover only uplo = :U - check uplo = :L as well?

[araujoms]

The julia-1.10 bug was that if you tried to take real(M) it would also crash on undefined elements. The workaround was to materialise the matrix. Funnily enough complex(M) works, so I switched to defining only the functions for the complex case, and everything else in terms of it.

I also added the specialisation for isbitstype, tests for both uplos, and genuinely complex matrices.

[devmotion]

Thanks, the performance regression, the fused closures, the complex eigenvalues and the missing uplo = :L tests are all fixed now.

A few remaining points:

• _maptri could be a single method: for Symmetric/Hermitian, eltype(A) == eltype(parent(A)), so

  # non-isbits storage may be #undef outside of the `uplo` triangle
  _maptri(f, A::Union{Symmetric,Hermitian}) = isbitstype(eltype(A)) ? map(f, parent(A)) : broadcast(f, A)

• The comment "the M[2, 1] element is undef" is only correct for uplo = :U.
• The mocks are type piracy of LinearAlgebra.eigvals/eigen, and they leak into all test files included after JacobianTest.jl. The closed-form 2x2 eigendecomposition is also quite a lot of code to review. I think it would be better to add GenericLinearAlgebra as a test dependency and remove the mocks. The tests would then also cover the code path from the MWE.
• Could we use an imaginary part that is independent of the real part, e.g. x[2] + im*x[1] as in the uplo testset above?
• It seems the new tests could be merged into the uplo testset above, or at least use named testsets (@testset "uplo = :$uplo" for ..., @testset "$name" for (name, wrap) in ...) as it does.
• The same issue exists for SymTridiagonal with length(ev) == length(dv) and an #undef last element of ev, e.g. ev = similar(x, 2); ev[1] = x[3]; eigvals(SymTridiagonal(x[1:2], ev)) with BigFloat inputs. That could be fixed separately, though.

[araujoms]

Importing GenericLinearAlgebra would make the piracy problem much worse, as it pirates not only eigen and eigvals but also several others. Here at least I'm only pirating it for BigFloat elements. If you insist I'll do it, though.

I don't think there's a point fixing it for SymTridiagonal, as support for length(ev) == length(dv) is an undocumented internal that has been removed in Julia 1.14: JuliaLang/LinearAlgebra.jl#1569

Everything else I did as requested.

[devmotion]

• The PR deletes the existing uplo = :$uplo testset. That's the only test that checks values and partials are read from the uplo triangle when the storage isn't symmetric, so it should be restored.
• I'd still prefer to drop the mocks. They only work for 2x2: for any other size eigvals silently returns wrong values (e.g. M[4] is entry (1, 2) of a 3x3 matrix), and eigen divides by M[1, 2]. They also stay defined for all test files included after JacobianTest.jl. The fix can be tested directly without any eigensolver, e.g. inside the restored uplo testset (passes with this PR, UndefRefError on master):

  @testset "#undef outside of the `uplo` triangle" begin
      x = ForwardDiff.Dual{Nothing}.(BigFloat[1, 2, 3], BigFloat[1, 0, 0], BigFloat[0, 1, 0])
      k = uplo === :U ? 3 : 2 # linear index of the stored off-diagonal element
      M = similar(x, 2, 2)
      M[1, 1], M[k], M[2, 2] = x[1], x[2], x[3]
      Mc = similar(x, Complex{eltype(x)}, 2, 2)
      Mc[1, 1], Mc[k], Mc[2, 2] = x[1], x[2] + im * x[1], x[3]
      s = uplo === :U ? 1 : -1 # sign of imag(A[1, 2])
      @testset "$name" for (name, A, value, ∂1, ∂2) in (
          ("Symmetric{<:Real}", Symmetric(M, uplo), [1 2; 2 3], [1 0; 0 0], [0 1; 1 0]),
          ("Hermitian{<:Real}", Hermitian(M, uplo), [1 2; 2 3], [1 0; 0 0], [0 1; 1 0]),
          ("Hermitian{<:Complex}", Hermitian(Mc, uplo), [1 2+s*im; 2-s*im 3], [1 s*im; -s*im 0], [0 1; 1 0]),
      )
          @test ForwardDiff._structured_value(A) == value
          @test ForwardDiff._structured_partials(A, 1) == ∂1
          @test ForwardDiff._structured_partials(A, 2) == ∂2
      end
  end

• The Float64 inputs in the new testset don't produce #undef elements, only uninitialized memory. So these tests run on garbage values that happen to be ignored.
• length(ev) == length(dv) is removed in 1.14, but ForwardDiff still supports 1.10–1.13, where it's accepted. The fix is simple (only map over the first length(dv) - 1 elements of ev), so I think it would be good to include it here.

[araujoms]

This is Kafkaesque. You asked me to merge the uplo test with the undef test, I did so, and then you complain that I deleted the uplo test.

Fine, I deleted everything, added the test you wanted, and now fixed SymTridiagonal as well.