L-curve analysis

This example shows how to perform L-curve analysis to estimate the optimal regularization weight for an optimization problem. We will use the Gauss-Newton method to solve a regularized optimization problem for a range of regularization weights, and then we will analyze the results to find the optimal weight.

Define the model function

First, we will define the model function that represents the curve we want to fit. In this example, we will use an exponential function defined as:

\[y = a \cdot e^{b \cdot x}\]

We will also generate synthetic data points based on this model function, and we will add some noise and bias to the data to make the problem more realistic.

import numpy as np
import matplotlib.pyplot as plt
import pysolvegn


np.random.seed(0)


# Define the model function
def model(params, x):
    a, b = params
    return a * np.exp(b * x)


# Generate synthetic data points with lot of noise
x_data = np.linspace(0, 3, 100)
true_params = [2.5, 0.5]  # True parameters for the curve: y = a * exp(b * x)
y_true = model(true_params, x_data)
y_data = (
    y_true + 2.0 * np.random.normal(size=y_true.shape) + 0.2
)  # Add noise and bias to the data

Create the residual and Jacobian functions

Secondly, we will define the residual function, and the Jacobian function. The residual function computes the difference between the observed data points and the model predictions, and the Jacobian function computes the derivatives of the residuals with respect to the parameters.

This equation is called a term in pysolve-gn, and it represents a single component of the optimization problem. See pysolvegn.Term for more details on how to define terms in pysolve-gn.

# Define the residual function
def residual_function(params, x, y):
    return model(params, x) - y


residual_func = lambda params: residual_function(params, x_data, y_data)


# Define the Jacobian function
def jacobian_function(params, x):
    a, b = params
    J = np.zeros((len(x), len(params)))
    J[:, 0] = np.exp(b * x)  # Derivative with respect to a
    J[:, 1] = a * x * np.exp(b * x)  # Derivative with respect to b
    return J


jacobian_func = lambda params: jacobian_function(params, x_data)

# Build the regularization term
estimated_params = [2.2, 0.55]
estimated_stds = [0.5, 0.2]
estimated_trust = [0.1, 0.05]

residual_reg_func, jacobian_reg_func = pysolvegn.build_soft_squared_regularization(
    means=estimated_params, stds=estimated_stds, thresholds=estimated_trust
)

data_terms = pysolvegn.Term.from_rJ(
    residual_func=residual_func, jacobian_func=jacobian_func, loss="linear", weight=1.0
)
reg_term = pysolvegn.Term.from_rJ(
    residual_func=residual_reg_func,
    jacobian_func=jacobian_reg_func,
    loss="linear",
    weight=1.0,
)

Perform L-curve analysis

Now we can perform L-curve analysis by solving the optimization problem for a range of regularization weights. We will use the solve function to solve the optimization problem for each weight, and we will store the results to analyze later.

weights = np.logspace(-2, 4, 50)  # Range of regularization weights to test

pysolvegn.perform_Lcurve_analysis(
    data_term=data_terms,
    reg_term=reg_term,
    x0=estimated_params,
    reg_weights=weights,
    max_iteration=10,
    xtol=1e-6,
    ftol=1e-6,
    verbosity=0,
    n_labels=20,
)