From 639daa5ecd4c95e37981a78c488c781460335ad6 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Thu, 12 Dec 2024 23:53:00 +0100 Subject: [PATCH 1/5] Assume unknown functions are non-linear in `hessian_sparsity` --- src/diff.jl | 26 ++++++---------------- src/linearity.jl | 22 +++--------------- test/diff.jl | 58 ++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 68 insertions(+), 38 deletions(-) diff --git a/src/diff.jl b/src/diff.jl index 891e2cde2..de2a4c091 100644 --- a/src/diff.jl +++ b/src/diff.jl @@ -646,24 +646,13 @@ let linearity_rules = [ @rule +(~~xs) => reduce(+, filter(isidx, ~~xs), init=_scalar) @rule *(~~xs) => reduce(*, filter(isidx, ~~xs), init=_scalar) - @rule (~f)(~x::(!isidx)) => _scalar - @rule (~f)(~x::isidx) => if haslinearity_1(~f) - combine_terms_1(linearity_1(~f), ~x) - else - error("Function of unknown linearity used: ", ~f) - end + @rule (~f)(~x) => isidx(~x) ? combine_terms_1(linearity_1(~f), ~x) : _scalar @rule (^)(~x::isidx, ~y) => ~y isa Number && isone(~y) ? ~x : (~x) * (~x) - @rule (~f)(~x, ~y) => begin - if haslinearity_2(~f) - a = isidx(~x) ? ~x : _scalar - b = isidx(~y) ? ~y : _scalar - combine_terms_2(linearity_2(~f), a, b) - else - error("Function of unknown linearity used: ", ~f) - end - end - @rule ~x::issym => 0] + @rule (~f)(~x, ~y) => combine_terms_2(linearity_2(~f), isidx(~x) ? ~x : _scalar, isidx(~y) ? ~y : _scalar) + + @rule ~x::issym => 0 + ] linearity_propagator = Fixpoint(Postwalk(Chain(linearity_rules); maketerm=basic_mkterm)) global hessian_sparsity @@ -696,9 +685,8 @@ let @assert !(expr isa AbstractArray) expr = value(expr) u = map(value, vars) - idx(i) = TermCombination(Set([Dict(i=>1)])) - dict = Dict(u .=> idx.(1:length(u))) - f = Rewriters.Prewalk(x->haskey(dict, x) ? dict[x] : x; maketerm=basic_mkterm)(expr) + dict = Dict(ui => TermCombination(Set([Dict(i=>1)])) for (i, ui) in enumerate(u)) + f = Rewriters.Prewalk(x-> get(dict, x, x); maketerm=basic_mkterm)(expr) lp = linearity_propagator(f) S = _sparse(lp, length(u)) S = full ? S : tril(S) diff --git a/src/linearity.jl b/src/linearity.jl index 3fab3c487..3eae625cc 100644 --- a/src/linearity.jl +++ b/src/linearity.jl @@ -1,61 +1,45 @@ using SpecialFunctions import Base.Broadcast - -const linearity_known_1 = IdDict{Function,Bool}() -const linearity_known_2 = IdDict{Function,Bool}() - const linearity_map_1 = IdDict{Function, Bool}() const linearity_map_2 = IdDict{Function, Tuple{Bool, Bool, Bool}}() # 1-arg - const monadic_linear = [deg2rad, +, rad2deg, transpose, -, conj] const monadic_nonlinear = [asind, log1p, acsch, erfc, digamma, acos, asec, acosh, airybiprime, acsc, cscd, log, tand, log10, csch, asinh, airyai, abs2, gamma, lgamma, erfcx, bessely0, cosh, sin, cos, atan, cospi, cbrt, acosd, bessely1, acoth, erfcinv, erf, dawson, inv, acotd, airyaiprime, erfinv, trigamma, asecd, besselj1, exp, acot, sqrt, sind, sinpi, asech, log2, tan, invdigamma, airybi, exp10, sech, erfi, coth, asin, cotd, cosd, sinh, abs, besselj0, csc, tanh, secd, atand, sec, acscd, cot, exp2, expm1, atanh, slog, ssqrt, scbrt] -# We store 3 bools even for 1-arg functions for type stability -const three_trues = (true, true, true) for f in monadic_linear - linearity_known_1[f] = true linearity_map_1[f] = true end for f in monadic_nonlinear - linearity_known_1[f] = true linearity_map_1[f] = false end # 2-arg for f in [+, rem2pi, -, >, isless, <, isequal, max, min, convert, <=, >=] - linearity_known_2[f] = true linearity_map_2[f] = (true, true, true) end for f in [*] - linearity_known_2[f] = true linearity_map_2[f] = (true, true, false) end for f in [/] - linearity_known_2[f] = true linearity_map_2[f] = (true, false, false) end for f in [\] - linearity_known_2[f] = true linearity_map_2[f] = (false, true, false) end for f in [hypot, atan, mod, rem, lbeta, ^, beta] - linearity_known_2[f] = true linearity_map_2[f] = (false, false, false) end -haslinearity_1(@nospecialize(f)) = get(linearity_known_1, f, false) -haslinearity_2(@nospecialize(f)) = get(linearity_known_2, f, false) - -linearity_1(@nospecialize(f)) = linearity_map_1[f] -linearity_2(@nospecialize(f)) = linearity_map_2[f] +# Fallback assumption: Function is not linear, i.e., derivatives are non-zero +linearity_1(@nospecialize(f)) = get(linearity_map_1, f, false) +linearity_2(@nospecialize(f)) = get(linearity_map_2, f, (false, false, false)) # TermCombination datastructure diff --git a/test/diff.jl b/test/diff.jl index d40fa1185..6f0c97160 100644 --- a/test/diff.jl +++ b/test/diff.jl @@ -407,3 +407,61 @@ let @test isequal(expand_derivatives(D(Symbolics.scbrt(1 + x ^ 2))), simplify((2x) / (3Symbolics.scbrt(1 + x^2)^2))) @test isequal(expand_derivatives(D(Symbolics.slog(1 + x ^ 2))), simplify((2x) / (1 + x ^ 2))) end + +# Hessian sparsity involving unknown functions +let + @variables x₁ x₂ p q[1:1] + expr = 3x₁^2 + 4x₁ * x₂ + @test Matrix(Symbolics.hessian_sparsity(expr, [x₁, x₂])) == [true true; true false] + + expr = 3x₁^2 + 4x₁ * x₂ + p + @test Matrix(Symbolics.hessian_sparsity(expr, [x₁, x₂])) == [true true; true false] + + # issue 643: example test2_num + expr = 3x₁^2 + 4x₁ * x₂ + q[1] + @test Matrix(Symbolics.hessian_sparsity(expr, [x₁, x₂])) == [true true; true false] + + # Custom function: By default assumed to be non-linear + myexp(x) = exp(x) + @register_symbolic myexp(x) + expr = 3x₁^2 + 4x₁ * x₂ + myexp(p) + @test Matrix(Symbolics.hessian_sparsity(expr, [x₁, x₂])) == [true true; true false] + expr = 3x₁^2 + 4x₁ * x₂ + myexp(x₂) + @test Matrix(Symbolics.hessian_sparsity(expr, [x₁, x₂])) == [true true; true true] + + mylogaddexp(x, y) = log(exp(x) + exp(y)) + @register_symbolic mylogaddexp(x, y) + expr = 3x₁^2 + 4x₁ * x₂ + mylogaddexp(p, 2) + @test Matrix(Symbolics.hessian_sparsity(expr, [x₁, x₂])) == [true true; true false] + expr = 3x₁^2 + 4x₁ * x₂ + mylogaddexp(3, p) + @test Matrix(Symbolics.hessian_sparsity(expr, [x₁, x₂])) == [true true; true false] + expr = 3x₁^2 + 4x₁ * x₂ + mylogaddexp(p, 2) + @test Matrix(Symbolics.hessian_sparsity(expr, [x₁, x₂])) == [true true; true false] + expr = 3x₁^2 + 4x₁ * x₂ + mylogaddexp(p, q[1]) + @test Matrix(Symbolics.hessian_sparsity(expr, [x₁, x₂])) == [true true; true false] + expr = 3x₁^2 + 4x₁ * x₂ + mylogaddexp(p, x₂) + @test Matrix(Symbolics.hessian_sparsity(expr, [x₁, x₂])) == [true true; true true] + expr = 3x₁^2 + 4x₁ * x₂ + mylogaddexp(x₂, 4) + @test Matrix(Symbolics.hessian_sparsity(expr, [x₁, x₂])) == [true true; true true] + + # Custom linear function: Possible to extend `Symbolics.linearity_1`/`Symbolics.linearity_2` + myidentity(x) = x + @register_symbolic myidentity(x) + Symbolics.linearity_1(::typeof(myidentity)) = true + expr = 3x₁^2 + 4x₁ * x₂ + myidentity(p) + @test Matrix(Symbolics.hessian_sparsity(expr, [x₁, x₂])) == [true true; true false] + expr = 3x₁^2 + 4x₁ * x₂ + myidentity(q[1]) + @test Matrix(Symbolics.hessian_sparsity(expr, [x₁, x₂])) == [true true; true false] + expr = 3x₁^2 + 4x₁ * x₂ + myidentity(x₂) + @test Matrix(Symbolics.hessian_sparsity(expr, [x₁, x₂])) == [true true; true false] + + mymul1plog(x, y) = x * (1 + log(y)) + @register_symbolic mymul1plog(x, y) + Symbolics.linearity_2(::typeof(mymul1plog)) = (true, false, false) + expr = 3x₁^2 + 4x₁ * x₂ + mymul1plog(p, q[1]) + @test Matrix(Symbolics.hessian_sparsity(expr, [x₁, x₂])) == [true true; true false] + expr = 3x₁^2 + 4x₁ * x₂ + mymul1plog(x₂, q[1]) + @test Matrix(Symbolics.hessian_sparsity(expr, [x₁, x₂])) == [true true; true false] + expr = 3x₁^2 + 4x₁ * x₂ + mymul1plog(q[1], x₂) + @test Matrix(Symbolics.hessian_sparsity(expr, [x₁, x₂])) == [true true; true true] +end From 72b39f63cf1e14cbf3863fc31a4375359ee6b61a Mon Sep 17 00:00:00 2001 From: David Widmann Date: Fri, 13 Dec 2024 00:16:01 +0100 Subject: [PATCH 2/5] Add test of #555 --- test/diff.jl | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/test/diff.jl b/test/diff.jl index 6f0c97160..2d7906582 100644 --- a/test/diff.jl +++ b/test/diff.jl @@ -465,3 +465,21 @@ let expr = 3x₁^2 + 4x₁ * x₂ + mymul1plog(q[1], x₂) @test Matrix(Symbolics.hessian_sparsity(expr, [x₁, x₂])) == [true true; true true] end + +# issue #555 +let + # first example + @variables p[1:1] x[1:1] + p = collect(p) + x = collect(x) + @test collect(Symbolics.sparsehessian(p[1] * x[1], x)) == [0;;] + @test isequal(collect(Symbolics.sparsehessian(p[1] * x[1]^2, x)), [2p[1];;]) + + # second example + @variables a[1:2] + a = collect(a) + ex = (a[1]+a[2])^2 + @test Symbolics.hessian(ex, [a[1]]) == [2;;] + @test collect(Symbolics.sparsehessian(ex, [a[1]])) == [2;;] + @test collect(Symbolics.sparsehessian(ex, a)) == fill(2, 2, 2) +end From a764b96c22d10ce25e7f3d262e99ebb82b833b19 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Fri, 13 Dec 2024 00:20:20 +0100 Subject: [PATCH 3/5] Add test for #847 --- test/diff.jl | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/test/diff.jl b/test/diff.jl index 2d7906582..3600a9ba4 100644 --- a/test/diff.jl +++ b/test/diff.jl @@ -483,3 +483,18 @@ let @test collect(Symbolics.sparsehessian(ex, [a[1]])) == [2;;] @test collect(Symbolics.sparsehessian(ex, a)) == fill(2, 2, 2) end + +# issue #847 +let + @variables x[1:2] y[1:2] + x = Symbolics.scalarize(x) + y = Symbolics.scalarize(y) + + z = (x[1] + x[2]) * (y[1] + y[2]) + @test Symbolics.islinear(z, x) + @test Symbolics.isaffine(z, x) + + z = (x[1] + x[2]) + @test Symbolics.islinear(z, x) + @test Symbolics.isaffine(z, x) +end From 565d81af6ca9c41de342f2bdb10acb503fb4abf2 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Fri, 13 Dec 2024 00:44:52 +0100 Subject: [PATCH 4/5] Add test for #790 --- test/diff.jl | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/test/diff.jl b/test/diff.jl index 3600a9ba4..7cb2afb02 100644 --- a/test/diff.jl +++ b/test/diff.jl @@ -498,3 +498,15 @@ let @test Symbolics.islinear(z, x) @test Symbolics.isaffine(z, x) end + +# issue #790 +let + c(x) = [sum(x) - 1] + @variables xs[1:2] ys[1:1] + w = Symbolics.scalarize(xs) + v = Symbolics.scalarize(ys) + expr = dot(v, c(w)) + @test !Symbolics.islinear(expr, w) + @test Symbolics.isaffine(expr, w) + @test collect(Symbolics.hessian_sparsity(expr, w)) == fill(false, 2, 2) +end From 3185d301e6cf487e298c17a336baf0df3d6fdb14 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Fri, 13 Dec 2024 00:49:58 +0100 Subject: [PATCH 5/5] Add test for #749 --- test/diff.jl | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/test/diff.jl b/test/diff.jl index 7cb2afb02..fc29f3da5 100644 --- a/test/diff.jl +++ b/test/diff.jl @@ -510,3 +510,13 @@ let @test Symbolics.isaffine(expr, w) @test collect(Symbolics.hessian_sparsity(expr, w)) == fill(false, 2, 2) end + +# issue #749 +let + @variables x y + @register_symbolic Base.FastMath.exp_fast(x, y) + expr = Base.FastMath.exp_fast(x, y) + @test !Symbolics.islinear(expr, [x, y]) + @test !Symbolics.isaffine(expr, [x, y]) + @test collect(Symbolics.hessian_sparsity(expr, [x, y])) == fill(true, 2, 2) +end