entra 0.4.0

Entropy conserving transformations using divergence free vector fields

entra

Divergence-free tensor basis functions for incompressible vector field representation and entropy-conserving transformations.

---

Overview

This package implements a method for constructing divergence-free vector basis functions by applying a differential operator to Gaussian radial basis functions (RBFs). The resulting basis functions are ideal for representing incompressible vector fields and entropy-conserving transformations.

Theoretical Foundation

Maximum Entropy Principle

A fundamental theorem in information theory states that among all distributions with a given covariance matrix, the Gaussian distribution has maximum entropy. This package exploits a corollary of this theorem: by minimizing the covariance determinant of a point distribution using volume-preserving transformations, we can transform any distribution towards a Gaussian.

Why Divergence-Free?

Divergence-free vector fields are volume-preserving (incompressible). When we use divergence-free basis functions to define a transformation:

• The Jacobian determinant equals 1 everywhere
• Total volume is conserved
• Entropy is conserved under the transformation

This allows us to iteratively minimize the covariance determinant while preserving the fundamental entropy of the distribution, effectively reshaping any distribution into a Gaussian form.

The Divergence-Free Operator

The construction of divergence-free vector fields from scalar RBFs uses the differential operator discovered by Lowitzsch [1]:

Ô = -I∇² + ∇∇ᵀ

When applied to a scalar function φ(x), this operator produces a D×D matrix where each column is a divergence-free vector field.

References

[1] S. Lowitzsch, Approximation and Interpolation Employing Divergence-Free Radial Basis Functions With Applications, PhD thesis, Department of Mathematics, Texas A&M University, 2002.

Algorithm Flowchart

┌─────────────────────────────────────────────────────────────────────────────┐
│                         DIVERGENCE-FREE BASIS PIPELINE                       │
└─────────────────────────────────────────────────────────────────────────────┘

    ┌─────────────────────────────────────────┐
    │  INPUT PARAMETERS                       │
    │  • D: dimension                         │
    │  • δx: grid spacing                     │
    │  • n: points per dimension              │
    │  • σ = 0.7 × δx: RBF width             │
    └────────────────────┬────────────────────┘
                         │
                         ▼
    ╔═════════════════════════════════════════╗
    ║  STEP 1: GRID SAMPLING                  ║
    ║  VectorSampler                          ║
    ╠═════════════════════════════════════════╣
    ║  Sample J = n^D points on regular grid  ║
    ║                                         ║
    ║  x_{i} = x_center + i × δx              ║
    ║                                         ║
    ║  Output: eval_points (J, D)             ║
    ╚════════════════════┬════════════════════╝
                         │
                         ▼
    ╔═════════════════════════════════════════╗
    ║  STEP 2: CHOOSE L CENTERS               ║
    ╠═════════════════════════════════════════╣
    ║  Select L center points c_1, ..., c_L   ║
    ║  for basis functions                    ║
    ║                                         ║
    ║  Output: centers (L, D)                 ║
    ╚════════════════════┬════════════════════╝
                         │
                         ▼
    ╔═════════════════════════════════════════╗
    ║  STEP 3: SCALAR BASIS FUNCTIONS         ║
    ║  ScalarBasis                            ║
    ╠═════════════════════════════════════════╣
    ║  Gaussian RBF for each center:          ║
    ║                                         ║
    ║  φ_l(x) = exp(-||x - c_l||² / 2σ²)     ║
    ║                                         ║
    ║  Output: φ values (J, L)                ║
    ╚════════════════════┬════════════════════╝
                         │
                         ▼
    ╔═════════════════════════════════════════╗
    ║  STEP 4: APPLY TENSOR OPERATOR          ║
    ║  TensorBasis                            ║
    ╠═════════════════════════════════════════╣
    ║  Apply Ô = -I∇² + ∇∇ᵀ to each φ_l      ║
    ║                                         ║
    ║  Φ_l(x) = Ô φ_l(x)                      ║
    ║                                         ║
    ║  Output: Φ values (J, L, D, D)          ║
    ║          ───────────────────            ║
    ║          D×D matrix at each point       ║
    ╚════════════════════┬════════════════════╝
                         │
                         ▼
    ╔═════════════════════════════════════════╗
    ║  STEP 5: EXTRACT COLUMN VECTOR FIELDS   ║
    ╠═════════════════════════════════════════╣
    ║  Each column d of Φ is a vector field:  ║
    ║                                         ║
    ║  V_d = Φ[:, :, d]  shape: (J, D)        ║
    ║                                         ║
    ║  For D=2:                               ║
    ║  • Column 0: V_0 = [Φ_00, Φ_10]ᵀ       ║
    ║  • Column 1: V_1 = [Φ_01, Φ_11]ᵀ       ║
    ╚════════════════════┬════════════════════╝
                         │
                         ▼
    ╔═════════════════════════════════════════╗
    ║  STEP 6: VERIFY DIVERGENCE-FREE         ║
    ╠═════════════════════════════════════════╣
    ║  Compute divergence of each column:     ║
    ║                                         ║
    ║  div(V_d) = ∂V_{d,x}/∂x + ∂V_{d,y}/∂y  ║
    ║                                         ║
    ║  Each column should satisfy:            ║
    ║  div(V_d) ≈ 0                           ║
    ╚════════════════════┬════════════════════╝
                         │
                         ▼
    ┌─────────────────────────────────────────┐
    │  OUTPUT                                 │
    │  • D divergence-free vector fields      │
    │    per center (L × D total)             │
    │  • Each field has J discrete values     │
    └─────────────────────────────────────────┘

Covariance Optimization Algorithm

The optimization proceeds in two stages to transform any distribution towards Gaussian form while conserving entropy.

┌─────────────────────────────────────────────────────────────────────────────┐
│                    COVARIANCE OPTIMIZATION PIPELINE                          │
└─────────────────────────────────────────────────────────────────────────────┘

    ┌─────────────────────────────────────────┐
    │  INPUT                                  │
    │  • Points: (J, D) array                 │
    │  • σ: RBF width parameter               │
    │  • Target: H(uniform) = ln(V)           │
    └────────────────────┬────────────────────┘
                         │
                         ▼
    ╔═════════════════════════════════════════╗
    ║  STAGE 1: TENSOR BASIS OPTIMIZATION     ║
    ║  (L × D parameters)                     ║
    ╠═════════════════════════════════════════╣
    ║                                         ║
    ║  1. Build tensor basis Φ (J, L, D, D)   ║
    ║                                         ║
    ║  2. Define transformation:              ║
    ║     y' = y + Σ_l Φ_l(y) · c_l           ║
    ║     where c_l ∈ ℝ^D (L×D params)        ║
    ║                                         ║
    ║  3. Levenberg-Marquardt optimization:   ║
    ║     minimize det(Cov(y'))               ║
    ║                                         ║
    ║  4. Collapse to effective basis:        ║
    ║     V_l = Φ_l · c_l  →  (J, L, D)       ║
    ║                                         ║
    ╚════════════════════┬════════════════════╝
                         │
                         ▼
    ╔═════════════════════════════════════════╗
    ║  STAGE 2: OUTER LOOP REFINEMENT         ║
    ║  (L parameters per round)               ║
    ╠═════════════════════════════════════════╣
    ║                                         ║
    ║  For each outer round (typically 5):    ║
    ║                                         ║
    ║  1. Define scalar-weighted transform:   ║
    ║     y' = y + Σ_l α_l V_l                ║
    ║     where α_l ∈ ℝ (L scalar params)     ║
    ║                                         ║
    ║  2. Optimize: minimize det(Cov(y'))     ║
    ║                                         ║
    ║  3. Update basis: V_l ← α_l V_l         ║
    ║                                         ║
    ║  4. Transform points: y ← y'            ║
    ║                                         ║
    ║  5. Repeat with updated points & basis  ║
    ║                                         ║
    ╚════════════════════┬════════════════════╝
                         │
                         ▼
    ┌─────────────────────────────────────────┐
    │  OUTPUT                                 │
    │  • Transformed points (J, D)            │
    │  • H(Gaussian) ≈ H(uniform) = Target    │
    │  • Distribution → Gaussian form         │
    └─────────────────────────────────────────┘

    ┌─────────────────────────────────────────┐
    │  CONVERGENCE CRITERION                  │
    │                                         │
    │  Gap = H(Gaussian) - H(uniform) → 0     │
    │                                         │
    │  When gap ≈ 0, the distribution has     │
    │  been successfully Gaussianized while   │
    │  preserving the original entropy.       │
    └─────────────────────────────────────────┘

Why Two Stages?

Stage 1 uses L × D parameters (matrix-valued coefficients) for maximum flexibility in finding the initial transformation direction. This explores a larger parameter space but is computationally expensive.

Stage 2 uses only L parameters (scalar coefficients) per round. Each round operates on the transformed points from the previous round, allowing incremental refinements that compound across rounds. The effective basis adapts at each round, capturing increasingly fine-grained displacement patterns.

Full Optimization Example

import numpy as np
from entra import (
    VectorSampler, TensorBasis, Transformation,
    CovarianceMinimizer, EffectiveTransformation,
    EffectiveCovarianceMinimizer,
    shannon_entropy_gaussian, shannon_entropy_uniform
)

# Generate 2D uniform distribution (20×20 = 400 points)
sampler = VectorSampler(center=[0.0, 0.0], delta_x=1, num_points_per_dim=20)
points = sampler.sample()

# Create centers along axes
centers = np.array([[i, 0], [-i, 0], [0, i], [0, -i]]
                   for i in range(10)).reshape(-1, 2)

# Target entropy
target_H = shannon_entropy_uniform(points)
print(f"Target H(uniform): {target_H:.4f} nats")

# ============================================================
# STAGE 1: Tensor Basis Optimization
# ============================================================
sigma = 4.0  # Optimized RBF width
basis = TensorBasis(centers, sigma=sigma)
transformation = Transformation(basis)
minimizer = CovarianceMinimizer(transformation, points)

# Get initial coefficients and optimize
x = transformation.get_coefficients_flat().copy()
n_params = len(x)
lam, eps, tol = 1.0, 1e-7, 1e-12

for iteration in range(1000):
    r = minimizer.residuals_for_lm(x)

    # Numerical Jacobian
    J_mat = np.zeros((len(r), n_params))
    for i in range(n_params):
        x_plus = x.copy()
        x_plus[i] += eps
        J_mat[:, i] = (minimizer.residuals_for_lm(x_plus) - r) / eps

    # LM update
    JTJ = J_mat.T @ J_mat
    JTr = J_mat.T @ r
    delta = np.linalg.solve(JTJ + lam * np.eye(n_params), -JTr)

    x_new = x + delta
    if minimizer.objective_logdet(x_new) < minimizer.objective_logdet(x):
        improvement = abs(minimizer.objective_logdet(x) - minimizer.objective_logdet(x_new))
        x = x_new
        lam *= 0.1
        if improvement < tol:
            break
    else:
        lam *= 10.0

transformation.set_coefficients_flat(x)
stage1_points = transformation.transform(points)

cov = np.cov(stage1_points, rowvar=False)
print(f"After Stage 1: H(Gaussian) = {shannon_entropy_gaussian(cov):.4f}")

# ============================================================
# STAGE 2: Outer Loop Refinement (5 rounds)
# ============================================================
current_basis = transformation.get_updated_basis(points)
current_points = stage1_points.copy()

for round_num in range(1, 6):
    eff_transform = EffectiveTransformation(current_basis, current_points)
    eff_minimizer = EffectiveCovarianceMinimizer(eff_transform)

    result = eff_minimizer.optimize(max_iterations=1000, tolerance=1e-10)

    current_points = eff_transform.transform()
    current_basis = eff_transform.get_updated_basis()

    cov = np.cov(current_points, rowvar=False)
    H = shannon_entropy_gaussian(cov)
    gap = H - target_H
    print(f"Round {round_num}: H(Gaussian) = {H:.4f}, Gap = {gap:+.4f}")

print(f"\nFinal gap: {gap:+.4f} nats (target: 0)")

Example Results

Initial uniform distribution (before transformation):

Final Gaussian distribution (after 5 outer loop rounds):

Optimization history showing convergence:

Mathematical Details

For the complete mathematical formulation including:

• Grid sampling formulas
• Gaussian RBF definitions
• Tensor operator derivation
• Proof of divergence-free property
• Discrete divergence computation

See docs/equations.rst

API Reference

Classes

Class	Description
VectorSampler	Sample points on a regular D-dimensional grid
ScalarBasis	Gaussian RBF scalar basis functions
TensorBasis	Apply tensor operator to produce D×D matrix basis

Functions

Function	Description
divergence(V, dx, grid_shape)	Compute divergence of vector field (signed sum)
divergence_components(V, dx, grid_shape)	Get individual ∂V_d/∂x_d terms (signed)
gradient_component(f, dx, axis, grid_shape)	Compute ∂f/∂x_axis for scalar field
is_divergence_free(V, dx, grid_shape)	Check if field is divergence-free
tensor_basis_column_divergence(Phi, dx, grid_shape)	Divergence of all tensor basis columns
verify_tensor_basis_divergence_free(Phi, dx, grid_shape)	Verify all columns are divergence-free

Note: All divergence computations preserve signs. Use divergence_components() to see how individual terms (which may be positive or negative) sum to give the total divergence.

Shape Conventions

Symbol	Meaning
D	Dimension of space
J	Number of evaluation points (grid points)
L	Number of basis function centers

Array	Shape	Description
eval_points	(J, D)	Evaluation grid points
centers	(L, D)	Basis function centers
ScalarBasis.evaluate()	(J, L)	Scalar basis values
TensorBasis.evaluate()	(J, L, D, D)	Tensor basis matrices
divergence()	(J,)	Divergence at each point

Examples

See the notebooks/ directory for Jupyter notebook examples:

• 2d_vector_field_example.ipynb - 2D visualization of tensor basis columns

See the scripts/ directory for standalone demos:

• basis_pipeline_demo.py - Full pipeline with shape verification

Installation

You can install entra via pip:

pip install entra

To install latest development version:

pip install git+https://github.com/Kapoorlabs-CAPED/entra.git

Contributing

Contributions are very welcome. Tests can be run with tox, please ensure the coverage at least stays the same before you submit a pull request.

License

Distributed under the terms of the BSD-3 license, "entra" is free and open source software

Issues

If you encounter any problems, please file an issue along with a detailed description.