jacobian(f, x) silently returns zeros if f is complex-valued

Aperta
#830 0 commenti 0 reazioni 0 assegnatari Vedi su GitHub

I maintainer di solito rispondono entro 1 giorno

Nessuno ha ancora preso questa issue.

Valutazione

Difficoltà
3/5
Tempo stimato
1-2 giorni
Idoneità per principianti
75/100
Tipo di issue
Bug
Chiarezza
Specificata chiaramente
Stato di attività
Tranquilla
Stack tecnologico
julia
Ambito
tooling

Direzione di ricerca

Inizia leggendo extract_jacobian!, extract_jacobian_chunk!, extract_derivative e valtype, concentrandoti su come vengono gestiti i valori Complex{Dual}. Esegui il test della modalità vettoriale, della modalità chunk e di jacobian! con il riproduttore dell’issue. Il lavoro è completato quando le funzioni a valori complessi producono una matrice jacobiana ComplexF64 corretta, mentre i percorsi a valori reali rimangono invariati.

Scritto dal modello di indicizzazione a partire dal testo della issue.

Descrizione

jacobian returns zeros, with Duals leaked into the eltype, instead of the Jacobian or an error:

julia> g(x) = [cis(x[1])*x[2], x[1] + 2im*x[2]];

julia> J = ForwardDiff.jacobian(g, [0.5, 1.0]);

julia> eltype(J)
Complex{ForwardDiff.Dual{ForwardDiff.Tag{typeof(g), Float64}, Float64, 2}}

julia> all(iszero, J)
true

extract_jacobian! extracts with partials(T, ydual, i), which for a Complex{Dual} hits the partials(x, i...) = zero(x) fallback. derivative(f, x) is fine since it goes through extract_derivative, which handles Complex{<:Dual}.

extract_derivative cannot be reused directly because it hardcodes the partial index, while the Jacobian needs index i. But generalising it to an index is enough, and keeps the knowledge of how a complex number carries derivative information in one place:

@inline extract_partial(::Type{T}, y::Real, i) where {T} = zero(y)
@inline extract_partial(::Type{T}, y::Complex, i) where {T} = zero(y)
@inline extract_partial(::Type{T}, y::Dual, i) where {T} = partials(T, y, i)
@inline extract_partial(::Type{T}, y::Complex{<:Dual}, i) where {T} =
    complex(partials(T, real(y), i), partials(T, imag(y), i))

extract_jacobian!/extract_jacobian_chunk! then use extract_partial(T, ·, i) instead of partials(T, ·, i), and the scalar extract_derivative methods become the i = 1 case.

One more piece is needed for the result eltype, since jacobian allocates similar(ydual, valtype(T, eltype(ydual)), ...):

@inline valtype(::Type{T}, ::Type{Complex{D}}) where {T,D} = Complex{valtype(T, D)}

With those, vector mode, chunk mode and jacobian! all give the correct ComplexF64 Jacobian and the real paths are unchanged.

If we would rather not support complex-valued f here, it should at least throw rather than return zeros.

Lingua principale
Julia
Stelle
1k
Fork
165
Merge medio
8g 22h
PR unite (30g)
3

Preparare l'ambiente

Questo progetto non fornisce container di sviluppo, Dockerfile né guida per i contributori, quindi l'ambiente è a tuo carico: parti dal suo README e consulta la nostra guida al primo contributo per i passaggi generali.

Come iniziare

  1. Leggi tutta la issue e poi la guida ai contributi del progetto.
  2. Commenta sulla issue per dire che te ne occupi tu — evita che due persone facciano lo stesso lavoro.
  3. Fai un fork del repository e lavora su un branch.
  4. Apri una pull request che faccia riferimento al numero della issue.

Altre issue di JuliaDiff/ForwardDiff.jl

Tutte le issue di JuliaDiff/ForwardDiff.jl

Issue simili

Altre issue su Julia

Ricevi le nuove issue nella tua casella

Un breve riepilogo di issue GitHub adatte ai principianti.