qualia_core_db/specialized_libs/
multivar_calculus.rs

1//! **Multivariable symbolic differentiation** — gradient, Jacobian, Hessian (Calculus
2//! plan §3, the ★★ standout). Built on the CAS's existing single-variable
3//! [`differentiate`](super::symbolic_algebra::differentiate), so every partial is a
4//! *symbolic, provenance-bearing* derivative (citable via the CAS's `to_quins`/
5//! `expr_citation_hash`) — honest math, not a black-box autodiff number.
6//!
7//! Why this is the highest-demand gap: the learning spine needs it *now* — IRLS
8//! (logistic/Poisson GLM) needs the gradient + Hessian, the Bayesian Laplace
9//! approximation needs the Hessian of the log-posterior, and second-order optimisers
10//! need a Hessian. The symbolic forms here are differentiated once, then evaluated
11//! numerically at a point ([`gradient_at`]/[`hessian_at`]) for those consumers.
12//!
13//! No hot kernel (symbolic); the *numeric evaluation* of a gradient at scale is the
14//! bridge's `DenseLinear`/`ElementwiseMap` case.
15
16use std::collections::HashMap;
17
18use super::symbolic_algebra::{differentiate, simplify, Expr};
19
20/// `∂expr/∂var` — a single partial derivative (simplified). Thin alias over the CAS.
21pub fn partial(expr: &Expr, var: &str) -> Expr {
22    simplify(&differentiate(expr, var))
23}
24
25/// The **gradient** `∇f = [∂f/∂x₁, …, ∂f/∂xₙ]` as one simplified expression per variable.
26pub fn gradient(expr: &Expr, vars: &[&str]) -> Vec<Expr> {
27    vars.iter().map(|v| partial(expr, v)).collect()
28}
29
30/// The **Jacobian** of a vector of expressions: row `i` is `∇fᵢ`. Shape
31/// `exprs.len() × vars.len()`.
32pub fn jacobian(exprs: &[Expr], vars: &[&str]) -> Vec<Vec<Expr>> {
33    exprs.iter().map(|f| gradient(f, vars)).collect()
34}
35
36/// The **Hessian** `H[i][j] = ∂²f/∂xᵢ∂xⱼ` as an `n×n` matrix of simplified expressions.
37/// Symmetric by Clairaut's theorem (computed both ways implicitly via repeated
38/// differentiation).
39pub fn hessian(expr: &Expr, vars: &[&str]) -> Vec<Vec<Expr>> {
40    let grad = gradient(expr, vars);
41    grad.iter().map(|gi| gradient(gi, vars)).collect()
42}
43
44/// Evaluate the gradient numerically at `point` (variable → value). `None` if any
45/// partial fails to evaluate there (e.g. a division by zero in the domain).
46pub fn gradient_at(expr: &Expr, vars: &[&str], point: &HashMap<String, f64>) -> Option<Vec<f64>> {
47    gradient(expr, vars).iter().map(|g| g.eval(point)).collect()
48}
49
50/// Evaluate the Hessian numerically at `point`. `None` if any entry fails to evaluate.
51pub fn hessian_at(
52    expr: &Expr,
53    vars: &[&str],
54    point: &HashMap<String, f64>,
55) -> Option<Vec<Vec<f64>>> {
56    hessian(expr, vars)
57        .iter()
58        .map(|row| {
59            row.iter()
60                .map(|h| h.eval(point))
61                .collect::<Option<Vec<f64>>>()
62        })
63        .collect()
64}
65
66#[cfg(test)]
67mod tests {
68    use super::*;
69    use crate::specialized_libs::symbolic_algebra::{add, c, mul, pow, var};
70
71    fn env(pairs: &[(&str, f64)]) -> HashMap<String, f64> {
72        pairs.iter().map(|&(k, v)| (k.to_string(), v)).collect()
73    }
74
75    #[test]
76    fn gradient_of_a_quadratic_form() {
77        // f = x² + x·y + y²  →  ∇f = [2x + y, x + 2y]
78        let f = add(
79            add(pow(var("x"), 2), mul(var("x"), var("y"))),
80            pow(var("y"), 2),
81        );
82        let g = gradient(&f, &["x", "y"]);
83        let p = env(&[("x", 3.0), ("y", 5.0)]);
84        // ∂f/∂x = 2·3 + 5 = 11 ; ∂f/∂y = 3 + 2·5 = 13
85        assert!((g[0].eval(&p).unwrap() - 11.0).abs() < 1e-9);
86        assert!((g[1].eval(&p).unwrap() - 13.0).abs() < 1e-9);
87    }
88
89    #[test]
90    fn hessian_of_a_quadratic_is_constant() {
91        // f = x² + x·y + y²  →  H = [[2, 1], [1, 2]], symmetric & constant.
92        let f = add(
93            add(pow(var("x"), 2), mul(var("x"), var("y"))),
94            pow(var("y"), 2),
95        );
96        let h = hessian_at(&f, &["x", "y"], &env(&[("x", 0.0), ("y", 0.0)])).unwrap();
97        assert!((h[0][0] - 2.0).abs() < 1e-9);
98        assert!((h[0][1] - 1.0).abs() < 1e-9);
99        assert!((h[1][0] - 1.0).abs() < 1e-9); // symmetry
100        assert!((h[1][1] - 2.0).abs() < 1e-9);
101    }
102
103    #[test]
104    fn jacobian_shape_and_values() {
105        // F = [x·y, x + y] → J = [[y, x], [1, 1]]
106        let f1 = mul(var("x"), var("y"));
107        let f2 = add(var("x"), var("y"));
108        let j = jacobian(&[f1, f2], &["x", "y"]);
109        let p = env(&[("x", 2.0), ("y", 7.0)]);
110        assert_eq!(j.len(), 2);
111        assert!((j[0][0].eval(&p).unwrap() - 7.0).abs() < 1e-9); // ∂(xy)/∂x = y = 7
112        assert!((j[0][1].eval(&p).unwrap() - 2.0).abs() < 1e-9); // ∂(xy)/∂y = x = 2
113        assert!((j[1][0].eval(&p).unwrap() - 1.0).abs() < 1e-9);
114        assert!((j[1][1].eval(&p).unwrap() - 1.0).abs() < 1e-9);
115    }
116
117    #[test]
118    fn gradient_at_evaluates_numerically() {
119        // f = 3·x²  →  ∂f/∂x = 6x ; at x=4 → 24
120        let f = mul(c(3.0), pow(var("x"), 2));
121        let g = gradient_at(&f, &["x"], &env(&[("x", 4.0)])).unwrap();
122        assert!((g[0] - 24.0).abs() < 1e-9);
123    }
124}
qualia_core_db/specialized_libs/multivar_calculus.rs

qualia_core_db/specialized_libs/
multivar_calculus.rs
