沿着山坡寻找最低点

import math
import random


def loss(x, y, curvature=8.0):
    """
    二维正定二次目标函数：

    f(x,y) = 0.5 * (x*x + curvature*y*y)
    """
    return 0.5 * (
        x*x + curvature*y*y
    )


def gradient(x, y, curvature=8.0):
    """
    ∂f/∂x = x
    ∂f/∂y = curvature*y
    """
    return x, curvature*y


def gradient_descent(
    start=(-2.4, 1.4),
    lr=0.1,
    steps=30,
    curvature=8.0,
):
    """
    固定学习率梯度下降。
    """

    x, y = start

    history = []

    for step in range(steps + 1):
        history.append({
            "step": step,
            "x": x,
            "y": y,
            "loss": loss(x, y, curvature),
        })

        gx, gy = gradient(
            x, y, curvature
        )

        x -= lr * gx
        y -= lr * gy

        if not (
            math.isfinite(x)
            and math.isfinite(y)
        ):
            break

    return history

for lr in (0.1, 0.25, 0.3):

    history = gradient_descent(
        start=(-2.4, 1.4),
        lr=lr,
        steps=15,
        curvature=8.0,
    )

    final = history[-1]

    print(
        f"lr={lr:.2f}",
        f"x={final['x']:.5f}",
        f"y={final['y']:.5f}",
        f"loss={final['loss']:.6f}",
    )

def exact_line_search(
    start=(-2.4, 1.4),
    steps=30,
    curvature=8.0,
):
    x, y = start
    history = []

    for step in range(steps + 1):
        history.append((
            step, x, y,
            loss(x, y, curvature)
        ))

        gx, gy = gradient(x, y, curvature)

        numerator = gx*gx + gy*gy

        if numerator < 1e-20:
            break

        denominator = (
            gx*gx + curvature*gy*gy
        )

        lr = numerator / denominator

        x -= lr * gx
        y -= lr * gy

    return history


result = exact_line_search()

for step, x, y, value in result[:8]:
    print(
        step,
        round(x, 5),
        round(y, 5),
        round(value, 7)
    )

def armijo_step(
    x,
    y,
    curvature=8.0,
    initial_lr=1.0,
    shrink=0.5,
    sufficient_decrease=1e-4,
):
    """
    在负梯度方向上执行Armijo回溯线搜索。

    使用：
        f(new) <= f(old) - c*lr*||g||^2
    """

    gx, gy = gradient(x, y, curvature)

    norm2 = gx*gx + gy*gy

    if norm2 < 1e-20:
        return 0.0

    old_loss = loss(x, y, curvature)

    lr = initial_lr

    for _ in range(60):

        nx = x - lr*gx
        ny = y - lr*gy

        new_loss = loss(
            nx, ny, curvature
        )

        threshold = (
            old_loss
            - sufficient_decrease*lr*norm2
        )

        if new_loss <= threshold:
            return lr

        lr *= shrink

    raise RuntimeError(
        "回溯线搜索未找到合适步长"
    )


def descent_with_armijo(
    start=(-2.4, 1.4),
    steps=30,
):
    x, y = start
    history = []

    for step in range(steps + 1):
        history.append(
            (step, x, y, loss(x, y))
        )

        lr = armijo_step(x, y)

        if lr == 0:
            break

        gx, gy = gradient(x, y)

        x -= lr*gx
        y -= lr*gy

    return history


history = descent_with_armijo()

assert all(
    history[i+1][3] <= history[i][3] + 1e-12
    for i in range(len(history)-1)
)

print("Armijo线搜索测试通过")

def gradient_descent_momentum(
    start=(-2.4, 1.4),
    lr=0.06,
    beta=0.8,
    steps=30,
    curvature=8.0,
):
    x, y = start

    vx = 0.0
    vy = 0.0

    history = []

    for step in range(steps + 1):

        history.append({
            "step": step,
            "x": x,
            "y": y,
            "loss": loss(x, y, curvature),
        })

        gx, gy = gradient(
            x, y, curvature
        )

        # 梯度的指数累积
        vx = beta * vx + gx
        vy = beta * vy + gy

        x -= lr * vx
        y -= lr * vy

    return history


plain = gradient_descent(
    lr=0.06,
    steps=30,
)

momentum = gradient_descent_momentum(
    lr=0.06,
    beta=0.8,
    steps=30,
)

print("普通GD最终损失:", plain[-1]["loss"])
print("Momentum最终损失:", momentum[-1]["loss"])

beta=0.5

beta=0.9

(1, 2)
(2, 4)
(3, 6)
(4, 8)

def train_sgd(
    samples,
    initial_weight=0.0,
    lr=0.02,
    steps=400,
    seed=42,
):
    rng = random.Random(seed)

    w = initial_weight

    history = []

    for step in range(steps):

        x, y = rng.choice(samples)

        prediction = w * x

        # 单样本平方损失梯度
        grad = (prediction - y) * x

        w -= lr * grad

        if step % 20 == 0:
            history.append((step, w))

    return w, history


samples = [
    (1.0, 2.0),
    (2.0, 4.0),
    (3.0, 6.0),
    (4.0, 8.0),
]

weight, history = train_sgd(samples)

print("训练得到的w:", weight)

for step, w in history:
    print(
        f"step={step:3d}, w={w:.6f}"
    )

assert abs(weight - 2.0) < 1e-3

print("SGD训练完成")

def numerical_gradient(
    x,
    y,
    curvature=8.0,
    h=1e-5,
):
    gx = (
        loss(x+h, y, curvature)
        - loss(x-h, y, curvature)
    ) / (2*h)

    gy = (
        loss(x, y+h, curvature)
        - loss(x, y-h, curvature)
    ) / (2*h)

    return gx, gy


x, y = -1.7, 0.8

analytic = gradient(x, y)
numeric = numerical_gradient(x, y)

print("解析梯度:", analytic)
print("数值梯度:", numeric)

assert abs(
    analytic[0] - numeric[0]
) < 1e-7

assert abs(
    analytic[1] - numeric[1]
) < 1e-7

print("梯度检查通过")