Add analyze_initialization_jacobian — diagnose rank-deficient initialization

lib/ModelingToolkitBase/src/problems/initializationproblem.jl

@@ -131,6 +131,9 @@ function overdetermined_initialization_message(neqs::Integer, nunknown::Integer,
     Initialization system is overdetermined. $neqs equations for $nunknown unknowns. \
     Initialization will default to using least squares. $(extra)
 
+    Call `analyze_initialization_jacobian(prob)` on the constructed problem to see which \
+    equations are redundant (and which unknowns, if any, remain underdetermined).
+
     To suppress this warning, pass `warn_initialize_determined = false`. To turn this \
     warning into an error, pass `fully_determined = true`.
     """
@@ -141,6 +144,9 @@ function underdetermined_initialization_message(neqs::Integer, nunknown::Integer
     Initialization system is underdetermined. $neqs equations for $nunknown unknowns. \
     Initialization will default to using least squares. $(extra)
 
+    Call `analyze_initialization_jacobian(prob)` on the constructed problem to see which \
+    unknowns are underdetermined (and which equations, if any, are redundant).
+
     To suppress this warning, pass `warn_initialize_determined = false`. To turn this \
     warning into an error, pass `fully_determined = true`.
     """

src/ModelingToolkit.jl

@@ -152,6 +152,7 @@ end
 @reexport using .StructuralTransformations
 
 export SemilinearODEFunction, SemilinearODEProblem
+export analyze_initialization_jacobian
 export alias_elimination
 export linearize, linearization_function,
     LinearizationProblem, linearization_ap_transform

src/initialization.jl

@@ -47,3 +47,165 @@ function MTKBase.get_initialization_problem_type(
         NonlinearLeastSquaresProblem
     end
 end
+
+"""
+    analyze_initialization_jacobian(prob; rtol = 1e-8, atol = 0.0, threshold = 1e-3, verbose = true)
+
+Diagnose rank deficiency of a system's initialization problem by inspecting the singular
+value decomposition of its residual Jacobian, and report both the **unknowns** that span
+its (right) null space and the **equations** that are locally redundant (its left null
+space).
+
+When the initialization system is rank deficient the solve can converge to different
+solutions, or fail, depending on the solver and on the (often nondeterministic) ordering of
+the assembled unknowns/equations — a common source of intermittent initialization failures:
+some orderings land on a singular realization and throw, others succeed. This utility makes
+the offending degrees of freedom explicit:
+
+  - **Underdetermined unknowns** (right null space): directions along which the
+    initialization is free to move. Pinning these unknowns (or adding equations that
+    constrain them) makes the initialization well-posed and deterministic.
+  - **Redundant equations** (left null space): combinations of initialization equations
+    whose Jacobian rows are linearly dependent at `u0`, i.e. equations that do not locally
+    constrain any additional degree of freedom. These explain why an initialization can
+    have more equations than unknowns yet still be underdetermined.
+
+It evaluates the Jacobian of the initialization residual at the initial guess `u0` (via
+`ForwardDiff`) and computes its SVD. `prob` may be a problem that carries initialization
+data (e.g. an `ODEProblem`/`DAEProblem` built from a `System`), or an initialization
+`NonlinearProblem`/`NonlinearLeastSquaresProblem` directly.
+
+# Keyword arguments
+  - `rtol`: a singular value `σ` is treated as zero when `σ ≤ max(atol, rtol * σmax)`,
+    where `σmax` is the largest singular value. Increase it to also surface near-singular
+    directions.
+  - `atol`: absolute singular-value threshold (see `rtol`).
+  - `threshold`: only unknowns/equations whose participation exceeds this value are
+    reported.
+  - `verbose`: print a human-readable report.
+
+# Returns
+
+A `NamedTuple` `(; jacobian, singular_values, rank, nullity, redundancy, underdetermined_unknowns, redundant_equations)`.
+
+  - `nullity = ncols - rank` is the number of underdetermined directions and `redundancy =
+    nrows - rank` the number of redundant equations.
+  - `underdetermined_unknowns` is a vector of `unknown => weight` pairs and
+    `redundant_equations` a vector of `equation => weight` pairs, each sorted by decreasing
+    `weight ∈ [0, 1]`. The weight is the diagonal of the corresponding null-space projector
+    (the squared row norm over an orthonormal basis of the right/left null space): how
+    strongly that unknown/equation participates in the underdetermined/redundant directions,
+    independent of the arbitrary basis chosen within the null space.
+
+A `nullity` of `0` means the Jacobian has full column rank at `u0` (no underdetermined
+unknowns); a `redundancy` of `0` means full row rank (no redundant equations).
+
+!!! note
+    The Jacobian is evaluated at a single point (`u0`), so this reports the *local* rank
+    structure there. A structurally well-posed initialization can still be numerically
+    rank deficient at a particular operating point, and vice versa.
+"""
+function analyze_initialization_jacobian(
+        prob; rtol = 1e-8, atol = 0.0, threshold = 1e-3, verbose = true)
+    empty_result = (; jacobian = nothing, singular_values = Float64[], rank = 0,
+        nullity = 0, redundancy = 0, underdetermined_unknowns = Pair[],
+        redundant_equations = Pair[])
+    iprob = _initialization_problem(prob)
+    if iprob === nothing
+        verbose &&
+            @info "No initialization problem to analyze: the system is fully determined by its initial conditions."
+        return empty_result
+    end
+    u0 = state_values(iprob)
+    if u0 === nothing || isempty(u0)
+        verbose && @info "The initialization problem has no unknowns to solve for."
+        return empty_result
+    end
+    u0 = collect(float.(u0))
+    p = parameter_values(iprob)
+    f = iprob.f
+    nresid = f.resid_prototype !== nothing ? length(f.resid_prototype) : length(u0)
+    residual = if SciMLBase.isinplace(f)
+        u -> (du = similar(u, nresid); f(du, u, p); du)
+    else
+        u -> f(u, p)
+    end
+    J = ForwardDiff.jacobian(residual, u0)
+    nrows, ncols = size(J)
+    fact = svd(J; full = true)
+    S = fact.S
+    σmax = isempty(S) ? zero(eltype(S)) : maximum(S)
+    tol = max(atol, rtol * σmax)
+    rank = count(>(tol), S)
+    # A singular vector is a null direction when its singular value is below the
+    # tolerance. Right vectors (columns of `V`) beyond `length(S)` (present when
+    # ncols > nrows) and left vectors (columns of `U`) beyond `length(S)` (present when
+    # nrows > ncols) have an implied singular value of zero and are always null.
+    col_isnull = [(j <= length(S) ? S[j] : zero(σmax)) <= tol for j in 1:ncols]
+    row_isnull = [(i <= length(S) ? S[i] : zero(σmax)) <= tol for i in 1:nrows]
+    nullity = count(col_isnull)
+    redundancy = count(row_isnull)
+    # Participation = diagonal of the null-space projector (squared row norm over an
+    # orthonormal null-space basis), invariant to the arbitrary basis within the null space.
+    col_w = nullity == 0 ? zeros(ncols) :
+            vec(sum(abs2, @view(fact.V[:, col_isnull]); dims = 2))
+    row_w = redundancy == 0 ? zeros(nrows) :
+            vec(sum(abs2, @view(fact.U[:, row_isnull]); dims = 2))
+    syms = variable_symbols(iprob)
+    eqs = _initialization_equations(iprob)
+    underdetermined_unknowns = sort(
+        [syms[i] => col_w[i] for i in 1:ncols if col_w[i] > threshold];
+        by = last, rev = true)
+    redundant_equations = sort(
+        [(eqs === nothing ? i : eqs[i]) => row_w[i] for i in 1:nrows if row_w[i] > threshold];
+        by = last, rev = true)
+    if verbose
+        println("Initialization Jacobian rank analysis")
+        println("  residual Jacobian: $(nrows)×$(ncols), rank ≈ $rank, nullity ≈ $nullity, redundancy ≈ $redundancy")
+        if !isempty(S)
+            k = min(5, length(S))
+            println("  smallest $k singular value(s): ",
+                round.(S[(end - k + 1):end], sigdigits = 3))
+        end
+        if nullity == 0
+            println("  Full column rank at u0 — no underdetermined unknowns.")
+        else
+            println("  Underdetermined unknowns (participation ∈ [0, 1]):")
+            for (s, w) in underdetermined_unknowns
+                println("    ", lpad(string(round(w, digits = 4)), 8), "  ", s)
+            end
+        end
+        if redundancy == 0
+            println("  Full row rank at u0 — no redundant equations.")
+        else
+            println("  Redundant equations (participation ∈ [0, 1]):")
+            for (e, w) in redundant_equations
+                println("    ", lpad(string(round(w, digits = 4)), 8), "  ", e)
+            end
+        end
+    end
+    return (; jacobian = J, singular_values = S, rank, nullity, redundancy,
+        underdetermined_unknowns, redundant_equations)
+end
+
+# Symbolic equations of an initialization problem, in residual order, or `nothing` if not
+# available (e.g. a problem not built from a `System`).
+function _initialization_equations(iprob)
+    f = iprob.f
+    hasproperty(f, :sys) || return nothing
+    sys = f.sys
+    sys === nothing && return nothing
+    return equations(sys)
+end
+
+# Return the initialization `NonlinearProblem`/`NonlinearLeastSquaresProblem` carried by
+# `prob`, or `prob` itself if it is already a nonlinear problem, or `nothing` if there is
+# no initialization problem to analyze.
+function _initialization_problem(prob)
+    prob isa SciMLBase.AbstractNonlinearProblem && return prob
+    f = prob.f
+    if hasproperty(f, :initialization_data) && f.initialization_data !== nothing
+        return f.initialization_data.initializeprob
+    end
+    return nothing
+end

test/group_initialization.jl

@@ -2,5 +2,6 @@ include("shared/mtktestset.jl")
 
 @mtktestset("Guess Propagation", "guess_propagation.jl")
 @safetestset "Hierarchical Initialization Equations" include("hierarchical_initialization_eqs.jl")
+@safetestset "Initialization Jacobian analysis" include("initialization_jacobian_analysis.jl")
 @mtktestset("InitializationSystem Test", "initializationsystem.jl")
 @mtktestset("Initial Values Test", "initial_values.jl")

test/initialization_jacobian_analysis.jl

@@ -0,0 +1,53 @@
+using ModelingToolkit, Test
+using ModelingToolkit: t_nounits as t, D_nounits as D
+
+# Underdetermined initialization: a point constrained to the unit circle, with no
+# initial condition pinning where on the circle (or with what velocity) it starts.
+# The initialization Jacobian is column-rank deficient, so the analysis must flag a
+# nonempty set of underdetermined unknowns.
+@variables x(t) y(t)
+@named sys = System([D(x) ~ y, 0 ~ x^2 + y^2 - 1], t)
+sys = mtkcompile(sys)
+prob = ODEProblem(sys, [], (0.0, 1.0); guesses = [x => 0.6, y => 0.8],
+    warn_initialize_determined = false)
+res = analyze_initialization_jacobian(prob; verbose = false)
+@test res.nullity ≥ 1
+@test !isempty(res.underdetermined_unknowns)
+@test res.jacobian !== nothing
+@test all(0 .≤ last.(res.underdetermined_unknowns) .≤ 1 + 1e-8)
+@test res.rank + res.nullity == size(res.jacobian, 2)
+
+# Overdetermined initialization with a redundant equation: two consistent initial
+# conditions on a single unknown. The analysis must report a redundant equation.
+@variables z(t)
+@named sys2 = System([D(z) ~ -z], t; initialization_eqs = [z ~ 1.0, z^2 ~ 1.0])
+sys2 = mtkcompile(sys2)
+prob2 = ODEProblem(sys2, [], (0.0, 1.0); guesses = [z => 0.9],
+    warn_initialize_determined = false)
+res2 = analyze_initialization_jacobian(prob2; verbose = false)
+@test res2.redundancy ≥ 1
+@test !isempty(res2.redundant_equations)
+@test all(0 .≤ last.(res2.redundant_equations) .≤ 1 + 1e-8)
+
+# Determined initialization: an algebraic unknown uniquely fixed by a pinned state has a
+# full-rank initialization Jacobian — no underdetermined unknowns and no redundant
+# equations.
+@variables u(t) w(t)
+@named sys3 = System([D(u) ~ -u, 0 ~ w - exp(u)], t)
+sys3 = mtkcompile(sys3)
+prob3 = ODEProblem(sys3, [u => 1.0], (0.0, 1.0); guesses = [w => 1.0])
+res3 = analyze_initialization_jacobian(prob3; verbose = false)
+@test res3.nullity == 0
+@test res3.redundancy == 0
+@test isempty(res3.underdetermined_unknowns)
+@test isempty(res3.redundant_equations)
+
+# A system fully determined by its initial conditions has no initialization problem to
+# analyze; the utility handles that gracefully rather than erroring.
+@variables a(t) b(t)
+@named sys4 = System([D(a) ~ -a, D(b) ~ -b], t)
+sys4 = mtkcompile(sys4)
+prob4 = ODEProblem(sys4, [a => 1.0, b => 2.0], (0.0, 1.0))
+res4 = analyze_initialization_jacobian(prob4; verbose = false)
+@test res4.nullity == 0
+@test res4.redundancy == 0