求助：从零推导正交线性回归梯度下降（基于NumPy）

Let’s break this down step by step—first, we’ll clarify the error function for orthogonal linear regression (it’s different from ordinary least squares!), then derive the gradients needed for gradient descent, and wrap it up with working NumPy code.

Ordinary least squares (OLS) minimizes the vertical distance from data points to the regression line (since it assumes only the y variable has noise). But orthogonal linear regression accounts for noise in both x and y, so we minimize the perpendicular distance from each point to the line y = W0 + W1*x.

For a single data point (x_i, y_i), the perpendicular distance to the line is:
$$d_i = \frac{|W1x_i - y_i + W0|}{\sqrt{W1^2 + 1}}$$

To make differentiation easier, we use the squared distance (absolute values are messy to take derivatives of). The total error function (scaled by 1/(2n) for cleaner gradient calculations) is:
$$J(W0, W1) = \frac{1}{2n} \sum_{i=1}^n \frac{(W1x_i - y_i + W0)^2}{W12 + 1}$$
The 1/2 cancels out the 2 from the power rule when taking derivatives, and n normalizes the loss over the number of samples.

We need the partial derivatives of J(W0, W1) with respect to W0 and W1—these tell us how to update our parameters to reduce the loss.

Partial Derivative with Respect to W0 #

Using the chain rule:
$$\frac{\partial J}{\partial W0} = \frac{1}{n} \sum_{i=1}^n \frac{W1x_i - y_i + W0}{W1^2 + 1}$$
The derivative of the numerator (W1x_i - y_i + W0)^2 with respect to W0 is 2(W1x_i - y_i + W0), and the 1/2 from the loss function cancels the 2. We divide by the denominator W1^2 +1 as-is.

Partial Derivative with Respect to W1 #

This one is trickier—we use the quotient rule (for f/g, derivative is (f’g - fg’)/g²):

• Let f = (W1x_i - y_i + W0)^2, so f’ = 2(W1x_i - y_i + W0)x_i
• Let g = W1^2 +1, so g’ = 2W1

Plugging into the quotient rule:
$$\frac{\partial}{\partial W1} \left( \frac{f}{g} \right) = \frac{2(W1x_i - y_i + W0)x_i(W1^2 +1) - 2W1(W1x_i - y_i + W0)^2}{(W12 +1)^2}$$

We can factor out 2(W1x_i - y_i + W0) from the numerator, then simplify the remaining terms:
$$2(W1x_i - y_i + W0) \left[ x_i(W1^2 +1) - W1(W1x_i - y_i + W0) \right]$$
Expanding the bracket:
$$x_iW1^2 + x_i - W1^2x_i + W1y_i - W1W0 = x_i + W1(y_i - W0)$$

Putting it all together (and canceling the 2 with the 1/2 from the loss function):
$$\frac{\partial J}{\partial W1} = \frac{1}{n} \sum_{i=1}^n \frac{(W1x_i - y_i + W0) \cdot (x_i + W1(y_i - W0))}{(W1^2 +1)^2}$$

Now we have everything to implement gradient descent:

• Initialize W0 and W1 (e.g., small random values)
• Choose a learning rate α (start small, like 0.01 or 0.05)
• Repeat until convergence (loss stops changing, or we hit a max number of iterations):
  1. Calculate the gradients ∂J/∂W0 and ∂J/∂W1 using all samples
  2. Update W0 = W0 - α * ∂J/∂W0
  3. Update W1 = W1 - α * ∂J/∂W1
  4. Check if the change in parameters (or loss) is below a small threshold (like 1e-6)

Here’s a complete Python code example with simulated data:

import numpy as np
import matplotlib.pyplot as plt

# Generate synthetic data (with noise in both x and y, perfect for orthogonal regression)
np.random.seed(42)
x = np.linspace(0, 10, 100)
true_W0 = 2
true_W1 = 3
# Add noise to both variables
y = true_W0 + true_W1 * x + np.random.normal(0, 1.5, size=x.shape)
x += np.random.normal(0, 0.5, size=x.shape)

def orthogonal_gradient_descent(x, y, learning_rate=0.01, epochs=10000, tol=1e-6):
    # Initialize parameters with small random values
    W0 = np.random.randn()
    W1 = np.random.randn()
    n = len(x)
    loss_history = []
    
    for _ in range(epochs):
        # Compute common terms to avoid redundant calculations
        numerator = W1 * x - y + W0
        denominator = W1**2 + 1
        
        # Calculate gradients
        dW0 = np.mean(numerator / denominator)
        term = x + W1 * (y - W0)
        dW1 = np.mean(numerator * term / (denominator**2))
        
        # Update parameters
        W0_new = W0 - learning_rate * dW0
        W1_new = W1 - learning_rate * dW1
        
        # Track loss
        loss = np.mean(numerator**2 / denominator) / 2
        loss_history.append(loss)
        
        # Check for convergence
        if np.abs(W0_new - W0) < tol and np.abs(W1_new - W1) < tol:
            break
        
        W0, W1 = W0_new, W1_new
    
    return W0, W1, loss_history

# Run the algorithm
W0_hat, W1_hat, loss_history = orthogonal_gradient_descent(x, y, learning_rate=0.05, epochs=20000)

# Print results
print(f"Estimated W0: {W0_hat:.4f}, Estimated W1: {W1_hat:.4f}")
print(f"True W0: {true_W0}, True W1: {true_W1}")

# Visualize the results
plt.figure(figsize=(12, 6))

# Plot data and regression lines
plt.subplot(1, 2, 1)
plt.scatter(x, y, alpha=0.6, label='Noisy Data')
plt.plot(x, W0_hat + W1_hat * x, color='red', linewidth=2, label='Orthogonal Regression Line')
plt.plot(x, true_W0 + true_W1 * x, color='green', linestyle='--', label='True Line')
plt.xlabel('X')
plt.ylabel('Y')
plt.legend()
plt.title('Orthogonal Regression Results')

# Plot loss over time
plt.subplot(1, 2, 2)
plt.plot(loss_history)
plt.xlabel('Epochs')
plt.ylabel('Loss')
plt.title('Loss Reduction During Gradient Descent')

plt.tight_layout()
plt.show()

When you run this, you’ll see the estimated parameters get very close to the true values, and the loss curve will flatten out as the algorithm converges.

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