Sparsity

Understanding the sparsity pattern of your Jacobian and Hessian matrices is crucial for high-performance optimization. Janus provides tools to inspect sparsity patterns from symbolic graphs, visualize them as ASCII/PDF/HTML spy plots, compile sparse derivative evaluators that return only structural nonzeros, and surface CasADi graph coloring metadata. Sparsity analysis works in symbolic mode; the NaN-propagation fallback also supports numeric mode for black-box functions.

Quick Start

#include <janus/janus.hpp>
 
auto x = janus::sym("x", 10);
auto f = janus::SymbolicScalar::vertcat({
    x(1) - x(0),
    x(2) - janus::sin(x(1)),
    x(3) - x(2)
});
 
// Extract Jacobian sparsity pattern
janus::SparsityPattern sp = janus::sparsity_of_jacobian(f, x);
 
// Print ASCII spy plot
std::cout << sp.to_string() << std::endl;
 
// Build a sparse Jacobian evaluator (returns only nonzeros)
auto J = janus::sparse_jacobian(f, x);
std::cout << "nnz = " << J.nnz() << "\n";

Core API

The core class is janus::SparsityPattern in <janus/core/Sparsity.hpp>.

Extracting Sparsity

Function	Description
janus::sparsity_of_jacobian(f, x)	Jacobian sparsity from symbolic expressions
janus::sparsity_of_hessian(f, x)	Hessian sparsity from symbolic expressions
janus::get_jacobian_sparsity(func, out_idx, in_idx)	Jacobian sparsity from a compiled janus::Function
janus::get_hessian_sparsity(func, out_idx, in_idx)	Hessian sparsity from a compiled janus::Function
janus::nan_propagation_sparsity(callable, n_in, n_out)	Black-box sparsity via NaN propagation
janus::nan_propagation_sparsity(func)	NaN-propagation sparsity from janus::Function

Querying a SparsityPattern

Method	Description
sp.nnz()	Number of structural nonzeros
sp.density()	Fraction of nonzeros
sp.n_rows() / sp.n_cols()	Dimensions
sp.get_triplet()	Row/column index vectors
sp.get_ccs()	Compressed column storage

Sparse Derivative Evaluators

Function	Description
janus::sparse_jacobian(f, x)	Build a sparse Jacobian evaluator from expressions
janus::sparse_hessian(phi, x)	Build a sparse Hessian evaluator from a scalar expression
janus::sparse_jacobian(func, out_idx, in_idx)	Build from a janus::Function block
janus::sparse_hessian(func, out_idx, in_idx)	Build from a janus::Function block

Visualization

Method	Description
sp.to_string()	ASCII spy plot for terminal output
sp.visualize_spy(filename)	Render spy plot to PDF (requires Graphviz)
sp.export_spy_html(filename, title)	Interactive HTML spy plot with pan/zoom

Usage Patterns

From Symbolic Expressions

auto x = janus::sym("x", 10);
auto f = ...; // some expression depending on x
 
// Jacobian sparsity (df/dx)
janus::SparsityPattern J_sp = janus::sparsity_of_jacobian(f, x);
 
// Hessian sparsity (d^2 f/dx^2)
janus::SparsityPattern H_sp = janus::sparsity_of_hessian(f, x);

From a Compiled Function

janus::Function func(inputs, outputs);
auto sp = janus::get_jacobian_sparsity(func);

For multi-input or multi-output functions, use explicit block selection:

auto J_sp = janus::get_jacobian_sparsity(func, output_idx, input_idx);
auto H_sp = janus::get_hessian_sparsity(func, scalar_output_idx, input_idx);

From CasADi Types Directly

casadi::Sparsity raw = ...;
janus::SparsityPattern sp(raw);
 
janus::SymbolicScalar expr = ...;
janus::SparsityPattern expr_sp(expr); // Extract from MX sparsity

Sparse Jacobian Pipeline

auto x = janus::sym("x", 6);
auto f = janus::SymbolicScalar::vertcat({
    x(1) - x(0),
    x(2) - janus::sin(x(1)),
    x(3) - x(2)
});
 
auto J = janus::sparse_jacobian(f, x);
 
std::cout << "nnz = " << J.nnz() << "\n";
std::cout << "forward colors = " << J.forward_coloring().n_colors() << "\n";
std::cout << "reverse colors = " << J.reverse_coloring().n_colors() << "\n";
std::cout << "preferred mode = "
          << (J.preferred_mode() == janus::SparseJacobianMode::Forward ? "forward" : "reverse")
          << "\n";
 
janus::NumericMatrix x_val(6, 1);
x_val << 0.0, 0.2, 0.5, 0.9, 0.0, 0.0;
 
janus::NumericMatrix jac_nz = J.values(x_val);

jac_nz is a column vector of derivative values in the same CCS ordering reported by J.sparsity().get_triplet() and J.sparsity().get_ccs(). That ordering is fixed, so the sparsity structure can be reused across many evaluations.

Sparse Hessian Pipeline

auto x = janus::sym("x", 5);
janus::SymbolicScalar phi = 0;
for (int k = 0; k < 4; ++k) {
    auto diff = x(k + 1) - x(k);
    phi = phi + diff * diff;
}
 
auto H = janus::sparse_hessian(phi, x);
 
std::cout << "nnz = " << H.nnz() << "\n";
std::cout << "star colors = " << H.coloring().n_colors() << "\n";
 
janus::NumericMatrix x_val(5, 1);
x_val << 0.0, 0.1, 0.3, 0.7, 1.0;
 
janus::NumericMatrix hess_nz = H.values(x_val);

For Hessians, Janus exposes CasADi's star coloring through H.coloring().

From Function Blocks

Sparse derivative evaluators can also be built from an already-compiled janus::Function, selecting a specific output block and input block:

janus::Function terms("terms", {x, u}, {defects, objective});
 
auto ddefects_dx = janus::sparse_jacobian(terms, 0, 0);
auto ddefects_du = janus::sparse_jacobian(terms, 0, 1);
auto hobjective_xx = janus::sparse_hessian(terms, 1, 0);

This is the most useful form for optimization pipelines where one compiled function already exposes multiple residual, constraint, and objective blocks.

Reconstructing a Dense Matrix

When debugging, it is often useful to reconstruct the dense matrix from sparse values:

auto nz = J.values(x_val);
auto [rows, cols] = J.sparsity().get_triplet();
 
janus::NumericMatrix dense =
    janus::NumericMatrix::Zero(J.sparsity().n_rows(), J.sparsity().n_cols());
 
for (Eigen::Index k = 0; k < nz.size(); ++k) {
    dense(rows[static_cast<size_t>(k)], cols[static_cast<size_t>(k)]) = nz(k);
}

In production you usually keep the structural ordering and pass the nonzero vector straight into a sparse solver or downstream callback.

Visualization Workflows

ASCII spy plot for quick terminal debugging:

std::cout << sp.to_string() << std::endl;

Output:

Sparsity: 10x10, nnz=28 (density=28.000%)
+----------+
|**........|
|***.......
|.***......|
...
+----------+

PDF rendering for reports:

sp.visualize_spy("my_pattern"); // Creates my_pattern.pdf

Interactive HTML for exploring large matrices:

sp.export_spy_html("my_pattern", "My Jacobian"); // Creates my_pattern.html

The HTML output includes pan/zoom, clickable cells with row/col details, axis labels, and a stats panel.

Advanced Usage

Reading Coloring Metadata

janus::GraphColoring exposes:

• n_entries() for the uncompressed derivative size
• n_colors() for the compressed directional count
• compression_ratio() for a quick summary
• colorvec() for the per-entry color assignment

This is useful when comparing derivative blocks and deciding whether sparse directional sweeps are worth it.

NaN-Propagation Sparsity Detection

Sometimes you have black-box functions where symbolic sparsity analysis is not possible (external library calls, non-traceable operations, functions with runtime branching). Janus provides nan_propagation_sparsity() for these cases.

How it works:

1. Evaluate f(x) at a reference point
2. For each input i: set x[i] = NaN, evaluate f(x)
3. If output[j] becomes NaN, then it depends on input i, so Jacobian(j, i) is nonzero

// For lambda/callable functions
auto sp = janus::nan_propagation_sparsity(
    [](const janus::NumericVector& x) {
        janus::NumericVector y(x.size());
        for (int i = 0; i < x.size(); ++i) y(i) = x(i) * x(i);
        return y;
    },
    n_inputs, n_outputs);
 
// For janus::Function
janus::Function fn(...);
auto sp = janus::nan_propagation_sparsity(fn);
 
// With custom options (reference point)
janus::NaNSparsityOptions opts;
opts.reference_point = janus::NumericVector{{1.0, 2.0, 3.0}};
auto sp = janus::nan_propagation_sparsity(fn, opts);

Verifying symbolic sparsity:

auto x = janus::sym("x", 4);
auto f = x * x;  // Element-wise square
 
// Symbolic sparsity
auto sp_symbolic = janus::sparsity_of_jacobian(f, x);
 
// NaN-propagation sparsity (black-box equivalent)
janus::Function fn({x}, {f});
auto sp_nan = janus::nan_propagation_sparsity(fn);
 
// They should match!
assert(sp_symbolic == sp_nan);

Remarks

Use NaN-propagation to validate that your symbolic expressions have the expected structure, or when interfacing with external numerical code.

Example Walkthrough: sparsity_intro.cpp

The example examples/intro/sparsity_intro.cpp demonstrates four common structures found in optimization:

1. Simple Jacobian – Shows how dependencies map to nonzeros. f_0 = x_0^2 depends only on x_0, so row 0 has a nonzero at column 0.
2. Chain Structure (Tridiagonal) – f = sum((x[i] - x[i+1])^2) creates a tridiagonal Hessian. This band structure is typical in trajectory optimization where each state depends only on its neighbors.
3. Independent Systems (Block Diagonal) – Two completely separate systems stacked together form a block-diagonal matrix. Solvers can parallelize this trivially.
4. 2D Laplacian (5-Point Stencil) – Typical in PDE constraints. Each node depends on itself and its 4 neighbors. Uses janus::sym_vec_pair for 2D indexing:

auto [x_vec, x_mx] = janus::sym_vec_pair("x", n_vars);
// Build equations using x_vec(k), create Function using raw x_mx
janus::Function f_pde({x_mx}, {janus::SymbolicScalar::vertcat(eqs)});

Example Walkthrough: sparse_derivative_pipeline.cpp

The example examples/math/sparse_derivative_pipeline.cpp shows the full workflow:

1. Build a trajectory-style residual with local coupling
2. Compile sparse Jacobian blocks and a sparse Hessian block
3. Print nnz versus dense size
4. Inspect forward, reverse, and star coloring counts
5. Reuse the same sparse kernels at two different evaluation points

Build and run:

ninja -C build sparse_derivative_pipeline
./build/examples/sparse_derivative_pipeline

See Also

• Graph Visualization Guide – Visualize the computational graph itself
• Structural Diagnostics Guide – Use sparsity for observability and identifiability analysis
• sparsity_intro.cpp – Introductory sparsity patterns example
• sparse_derivative_pipeline.cpp – Full sparse derivative pipeline example
• Sparsity.hpp – API reference