From fc8320727bf460772c6a37040776e7138ae11f1d Mon Sep 17 00:00:00 2001 From: Orjan Ameye Date: Sat, 25 Jul 2026 10:14:03 +0200 Subject: [PATCH 1/3] perf: solve symbolic systems with Bareiss elimination `Symbolics.symbolic_linear_solve` returns nested fractions whose denominators carry the coefficient determinant. For a trigonometric ansatz those denominators hold identities like cos(x)^2 + sin(x)^2 = 1 that SymbolicUtils 4 no longer collapses, so `rearrange!` grew expressions without bound on larger systems (a van der Pol oscillator with two harmonics exhausted 16 GB and never finished). Solve the system with fraction-free Bareiss elimination instead, falling back to `symbolic_linear_solve` when the system is not polynomial in the unknowns. --- Project.toml | 4 + src/DifferentialEquation.jl | 13 ++- src/HarmonicEquation.jl | 19 ++-- src/QuestBase.jl | 3 + src/Symbolics/linear_solve.jl | 182 ++++++++++++++++++++++++++++++++++ test/DifferentialEquations.jl | 28 ++++++ test/symbolics.jl | 49 ++++++++- 7 files changed, 281 insertions(+), 17 deletions(-) create mode 100644 src/Symbolics/linear_solve.jl diff --git a/Project.toml b/Project.toml index 32a0639..ddafdb8 100644 --- a/Project.toml +++ b/Project.toml @@ -5,18 +5,22 @@ version = "0.4.0" [deps] DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" +DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3" OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" [compat] DocStringExtensions = "0.9.4" +DynamicPolynomials = "0.6" SymbolicUtils = "4" Symbolics = "7" julia = "1.12" Random = "1.10" LinearAlgebra = "1.10" +MultivariatePolynomials = "0.5" Test = "1.10" OrderedCollections = "1.8, 2" Aqua = "0.8.11" diff --git a/src/DifferentialEquation.jl b/src/DifferentialEquation.jl index b09e955..df764fb 100644 --- a/src/DifferentialEquation.jl +++ b/src/DifferentialEquation.jl @@ -182,9 +182,16 @@ the left-hand sides. Uses symbolic linear solving to determine the right-hand si the equations in place. """ function rearrange!(eom::DifferentialEquation, new_lhs::Vector{Num}) - soln = Symbolics.symbolic_linear_solve( - get_equations(eom), new_lhs; simplify=false, check=true - ) + soln = try + fraction_free_linear_solve(get_equations(eom), new_lhs) + catch error + error isa BareissFailure || rethrow() + Num.( + Symbolics.symbolic_linear_solve( + get_equations(eom), new_lhs; simplify=false, check=true + ), + ) + end # Use original variable keys (not extracted from new_lhs, as v7's get_variables # treats derivatives as variables) eom.equations = OrderedDict(zip(collect(keys(eom.equations)), new_lhs .~ soln)) diff --git a/src/HarmonicEquation.jl b/src/HarmonicEquation.jl index dbebc9b..fce4503 100644 --- a/src/HarmonicEquation.jl +++ b/src/HarmonicEquation.jl @@ -171,15 +171,16 @@ end "Rearrange an equation system such that the field equations is equal to the vector specified in new_lhs" function rearrange!(eom::HarmonicEquation, new_rhs::Vector{Num}) - soln = Symbolics.symbolic_linear_solve( - eom.equations, new_rhs; simplify=false, check=true - ) - # `symbolic_linear_solve` returns nested fractions whose denominators hold the - # coefficient determinant; for a trigonometric ansatz it contains identities like - # cos(ωt)² + sin(ωt)² = 1 that SymbolicUtils 4's `simplify` no longer collapses. - # Flatten the nest (polynomial-level, no AC matching) and reduce the denominator, - # so no consumer works with superficially time-dependent denominators. - soln = reduce_denominator.(Symbolics.simplify_fractions.(Num.(soln))) + soln = try + fraction_free_linear_solve(eom.equations, new_rhs) + catch error + error isa BareissFailure || rethrow() + fallback = Symbolics.symbolic_linear_solve( + eom.equations, new_rhs; simplify=false, check=true + ) + Symbolics.simplify_fractions.(Num.(fallback)) + end + soln = reduce_denominator.(soln) eom.equations = soln .~ new_rhs return nothing end diff --git a/src/QuestBase.jl b/src/QuestBase.jl index db05c1b..865caf0 100644 --- a/src/QuestBase.jl +++ b/src/QuestBase.jl @@ -1,6 +1,8 @@ module QuestBase using DocStringExtensions +import DynamicPolynomials as DP +import MultivariatePolynomials as MP using OrderedCollections: OrderedCollections, OrderedDict, OrderedSet using LinearAlgebra: LinearAlgebra @@ -40,6 +42,7 @@ include("Symbolics/Symbolics_utils.jl") include("Symbolics/exponentials.jl") include("Symbolics/fourier.jl") include("Symbolics/drop_powers.jl") +include("Symbolics/linear_solve.jl") include("DifferentialEquation.jl") include("Variables.jl") include("HarmonicVariable.jl") diff --git a/src/Symbolics/linear_solve.jl b/src/Symbolics/linear_solve.jl new file mode 100644 index 0000000..3137348 --- /dev/null +++ b/src/Symbolics/linear_solve.jl @@ -0,0 +1,182 @@ +struct BareissFailure <: Exception + message::String +end + +const _BAREISS_VARIABLE = only(DP.@polyvar __questbase_bareiss_variable) + +function Base.showerror(io::IO, error::BareissFailure) + return print(io, error.message) +end + +""" + fraction_free_linear_solve(equations, variables) + +Solve a square symbolic linear system using pivoted Bareiss–Jordan elimination. + +The coefficients are converted together into SymbolicUtils' sparse polynomial +representation. Bareiss elimination then performs only exact polynomial divisions +and produces one determinant denominator per solution, instead of the nested +fractions produced by symbolic LU. +""" +function fraction_free_linear_solve(equations, variables) + coefficients, offsets, islinear = Symbolics.linear_expansion(equations, variables) + islinear || throw(ArgumentError("the equation system is not linear in the variables")) + + rows, columns = size(coefficients) + rows == columns == length(offsets) || + throw(DimensionMismatch("fraction-free linear solve requires a square system")) + + return _bareiss_jordan_solve(coefficients, -offsets) +end + +function _bareiss_jordan_solve(coefficients, right_hand_side) + n = size(coefficients, 1) + symbolic_type = eltype(coefficients) + poly_variable_type = typeof(_BAREISS_VARIABLE) + poly_to_symbolic = Dict{poly_variable_type,symbolic_type}() + symbolic_to_poly = Dict{symbolic_type,poly_variable_type}() + augmented = Matrix{Any}(undef, n, n + 1) + + for column in 1:n, row in 1:n + augmented[row, column] = _to_polynomial( + coefficients[row, column], symbolic_to_poly, poly_to_symbolic + ) + end + for row in 1:n + augmented[row, n + 1] = _to_polynomial( + right_hand_side[row], symbolic_to_poly, poly_to_symbolic + ) + end + + previous_pivot = 1 + for pivot_index in 1:n + pivot_row = _find_bareiss_pivot(augmented, pivot_index, n) + isnothing(pivot_row) && throw(LinearAlgebra.SingularException(pivot_index)) + if pivot_row != pivot_index + for column in axes(augmented, 2) + augmented[pivot_index, column], augmented[pivot_row, column] = + augmented[pivot_row, column], augmented[pivot_index, column] + end + end + + old = copy(augmented) + pivot = old[pivot_index, pivot_index] + for row in 1:n + row == pivot_index && continue + for column in axes(augmented, 2) + column == pivot_index && continue + numerator = + pivot * old[row, column] - + old[row, pivot_index] * old[pivot_index, column] + augmented[row, column] = if pivot_index == 1 + numerator + else + _exact_polynomial_division(numerator, previous_pivot) + end + end + augmented[row, pivot_index] = 0 + end + previous_pivot = pivot + end + + solution = Vector{Num}(undef, n) + for row in 1:n + numerator_poly = augmented[row, n + 1] + denominator_poly = augmented[row, row] + leading_coefficient = if denominator_poly isa Number + denominator_poly + else + MP.leading_coefficient(denominator_poly) + end + if leading_coefficient isa Real && leading_coefficient < 0 + numerator_poly = -numerator_poly + denominator_poly = -denominator_poly + end + numerator = _from_polynomial(numerator_poly, poly_to_symbolic) + denominator = _from_polynomial(denominator_poly, poly_to_symbolic) + solution[row] = numerator / denominator + end + return solution +end + +_to_polynomial(expression::Num, symbolic_to_poly, poly_to_symbolic) = _to_polynomial( + unwrap(expression), symbolic_to_poly, poly_to_symbolic +) +_to_polynomial(expression::Number, _, _) = expression +function _to_polynomial(expression::BasicSymbolic, symbolic_to_poly, poly_to_symbolic) + if SymbolicUtils.isconst(expression) + return unwrap_const(expression) + elseif isadd(expression) + result = 0 + for argument in arguments(expression) + result += _to_polynomial(argument, symbolic_to_poly, poly_to_symbolic) + end + return result + elseif ismul(expression) + result = 1 + for argument in arguments(expression) + result *= _to_polynomial(argument, symbolic_to_poly, poly_to_symbolic) + end + return result + elseif ispow(expression) + base, exponent = arguments(expression) + if SymbolicUtils.isconst(exponent) + value = unwrap_const(exponent) + if value isa Integer && value >= 0 + return _to_polynomial(base, symbolic_to_poly, poly_to_symbolic)^value + end + end + end + + variable = get!(symbolic_to_poly, expression) do + name = Symbol(:questbase_atom_, hash(expression)) + MP.similar_variable(_BAREISS_VARIABLE, name) + end + get!(poly_to_symbolic, variable, expression) + return variable +end + +_from_polynomial(value::Number, _) = Num(value) +function _from_polynomial(variable::typeof(_BAREISS_VARIABLE), poly_to_symbolic) + return Num(poly_to_symbolic[variable]) +end +function _from_polynomial(polynomial, poly_to_symbolic) + variables = MP.variables(polynomial) + result = Num(0) + for term in MP.terms(polynomial) + expression = Num(MP.coefficient(term)) + exponents = MP.exponents(MP.monomial(term)) + for (variable, exponent) in zip(variables, exponents) + iszero(exponent) && continue + expression *= Num(poly_to_symbolic[variable])^exponent + end + result += expression + end + return result +end + +function _find_bareiss_pivot(augmented, pivot_index, n) + best_row = nothing + best_size = typemax(Int) + for row in pivot_index:n + candidate = augmented[row, pivot_index] + iszero(candidate) && continue + candidate_size = candidate isa Number ? 1 : length(MP.terms(candidate)) + if candidate_size < best_size + best_row = row + best_size = candidate_size + end + end + return best_row +end + +function _exact_polynomial_division(numerator, denominator) + quotient = try + MP.div_multiple(numerator, denominator) + catch error + throw(BareissFailure("Bareiss exact division failed: $(sprint(showerror, error))")) + end + iszero(quotient * denominator - numerator) || + throw(BareissFailure("Bareiss division produced a nonzero remainder")) + return quotient +end diff --git a/test/DifferentialEquations.jl b/test/DifferentialEquations.jl index cb9955d..2a2a39d 100644 --- a/test/DifferentialEquations.jl +++ b/test/DifferentialEquations.jl @@ -87,6 +87,34 @@ end @eqtest new_eq != diff_eq end +@testset "State-dependent derivative mass matrix" begin + @variables α β μ γx γy kx ky Fx Fy + ẍ = d(x, t, 2) + ÿ = d(y, t, 2) + ẋ = d(x, t) + ẏ = d(y, t) + + m11 = 1 + α * x^2 + m12 = μ * x + m21 = μ * y + m22 = 1 + β * y^2 + force_x = Fx * cos(ω * t) - γx * ẋ - kx * x + force_y = Fy * sin(ω * t) - γy * ẏ - ky * y + equations = [m11 * ẍ + m12 * ÿ ~ force_x, m21 * ẍ + m22 * ÿ ~ force_y] + system = DifferentialEquation(equations, [x, y]) + + rearrange_standard!(system) + rearranged = collect(values(system.equations)) + determinant = m11 * m22 - m12 * m21 + + @eqtest rearranged[1].lhs == ẍ + @eqtest rearranged[2].lhs == ÿ + @eqtest rearranged[1].rhs == + Symbolics.expand(m22 * force_x - m12 * force_y) / Symbolics.expand(determinant) + @eqtest rearranged[2].rhs == + Symbolics.expand(m11 * force_y - m21 * force_x) / Symbolics.expand(determinant) +end + @testset "Error Cases" begin @testset "Equation{Vector}" begin # define equation of motion diff --git a/test/symbolics.jl b/test/symbolics.jl index 95b7c05..f9c08f1 100644 --- a/test/symbolics.jl +++ b/test/symbolics.jl @@ -173,12 +173,51 @@ end @eqtest reduce_denominator(a / (ω^2 + 1)) == a / (ω^2 + 1) end +@testset "fraction-free linear solve" begin + using LinearAlgebra + using QuestBase: fraction_free_linear_solve + + @variables a b c x y z + equations = [2x + y - z ~ a, x + 3y + z ~ b, 2x - y + 4z ~ c] + solution = fraction_free_linear_solve(equations, [x, y, z]) + @eqtest solution[1] == (13a - 3b + 4c) / 31 + @eqtest solution[2] == (-2a + 10b - 3c) / 31 + @eqtest solution[3] == (-7a + 4b + 5c) / 31 + + # A zero diagonal requires a row pivot. + pivoted = fraction_free_linear_solve([y ~ a, x + y ~ b], [x, y]) + @eqtest pivoted[1] == b - a + @eqtest pivoted[2] == a + + @test_throws LinearAlgebra.SingularException fraction_free_linear_solve( + [x + y ~ a, 2x + 2y ~ b], [x, y] + ) + + # Compare arbitrary nonsingular integer systems against exact rational arithmetic. + rng = Random.MersenneTwister(0x5eed) + for dimension in 2:5 + matrix = rand(rng, -4:4, dimension, dimension) + while iszero(det(matrix)) + matrix = rand(rng, -4:4, dimension, dimension) + end + rhs = rand(rng, -4:4, dimension) + symbolic_variables = only(@variables q[1:dimension]) + random_equations = [ + sum(matrix[row, column] * symbolic_variables[column] for column in 1:dimension) ~ + rhs[row] for row in 1:dimension + ] + actual = fraction_free_linear_solve(random_equations, symbolic_variables) + expected = Rational{Int}.(matrix) \ Rational{Int}.(rhs) + for index in eachindex(actual) + @eqtest actual[index] == expected[index] + end + end +end + @testset "rearrange! reduces determinant denominators" begin - # `symbolic_linear_solve` returns nested fractions whose denominators hold the - # coefficient determinant; for a trigonometric ansatz it contains cos² + sin² = 1. - # `rearrange!` must flatten and collapse it, otherwise every consumer works with - # superficially time-dependent denominators (this stalled the Krylov-Bogoliubov - # order-2 averaging for hours). + # Fraction-free elimination produces one determinant denominator instead of nested + # LU fractions. For a trigonometric ansatz that determinant contains + # cos² + sin² = 1, which `rearrange!` must collapse before downstream averaging. using QuestBase: HarmonicEquation, HarmonicVariable, DifferentialEquation using QuestBase: rearrange!, d @variables ω t T u1(T) v1(T) x(t) From 332428a80f69840c1aba302bfe2e9049b961e3b6 Mon Sep 17 00:00:00 2001 From: Orjan Ameye Date: Sat, 25 Jul 2026 10:14:03 +0200 Subject: [PATCH 2/3] release: QuestBase v0.4.1 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index ddafdb8..4735438 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "QuestBase" uuid = "7e80f742-43d6-403d-a9ea-981410111d43" authors = ["Orjan Ameye ", "Jan Kosata ", "Javier del Pino "] -version = "0.4.0" +version = "0.4.1" [deps] DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" From f161b8e8c7fadabc3a83f9c9961b410f532084a5 Mon Sep 17 00:00:00 2001 From: Orjan Ameye Date: Sat, 25 Jul 2026 10:22:19 +0200 Subject: [PATCH 3/3] style: apply JuliaFormatter to the Bareiss solver --- src/Symbolics/linear_solve.jl | 13 +++++++------ test/symbolics.jl | 5 +++-- 2 files changed, 10 insertions(+), 8 deletions(-) diff --git a/src/Symbolics/linear_solve.jl b/src/Symbolics/linear_solve.jl index 3137348..4de9a4c 100644 --- a/src/Symbolics/linear_solve.jl +++ b/src/Symbolics/linear_solve.jl @@ -54,8 +54,9 @@ function _bareiss_jordan_solve(coefficients, right_hand_side) isnothing(pivot_row) && throw(LinearAlgebra.SingularException(pivot_index)) if pivot_row != pivot_index for column in axes(augmented, 2) - augmented[pivot_index, column], augmented[pivot_row, column] = - augmented[pivot_row, column], augmented[pivot_index, column] + swapped = augmented[pivot_index, column] + augmented[pivot_index, column] = augmented[pivot_row, column] + augmented[pivot_row, column] = swapped end end @@ -99,9 +100,9 @@ function _bareiss_jordan_solve(coefficients, right_hand_side) return solution end -_to_polynomial(expression::Num, symbolic_to_poly, poly_to_symbolic) = _to_polynomial( - unwrap(expression), symbolic_to_poly, poly_to_symbolic -) +function _to_polynomial(expression::Num, symbolic_to_poly, poly_to_symbolic) + return _to_polynomial(unwrap(expression), symbolic_to_poly, poly_to_symbolic) +end _to_polynomial(expression::Number, _, _) = expression function _to_polynomial(expression::BasicSymbolic, symbolic_to_poly, poly_to_symbolic) if SymbolicUtils.isconst(expression) @@ -130,7 +131,7 @@ function _to_polynomial(expression::BasicSymbolic, symbolic_to_poly, poly_to_sym variable = get!(symbolic_to_poly, expression) do name = Symbol(:questbase_atom_, hash(expression)) - MP.similar_variable(_BAREISS_VARIABLE, name) + return MP.similar_variable(_BAREISS_VARIABLE, name) end get!(poly_to_symbolic, variable, expression) return variable diff --git a/test/symbolics.jl b/test/symbolics.jl index f9c08f1..b02b3f5 100644 --- a/test/symbolics.jl +++ b/test/symbolics.jl @@ -203,8 +203,9 @@ end rhs = rand(rng, -4:4, dimension) symbolic_variables = only(@variables q[1:dimension]) random_equations = [ - sum(matrix[row, column] * symbolic_variables[column] for column in 1:dimension) ~ - rhs[row] for row in 1:dimension + sum( + matrix[row, column] * symbolic_variables[column] for column in 1:dimension + ) ~ rhs[row] for row in 1:dimension ] actual = fraction_free_linear_solve(random_equations, symbolic_variables) expected = Rational{Int}.(matrix) \ Rational{Int}.(rhs)