Inlined linear SCC of three unknowns: the run-time solve dominates the cost of the right-hand side

Aperta
  • Qualcosa non va?
#173 3 commenti 0 reazioni 0 assegnatari Vedi su GitHub

Nessuno ha ancora preso questa issue.

Valutazione

Difficoltà
5/5
Tempo stimato
Più di una settimana
Idoneità per principianti
35/100
Tipo di issue
Funzionalità
Chiarezza
Abbastanza chiara
Stato di attività
Attiva
Stack tecnologico
julia
Ambito
performance

Direzione di ricerca

Start by reading the issue’s benchmark and the inline_linear_sccs and analytical_linear_scc_limit options in ModelingToolkit’s generated right-hand side path. Compare the current inlined solve with the StaticArrays example and check how DyadCompilerPasses handles solve buffers. Done means determining whether a supported option can generate a static solve for small SCCs, including the singular-matrix tradeoff, and where that option should live.

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

Descrizione

Goal. A fast generated right-hand side for models with small linear algebraic loops. These are thermal networks with nodes without heat capacity, where the loop matrix depends on parameters and inputs. Such a model is evaluated millions of times in a Kalman filter during parameter estimation, so the cost of f determines the cost of the fit.

What we observe. With inline_linear_sccs = true, a loop of three unknowns is solved at run time through a general linear solve with cached work buffers. This costs about 26 times as much as solving the same 3 × 3 system with StaticArrays:

using ModelingToolkit, LinearAlgebra, StaticArrays, BenchmarkTools
using ModelingToolkit: t_nounits as t, D_nounits as D

@variables x(t) = 20.0 T1(t) T2(t) T3(t)
@parameters g1 = 2.0 g2 = 3.0 g3 = 1.5 g12 = 0.7 g23 = 0.4 g13 = 0.2 Ta = 0.0 C = 10.0
eqs = [
    C * D(x) ~ g1 * (T1 - x) + g2 * (T2 - x),
    0 ~ g1 * (x - T1) + g12 * (T2 - T1) + g13 * (T3 - T1),
    0 ~ g2 * (x - T2) + g12 * (T1 - T2) + g23 * (T3 - T2),
    0 ~ g3 * (Ta - T3) + g13 * (T1 - T3) + g23 * (T2 - T3),
]
@named sys = System(eqs, t)
ssys = mtkcompile(sys; reassemble_alg = ModelingToolkit.StructuralTransformations.DefaultReassembleAlgorithm(;
    inline_linear_sccs = true, analytical_linear_scc_limit = 1))
prob = ODEProblem(ssys, [], (0.0, 1.0))

function static_rhs!(du, u, g)  # the same right-hand side, the loop as a static 3 × 3 system
    g1, g2, g3, g12, g23, g13, Ta, C = g
    A = @SMatrix [-(g1 + g12 + g13) g12 g13; g12 -(g2 + g12 + g23) g23; g13 g23 -(g3 + g13 + g23)]
    T = A \ (-@SVector [g1 * u[1], g2 * u[1], g3 * Ta])
    du[1] = (g1 * (T[1] - u[1]) + g2 * (T[2] - u[1])) / C
    du
end
du = zeros(1); u = [20.0]; g = (2.0, 3.0, 1.5, 0.7, 0.4, 0.2, 0.0, 10.0)
@btime $(prob.f)($du, $u, $(prob.p), 0.0)   # 164.7 ns, 3 allocations (128 bytes)
@btime static_rhs!($du, $u, $g)             # 6.3 ns, 0 allocations

(ModelingToolkit 11.45.5, ModelingToolkitBase 1.77.5, ModelingToolkitTearing 1.20.7, StateSelection 1.11.1, LinearSolve 5.18.2, Julia 1.13.1.)

In our model of a house (state dimension 9, one inlined 3 × 3 loop) the solve and its buffers are a large part of the cost of f. Under ForwardDiff with nested dual numbers they also cause most of the allocations of the Jacobian, about 35–75 GB per Jacobian of the likelihood.

What we tried.

  • analytical_linear_scc_limit = 3: no effect. The numeric solve is kept, presumably because the loop matrix contains parameters or inputs that may be zero.
  • The optimize levels of DyadCompilerPasses (:basic, :aggressive) at code generation: no measurable change. The ldiv rule replaces the solve by another cached solve. The static-array rule is not reached, because the system is filled into preallocated buffers rather than written as a literal array.
  • Solving the inlined system with StaticArrays, by overwriting the solve method in our code: about 1.8 times faster residuals in the full model, identical results (6e-15). We do not want to keep a method overwrite.

Question. Is there an intended option to generate a static (or fully unrolled) solve for small inlined linear SCCs, or to let the analytical solution apply when the matrix may be singular for some parameter values? If not, would such an option belong here or in DyadCompilerPasses?

Host: battomenstation · session https://claude.ai/code/session_01LxnzqRyF2ZMkgWPTEPb44g · working directory ~/.julia/dev/HomeAssistantBot.jl/.claude/worktrees/precious-inventing-map

Lingua principale
Julia
Stelle
8
Fork
14
Merge medio
8g 8h
PR unite (30g)
1

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 JuliaComputing/StateSelection.jl

Tutte le issue di JuliaComputing/StateSelection.jl

Issue simili

Altre issue su Julia

Ricevi le nuove issue nella tua casella

Un breve riepilogo di issue GitHub adatte ai principianti.