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

オープン
  • 何か問題がありますか？
#173 コメント 3 件 リアクション 0 件 担当者 0 名 GitHub で見る

まだ誰も着手していません。

評価

難易度
5/5
見積もり時間
1週間以上
初心者へのやさしさ
35/100
issue の種類
機能追加
明瞭さ
おおむね明確
活発さ
活発
技術スタック
julia
領域
performance

調査の方向性

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.

索引モデルが issue の本文から書いたものです。

説明

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

主要言語
Julia
スター
8
フォーク
14
平均マージ
8日 8時間
マージ済み PR（30日）
1

環境構築

このプロジェクトには開発コンテナ、Dockerfile、コントリビューションガイドがありません。まず README を読み、一般的な手順ははじめてのコントリビューションガイドを参照してください。

はじめの一歩

  1. issue を最後まで読み、次にプロジェクトのコントリビューションガイドを読みます。
  2. 着手することを issue にコメントします — 二人が同じ作業をするのを防げます。
  3. リポジトリをフォークし、ブランチを切って変更します。
  4. issue 番号を参照したプルリクエストを送ります。

JuliaComputing/StateSelection.jl のほかの issue

JuliaComputing/StateSelection.jl の issue をすべて見る

似ている issue

Julia の issue をもっと見る

新しい issue をメールで受け取る

初心者向けの GitHub issue を短くまとめたダイジェスト。