在Python中实现类Matlab fmincon SQP的定向步长适配scipy basinhopping

Great question! I’ve dug into both SciPy’s basinhopping and MATLAB’s fmincon SQP implementations, so let’s walk through how to get those critical Zk and Gk parameters and build your custom take_step routine.

First: Clarify Zk and Gk #

From MATLAB’s constrained optimization docs:

• Zk: Null-space basis of the Jacobian of active constraints. This matrix defines directions that keep you within the feasible region (since they’re orthogonal to the gradients of constraints that are currently "active"—i.e., at their boundary or nearly violated).
• Gk: The reduced gradient of the Lagrangian function, which is the Lagrangian gradient projected onto Zk. This tells you the feasible direction to move to reduce the objective function.

How to Extract These in SciPy’s SLSQP #

SciPy’s SLSQP optimizer (the right local solver to pair with basinhopping for SQP-like constrained behavior) computes these matrices internally, but they’re not exposed in the public API. You can access them via the optimizer’s private state (note: pin your SciPy version to avoid breaking changes—this works reliably on v1.7+):

Step 1: Wrap the Local Optimizer to Capture State #

Create a custom local minimizer that runs SLSQP and saves Zk and Gk after each local optimization:

import numpy as np
from scipy.optimize import minimize, basinhopping

# Store local optimization state between basinhopping iterations
local_opt_state = {}

def custom_slsqp_minimizer(x0, **kwargs):
    # Run SLSQP local optimization
    res = minimize(x0=x0, method='SLSQP', **kwargs)
    
    # Access the internal SLSQP optimizer object
    slsqp = res.optimizer
    
    # Extract Zk: Null-space basis of active constraints (from SLSQP's QR decomposition)
    local_opt_state['Zk'] = slsqp.Z
    
    # Compute Gk: Reduced Lagrangian gradient (projected onto Zk)
    # Lagrangian gradient = objective gradient + sum(lambda_i * constraint gradients)
    lagrange_grad = slsqp.grad + np.dot(slsqp.lambda_, slsqp.jac)
    # Project onto null space to get the reduced gradient
    local_opt_state['Gk'] = np.dot(slsqp.Z.T, lagrange_grad)
    
    return res

Step 2: Build the Custom take_step Routine #

Now use the captured Zk and Gk to compute a directional step (replacing basinhopping’s default random perturbation):

def sqp_take_step(xk):
    # Fallback to random step if we don't have state yet (first iteration)
    if 'Zk' not in local_opt_state or 'Gk' not in local_opt_state:
        return xk + np.random.normal(0, 0.1, size=xk.shape)
    
    Zk = local_opt_state['Zk']
    Gk = local_opt_state['Gk']
    
    # Compute feasible descent direction: negative reduced gradient projected back to original space
    step_dir = np.dot(Zk, -Gk)
    
    # Optional: Add line search (like MATLAB's merit function line search)
    # For simplicity, we'll use a fixed step size here—replace with Armijo/Wolfe line search for better performance
    step_size = 0.1
    
    # Ensure the new point stays within bounds (if you defined them)
    new_xk = xk + step_size * step_dir
    if 'bounds' in local_opt_state:
        bounds = np.array(local_opt_state['bounds'])
        new_xk = np.clip(new_xk, bounds[:, 0], bounds[:, 1])
    
    return new_xk

Step 3: Run basinhopping with Your Custom Routine #

Hook everything up and test with an example objective/constraint:

# Example objective function
def objective(x):
    return (x[0]-1)**2 + (x[1]-2.5)**2

# Example equality constraint
def eq_constraint(x):
    return x[0]**2 + x[1]**2 - 4

# Initial point and variable bounds
x0 = np.array([0, 0])
bounds = ((-2, 2), (-2, 2))

# Store bounds in state for step clipping
local_opt_state['bounds'] = bounds

# Run basinhopping
result = basinhopping(
    objective,
    x0,
    take_step=sqp_take_step,
    local_minimizer=custom_slsqp_minimizer,
    minimizer_kwargs={
        'fun': objective,
        'constraints': {'type': 'eq', 'fun': eq_constraint},
        'bounds': bounds
    }
)

print("Basinhopping result:\n", result)

Critical Caveats to Keep in Mind #

• Private Attributes: Accessing slsqp.Z, slsqp.grad, etc., relies on SciPy’s internal implementation. If you upgrade SciPy, these names might change—check the SLSQP source code for your version if things break.
• Line Search: MATLAB’s SQP uses a merit function with line search to find the optimal step size. For better performance, replace the fixed step_size with a line search that balances objective reduction and constraint feasibility.
• Feasibility: The step direction from Zk ensures feasibility for linearized constraints, but you may need to adjust the new point to stay strictly feasible for non-linear constraints.
• Active Constraints: SciPy’s SLSQP automatically tracks active constraints, so the Z matrix you extract already matches MATLAB’s Zk (focused on the active constraint set).

内容的提问来源于stack exchange，提问作者Asher11