线性回归与梯度下降

# 两条样本 (x=0, y=1)、(x=1, y=3)，损失 = Σ(ŷ−y)²/2，w=b=0
SAMPLES = [(0, 1), (1, 3)]


def grad(w, b):
    dw = 0.0
    db = 0.0
    for x, y in SAMPLES:
        err = b + w * x - y          # ŷ − y
        dw += err * x                # 对 w：误差乘特征值
        db += err                    # 对 b：误差本身
    return dw, db


def loss(w, b):
    return sum((b + w * x - y) ** 2 for x, y in SAMPLES) / 2


dw, db = grad(0, 0)
assert dw == -3 and db == -4
assert abs((loss(0, 0.001) - loss(0, 0)) / 0.001 - db) < 0.01     # 数值梯度 −3.999
assert abs((loss(0.001, 0) - loss(0, 0)) / 0.001 - dw) < 0.01     # 数值梯度 −2.9995
# 沿负梯度走一步（eta=0.1）：w=0.3、b=0.4，损失 5.0 → 2.825（浮点乘法有尾差，用容差比较）
w1, b1 = 0 - 0.1 * dw, 0 - 0.1 * db
assert abs(w1 - 0.3) < 1e-12 and abs(b1 - 0.4) < 1e-12
assert abs(loss(w1, b1) - 2.825) < 1e-12

断言覆盖第 04 节第一张表（−3、−4）、两个数值梯度（−2.9995、−3.999）和走一步后的 w=0.3、b=0.4、损失 2.825。

import sys


def fmt(v):
    if abs(v) < 0.00005:            # 避免打出 -0.0000
        v = 0.0
    return f"{v:.4f}"


def main():
    data = sys.stdin.read().split()
    if not data:
        return
    pos = 0
    n, d, m, epochs = (int(data[i]) for i in range(4))
    pos = 4
    eta = float(data[pos]); pos += 1
    xs, ys = [], []
    for _ in range(n):
        xs.append([float(v) for v in data[pos:pos + d]]); pos += d
        ys.append(float(data[pos])); pos += 1
    queries = []
    for _ in range(m):
        queries.append([float(v) for v in data[pos:pos + d]]); pos += d

    # fit：只用训练集算每列最小值与最大值
    mins = [min(row[j] for row in xs) for j in range(d)]
    maxs = [max(row[j] for row in xs) for j in range(d)]

    def transform(row):                 # 常数列归 0；超出范围不截断
        return [0.0 if maxs[j] == mins[j] else (row[j] - mins[j]) / (maxs[j] - mins[j])
                for j in range(d)]

    train = [transform(row) for row in xs]   # 训练集只归一化一次

    w = [0.0] * d
    b = 0.0
    for _ in range(epochs):
        gw = [0.0] * d
        gb = 0.0
        for row, y in zip(train, ys):        # 先把全部样本的梯度累加完
            err = b + sum(w[j] * row[j] for j in range(d)) - y
            gb += err
            for j in range(d):
                gw[j] += err * row[j]
        scale = 2.0 / n                      # 均方误差的求导系数
        b -= eta * scale * gb                # 再一起更新（同步更新）
        for j in range(d):
            w[j] -= eta * scale * gw[j]

    out = []
    for q in queries:
        z = transform(q)
        out.append(fmt(b + sum(w[j] * z[j] for j in range(d))))
    print("\n".join(out))


main()

transform 同时服务训练集和查询：常数列归 0、超出范围不截断都在这一个函数里。用题面示例 → 0.4000 / 0.7000 和第 05 节的三条规则用例核对。

import sys


def solve_linear(a, rhs):
    """高斯消元（列主元）解 a·x = rhs，返回 x。"""
    size = len(rhs)
    for col in range(size):
        pivot = max(range(col, size), key=lambda r: abs(a[r][col]))
        a[col], a[pivot] = a[pivot], a[col]
        rhs[col], rhs[pivot] = rhs[pivot], rhs[col]
        div = a[col][col]
        for j in range(col, size):
            a[col][j] /= div
        rhs[col] /= div
        for r in range(size):
            if r == col or a[r][col] == 0:
                continue
            factor = a[r][col]
            for j in range(col, size):
                a[r][j] -= factor * a[col][j]
            rhs[r] -= factor * rhs[col]
    return rhs


def fmt(v):
    if abs(v) < 0.00005:
        v = 0.0
    return f"{v:.4f}"


def main():
    data = sys.stdin.read().split()
    if not data:
        return
    n, d, m = int(data[0]), int(data[1]), int(data[2])
    lam = float(data[3])
    pos = 4
    size = d + 1                        # 第 0 列是恒为 1 的偏置列
    a = [[0.0] * size for _ in range(size)]
    rhs = [0.0] * size
    for _ in range(n):
        z = [1.0] + [float(v) for v in data[pos:pos + d]]
        pos += d
        y = float(data[pos]); pos += 1
        for r in range(size):           # 累加 ZᵀZ 与 Zᵀy
            rhs[r] += z[r] * y
            for c in range(size):
                a[r][c] += z[r] * z[c]
    for j in range(1, size):            # λ 只加在权重的对角元上，偏置不加
        a[j][j] += lam
    w = solve_linear(a, rhs)            # w[0] 是偏置，w[1..d] 是各列权重

    out = []
    for _ in range(m):
        row = []
        missing = -1
        for j in range(d):
            token = data[pos]; pos += 1
            if token == "?":
                missing = j
                row.append(0.0)
            else:
                row.append(float(token))
        target = float(data[pos]); pos += 1
        coef = w[missing + 1]
        if abs(coef) < 1e-10:
            out.append("NA")
            continue
        known = w[0] + sum(w[j + 1] * row[j] for j in range(d) if j != missing)
        out.append(fmt((target - known) / coef))
    print("\n".join(out))


main()

solve_linear 是列主元高斯消元，返回的 w[0] 是偏置。用题面两组示例（2.0030 / 2.9917、NA）核对。

import math
import sys


def sigmoid(z):
    if z >= 0:                       # 按符号分支：exp 的参数始终非正，不会溢出
        return 1.0 / (1.0 + math.exp(-z))
    e = math.exp(z)
    return e / (1.0 + e)


def fmt(v):
    if abs(v) < 0.00005:
        v = 0.0
    return f"{v:.4f}"


def main():
    data = sys.stdin.read().split()
    if not data:
        return
    n, d, m, epochs = (int(data[i]) for i in range(4))
    eta = float(data[4])
    pos = 5
    xs, ys = [], []
    for _ in range(n):
        xs.append([float(v) for v in data[pos:pos + d]]); pos += d
        ys.append(float(data[pos])); pos += 1
    queries = []
    for _ in range(m):
        queries.append([float(v) for v in data[pos:pos + d]]); pos += d

    mins = [min(row[j] for row in xs) for j in range(d)]
    maxs = [max(row[j] for row in xs) for j in range(d)]

    def transform(row):
        return [0.0 if maxs[j] == mins[j] else (row[j] - mins[j]) / (maxs[j] - mins[j])
                for j in range(d)]

    train = [transform(row) for row in xs]
    w = [0.0] * d
    b = 0.0
    for _ in range(epochs):
        gw = [0.0] * d
        gb = 0.0
        for row, y in zip(train, ys):
            err = sigmoid(b + sum(w[j] * row[j] for j in range(d))) - y   # p − y
            gb += err
            for j in range(d):
                gw[j] += err * row[j]
        b -= eta * gb / n                # 平均梯度，同步更新
        for j in range(d):
            w[j] -= eta * gw[j] / n

    out = [" ".join(fmt(v) for v in [b] + w)]
    for q in queries:
        z = transform(q)
        out.append(fmt(sigmoid(b + sum(w[j] * z[j] for j in range(d)))))
    print("\n".join(out))


main()

与 AI002 只差三处：预测值套 sigmoid、梯度乘 1/n 而不是 2/n、第一行先输出 b 与 w。sigmoid 按符号分支——|z| 达到 10⁵ 时直接写 1/(1+exp(-z)) 会溢出（第 09 节错误表）。

import math
import sys


def fmt(v, digits):
    if abs(v) < 0.5 * 10 ** (-digits):
        v = 0.0
    return f"{v:.{digits}f}"


def main():
    data = sys.stdin.read().split()
    if not data:
        return
    r, c = int(data[0]), int(data[1])
    pos = 2
    rows = []
    for _ in range(r):
        rows.append([float(v) for v in data[pos:pos + c]]); pos += c
    labels = [int(v) for v in data[pos:pos + r]]

    out = []
    total = 0.0
    for row, label in zip(rows, labels):
        mx = max(row)
        exps = [math.exp(z - mx) for z in row]     # 减最大值：最大项 exp(0)=1，不会上溢
        s = sum(exps)
        out.append(" ".join(fmt(e / s, 4) for e in exps))
        total += math.log(s) - (row[label] - mx)    # −ln p = ln Σexp(z−max) − (z_label − max)
    out.append(fmt(total / r, 6))
    print("\n".join(out))


main()

交叉熵按 ln Σexp(z−max) − (z_label − max) 计算：1 3 / 1000 0 -1000 / 2 真实类别的概率下溢成 0，对它取对数会报数学域错误，正确答案是 2000.000000。

import sys


def fmt(v):
    if abs(v) < 0.00005:
        v = 0.0
    return f"{v:.4f}"


def main():
    data = sys.stdin.read().split()
    if not data:
        return
    d, h, m = int(data[0]), int(data[1]), int(data[2])
    pos = 3
    hidden = []                                   # 每项：(d 个权重, 偏置)
    for _ in range(h):
        ws = [float(v) for v in data[pos:pos + d]]
        pos += d
        hidden.append((ws, float(data[pos]))); pos += 1
    vs = [float(v) for v in data[pos:pos + h]]; pos += h
    c = float(data[pos]); pos += 1

    out = []
    for _ in range(m):
        x = [float(v) for v in data[pos:pos + d]]; pos += d
        y = c
        for (ws, bias), v in zip(hidden, vs):     # 隐藏神经元按下标从小到大
            z = bias
            for wj, xj in zip(ws, x):             # 求和按下标从小到大
                z += wj * xj
            a = z if z > 0 else 0.0               # ReLU
            y += v * a
        out.append(fmt(y))
    print("\n".join(out))


main()

两层循环都按下标从小到大求和，与题面一致。用第 10 节练习 5 的题面示例（1.5000 / 0.5000）核对。

# 两条样本 (x=0, y=1)、(x=1, y=3)，损失 = Σ(ŷ−y)²/2，w=b=0
SAMPLES = [(0, 1), (1, 3)]


def grad(w, b):
    dw = 0.0
    db = 0.0
    for x, y in SAMPLES:
        err = b + w * x - y          # ŷ − y
        dw += err * x                # 对 w：误差乘特征值
        db += err                    # 对 b：误差本身
    return dw, db


def loss(w, b):
    return sum((b + w * x - y) ** 2 for x, y in SAMPLES) / 2


dw, db = grad(0, 0)
assert dw == -3 and db == -4
assert abs((loss(0, 0.001) - loss(0, 0)) / 0.001 - db) < 0.01     # 数值梯度 −3.999
assert abs((loss(0.001, 0) - loss(0, 0)) / 0.001 - dw) < 0.01     # 数值梯度 −2.9995
# 沿负梯度走一步（eta=0.1）：w=0.3、b=0.4，损失 5.0 → 2.825（浮点乘法有尾差，用容差比较）
w1, b1 = 0 - 0.1 * dw, 0 - 0.1 * db
assert abs(w1 - 0.3) < 1e-12 and abs(b1 - 0.4) < 1e-12
assert abs(loss(w1, b1) - 2.825) < 1e-12

见第 08 节展开区（同一份代码）。