SciPy Optimization Bridge¶

In scientific machine learning and engineering applications, we frequently encounter inverse design and constrained optimization problems. For example, in concrete mix optimization, we want to find the precise proportions of cement, water, and aggregates that yield a target compressive strength while satisfying physical bounds (e.g., ingredient ratios between 0 and 1) and minimizing ingredient costs.

While sorix provides native neural network optimizers (like Adam and SGD) designed for updating weights, mathematical libraries like SciPy feature highly optimized, classical optimization algorithms that support complex bounds and non-linear constraints (such as L-BFGS-B and SLSQP).

The ScipyBridge utility bridges the gap between sorix autograd and SciPy's solvers, combining exact analytical gradients with classical constrained optimization.

1. Why use ScipyBridge?¶

When using SciPy's optimizers with machine learning models, developers typically face two major challenges:

1. Numerical Finite Differences are Slow and Unstable: By default, SciPy approximates gradients using finite differences (evaluating the objective function $2D$ times for $D$ variables). This is computationally expensive and numerically unstable (highly sensitive to step size $h$). ScipyBridge solves this by calculating exact analytical gradients via sorix backpropagation.
2. CPU-GPU Memory Boundary: SciPy is a CPU-bound library that expects CPU NumPy arrays. If your sorix model runs on GPU (via CuPy/CUDA), ScipyBridge automatically handles copying candidates onto the GPU, evaluating the model, and returning gradients back to CPU NumPy arrays seamlessly.

In [1]:

# Uncomment to install sorix
#!pip install 'sorix @ git+https://github.com/Mitchell-Mirano/sorix.git@main'

# Uncomment to install sorix #!pip install 'sorix @ git+https://github.com/Mitchell-Mirano/sorix.git@main'

2. Example 1: Basic Function Minimization¶

Let's start with a simple optimization problem: minimizing a 2D quadratic bowl function: $$f(x, y) = (x - 2)^2 + (y - 3)^2 + 5$$

We will initialize our variables at [0.0, 0.0] and use SciPy's L-BFGS-B algorithm to find the optimum at [2.0, 3.0].

In [2]:

import numpy as np
from sorix import tensor
from sorix.optim import ScipyBridge
import scipy.optimize

# 1. Define the parameters to optimize (requires_grad must be True)
params = tensor([0.0, 0.0], requires_grad=True)

# 2. Define the objective function (must return a scalar loss Tensor)
def loss_fn():
    return (params[0] - 2.0)**2 + (params[1] - 3.0)**2 + 5.0

# 3. Wrap with ScipyBridge
bridge = ScipyBridge(params, loss_fn)

# 4. Optimize using SciPy L-BFGS-B
res = scipy.optimize.minimize(
    bridge.objective,
    bridge.get_x(),
    jac=True,  # Tell SciPy that the objective function returns both (value, gradient)
    method='L-BFGS-B'
)

print("Optimization Successful:", res.success)
print("Optimal x, y:", res.x)
print("Minimum loss value:", res.fun)

import numpy as np from sorix import tensor from sorix.optim import ScipyBridge import scipy.optimize # 1. Define the parameters to optimize (requires_grad must be True) params = tensor([0.0, 0.0], requires_grad=True) # 2. Define the objective function (must return a scalar loss Tensor) def loss_fn(): return (params[0] - 2.0)**2 + (params[1] - 3.0)**2 + 5.0 # 3. Wrap with ScipyBridge bridge = ScipyBridge(params, loss_fn) # 4. Optimize using SciPy L-BFGS-B res = scipy.optimize.minimize( bridge.objective, bridge.get_x(), jac=True, # Tell SciPy that the objective function returns both (value, gradient) method='L-BFGS-B' ) print("Optimization Successful:", res.success) print("Optimal x, y:", res.x) print("Minimum loss value:", res.fun)

Optimization Successful: True
Optimal x, y: [2.00000022 2.99999989]
Minimum loss value: 5.0

3. Example 2: Bounded Concrete Mix Design Optimization¶

In physical engineering, concrete mix design involves combining three main components:

1. Cement ($x_0$): High cost, high carbon footprint, but provides the primary binding strength.
2. Water ($x_1$): Low cost, but excessive water decreases compressive strength.
3. Aggregates ($x_2$): Low cost, acts as structural filler.

Step 2.1: Generate Synthetic Concrete Data¶

Let's assume we collect historical lab data where we know the mix ratios of cement, water, and aggregates, and the corresponding measured Compressive Strength (in MPa). We simulate this data using a known linear combination with random noise.

In [3]:

import numpy as np
from sorix import tensor

np.random.seed(42)
num_samples = 200

# Generate mix fractions
cement = np.random.uniform(0.2, 0.6, (num_samples, 1))
water = np.random.uniform(0.1, 0.4, (num_samples, 1))
aggregates = np.random.uniform(0.3, 0.7, (num_samples, 1))

X_data = np.hstack([cement, water, aggregates])
# Simulating strength: increases with cement, decreases with water, slight aggregate impact
y_data = (40.0 * cement - 15.0 * water + 5.0 * aggregates + np.random.normal(0, 0.5, (num_samples, 1)))

X = tensor(X_data)
y = tensor(y_data)

import numpy as np from sorix import tensor np.random.seed(42) num_samples = 200 # Generate mix fractions cement = np.random.uniform(0.2, 0.6, (num_samples, 1)) water = np.random.uniform(0.1, 0.4, (num_samples, 1)) aggregates = np.random.uniform(0.3, 0.7, (num_samples, 1)) X_data = np.hstack([cement, water, aggregates]) # Simulating strength: increases with cement, decreases with water, slight aggregate impact y_data = (40.0 * cement - 15.0 * water + 5.0 * aggregates + np.random.normal(0, 0.5, (num_samples, 1))) X = tensor(X_data) y = tensor(y_data)

Step 2.2: Define and Train a Neural Network¶

Now we train a simple Multilayer Perceptron (MLP) in sorix to learn the mapping from ingredients [cement, water, aggregates] to predicted strength.

In [4]:

from sorix.nn import Linear, Sequential, MSELoss
from sorix.optim import Adam

# 1. Define model
model = Sequential(
    Linear(3, 8),
    Linear(8, 1)
)

# 2. Train the model
optimizer = Adam(model.parameters(), lr=0.01)
criterion = MSELoss()

model.train()
for epoch in range(101):
    optimizer.zero_grad()
    pred = model(X)
    loss = criterion(pred, y)
    loss.backward()
    optimizer.step()
    if epoch % 20 == 0:
        print(f"Epoch {epoch:3d} | Training Loss: {loss.item():.4f}")

from sorix.nn import Linear, Sequential, MSELoss from sorix.optim import Adam # 1. Define model model = Sequential( Linear(3, 8), Linear(8, 1) ) # 2. Train the model optimizer = Adam(model.parameters(), lr=0.01) criterion = MSELoss() model.train() for epoch in range(101): optimizer.zero_grad() pred = model(X) loss = criterion(pred, y) loss.backward() optimizer.step() if epoch % 20 == 0: print(f"Epoch {epoch:3d} | Training Loss: {loss.item():.4f}")

Epoch   0 | Training Loss: 202.2552
Epoch  20 | Training Loss: 129.1668
Epoch  40 | Training Loss: 62.2291
Epoch  60 | Training Loss: 24.8746
Epoch  80 | Training Loss: 20.1260
Epoch 100 | Training Loss: 19.5862

Step 2.3: Bounded Optimization with Multi-Objective Loss¶

Now, we lock the model weights (model.eval()) and optimize the inputs to find a concrete mix that:

1. Hits a target strength of 20.0 MPa.
2. Minimizes ingredient costs using pricing coefficients:
   • Cement: $0.15/unit
   • Water: $0.02/unit
   • Aggregates: $0.05/unit
3. Satisfies physical box bounds for each ingredient:
   • Cement: $[0.2, 0.6]$
   • Water: $[0.1, 0.4]$
   • Aggregates: $[0.3, 0.7]$

In [5]:

# 1. Set model to evaluation mode and freeze parameters
model.eval()
for p in model.parameters():
    p.requires_grad = False

# 2. Define input decision variable (2D tensor for batch prediction)
x_input = tensor([[0.4, 0.25, 0.5]], requires_grad=True)

# 3. Define target strength & cost coefficient vector C
target_strength = 20.0
cost_coeffs = tensor([0.15, 0.02, 0.05])

# 4. Define the multi-objective loss function
def objective_fn():
    pred = model(x_input)
    # Target strength loss (scalarized via .sum())
    strength_loss = ((pred - target_strength) ** 2).sum()
    
    # Vectorized ingredient cost calculation: sum(C * X)
    cost = (x_input * cost_coeffs).sum()
    # Combine losses (1.0 weight on cost penalty)
    return strength_loss + 1.0 * cost

# 5. Define boundary constraints
bounds = [
    (0.2, 0.6),  # Cement bounds
    (0.1, 0.4),  # Water bounds
    (0.3, 0.7)   # Aggregates bounds
]

# 6. Run optimization via ScipyBridge
bridge = ScipyBridge(x_input, objective_fn)
res = scipy.optimize.minimize(
    bridge.objective,
    bridge.get_x(),
    jac=True,
    bounds=bounds,
    method='L-BFGS-B'
)

print("Optimization Successful:", res.success)
print("Optimal Mix (Cement, Water, Aggregates):", res.x)
print("Predicted Compressive Strength (MPa):", model(tensor([res.x])).item())
print("Optimal Mix Cost ($):", (res.x * cost_coeffs.numpy()).sum())

# 1. Set model to evaluation mode and freeze parameters model.eval() for p in model.parameters(): p.requires_grad = False # 2. Define input decision variable (2D tensor for batch prediction) x_input = tensor([[0.4, 0.25, 0.5]], requires_grad=True) # 3. Define target strength & cost coefficient vector C target_strength = 20.0 cost_coeffs = tensor([0.15, 0.02, 0.05]) # 4. Define the multi-objective loss function def objective_fn(): pred = model(x_input) # Target strength loss (scalarized via .sum()) strength_loss = ((pred - target_strength) ** 2).sum() # Vectorized ingredient cost calculation: sum(C * X) cost = (x_input * cost_coeffs).sum() # Combine losses (1.0 weight on cost penalty) return strength_loss + 1.0 * cost # 5. Define boundary constraints bounds = [ (0.2, 0.6), # Cement bounds (0.1, 0.4), # Water bounds (0.3, 0.7) # Aggregates bounds ] # 6. Run optimization via ScipyBridge bridge = ScipyBridge(x_input, objective_fn) res = scipy.optimize.minimize( bridge.objective, bridge.get_x(), jac=True, bounds=bounds, method='L-BFGS-B' ) print("Optimization Successful:", res.success) print("Optimal Mix (Cement, Water, Aggregates):", res.x) print("Predicted Compressive Strength (MPa):", model(tensor([res.x])).item()) print("Optimal Mix Cost ($):", (res.x * cost_coeffs.numpy()).sum())

Optimization Successful: True
Optimal Mix (Cement, Water, Aggregates): [0.6 0.4 0.7]
Predicted Compressive Strength (MPa): 19.317100524902344
Optimal Mix Cost ($): 0.1330000039190054

4. Key Performance Guidelines¶

When using ScipyBridge, keep these important design tips in mind:

1. Call model.eval(): Always set your neural network to evaluation mode before running optimization. This disables dropout noise and freezes batch normalization running statistics, ensuring a smooth, noise-free optimization surface.
2. Freeze Weights: If you are optimizing the input features, remember to iterate over model.parameters() and set requires_grad = False to prevent computing gradients for network weights, saving memory and runtime.
3. Memory Management: ScipyBridge automatically calls backward(retain_graph=False) on the loss, freeing autograd graph memory at each step. This allows optimization runs to execute for hundreds of steps with constant memory usage.