Rust图像处理第21节-最小二乘回归：用矩阵求逆解"拟合"问题📌 写在前面 PCA 找"数据的主方向"（无监督）——

🦀 Rust + WASM 实战系列第 21 篇阅读时间：约 5 分钟 | 实战可运行

📌 写在前面

PCA 找"数据的主方向"（无监督）——这一篇做"监督拟合"：给一些 (x, y) 数据点，找一条直线最佳逼近。

公式超简洁：

\beta = (X^T X)^{-1} X^T y

一行矩阵公式搞定"线性回归"——但这行公式是怎么来的？§ 一从最小二乘法的直觉出发，一步一步把它推导出来。代码侧用 nalgebra 库调用 try_inverse() 求逆就行。

🚀 TL;DR

两个真实图像应用：

应用	场景	公式
图像去噪	扫描旧照片有"上面亮、下面暗"的条纹	拟合 $y = a \cdot \text{row} + b$ ，逐行减去
颜色校准	显示器偏色（拍标准色卡拟合修正）	拟合 $\text{actual} = m \cdot \text{target} + b$ ，反向应用

核心算法就是 β = (X^T X)⁻¹ X^T y——矩阵求逆。

一、核心公式：从最小二乘到正规方程

这一节完整推导 $\beta = (X^T X)^{-1} X^T y$ ——从"最小二乘是什么"到"矩阵形式"，一步一步。

1.1 最小二乘是什么？——拟合的直觉

问题：给 $n$ 个数据点 $(x_1, y_1), (x_2, y_2), \dots, (x_n, y_n)$ ，找一条直线 $y = a x + b$ 最佳逼近。

为什么不要求直线穿过所有点？因为数据有噪点——同一 $x$ 对应的 $y$ 不会完全相同（拍照亮度会浮动、传感器有误差、人为测量也会偏）。强行穿过每个点，反而会被噪点带偏。

什么叫"最佳"？——需要一个量化标准来衡量"偏离程度"。

最小二乘的标准：让所有点到直线的纵向距离的平方和最小。

为什么是"纵向距离"（沿 $y$ 方向）而不是"垂直距离"？因为我们要预测 $y$ ——给一个 $x$ ，希望预测的 $\hat y$ 和真实 $y$ 接近。垂直距离在数学上更对称，但不直接对应预测误差。

1.2 优化目标：残差平方和

对每个点定义残差（真实值减预测值）：

e_i = y_i - (a x_i + b)

优化目标：

\min_{a, b} S(a, b) = \sum_{i=1}^{n} e_i^2 = \sum_{i=1}^{n} (y_i - a x_i - b)^2

为什么用"平方"而不是"绝对值"或"立方"？

可微：平方处处可微，绝对值在 $e_i = 0$ 处不可微 → 平方更容易求极值
放大离群点：平方让远离直线的点被"惩罚得更狠"——直线倾向于贴近主体数据，不被极端点带偏
闭式解：最小化平方和 → 线性方程组（直接有解析解，不用迭代优化）

1.3 怎么求最小值？——求导

$S(a, b)$ 是关于 $a$ 和 $b$ 的二元连续可微函数。极值点 ⟺ 偏导数都为 0（微积分基本定理）：

\frac{\partial S}{\partial a} = 0, \quad \frac{\partial S}{\partial b} = 0

只要解出这两个方程，就能找到极值点。这里是极小值而非极大值，因为 $S(a, b)$ 是 $a$ 、 $b$ 的凸二次函数（开口向上的抛物面）——凸函数的局部极小 = 全局最小。

1.4 手算：对 a 和 b 求偏导

把 $S$ 展开，链式法则：

\frac{\partial S}{\partial a} = \sum_i 2(y_i - a x_i - b) \cdot (-x_i) = -2 \sum_i x_i (y_i - a x_i - b)

\frac{\partial S}{\partial b} = \sum_i 2(y_i - a x_i - b) \cdot (-1) = -2 \sum_i (y_i - a x_i - b)

令两个偏导数都等于 0（除以 $-2$ ）：

\sum_i x_i (y_i - a x_i - b) = 0 \quad \cdots (1)

\sum_i (y_i - a x_i - b) = 0 \quad \cdots (2)

展开（把 $a$ 和 $b$ 移到右边）：

\sum_i x_i y_i = a \sum_i x_i^2 + b \sum_i x_i \quad \cdots (1')

\sum_i y_i = a \sum_i x_i + b \cdot n \quad \cdots (2')

两个方程、两个未知数 $a$ 、 $b$ ——理论上能解。但还有 5 个累加和散落在公式里，看着乱。下面把它们打包成矩阵。

1.5 写成矩阵形式——正规方程

把 (1')(2') 写成矩阵：

\begin{bmatrix} \sum x_i^2 & \sum x_i \\ \sum x_i & n \end{bmatrix} \begin{bmatrix} a \\ b \end{bmatrix} = \begin{bmatrix} \sum x_i y_i \\ \sum y_i \end{bmatrix}

关键观察：左边矩阵的元素恰好等于一个 $X^T X$ ，右边向量恰好是 $X^T y$ ——其中 设计矩阵：

X = \begin{bmatrix} x_1 & 1 \\ x_2 & 1 \\ \vdots & \vdots \\ x_n & 1 \end{bmatrix}, \quad y = \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix}

验证 $X^T X$ （一个 $2 \times n$ 矩阵乘 $n \times 2$ 矩阵）：

X^T X = \begin{bmatrix} x_1 & x_2 & \cdots & x_n \\ 1 & 1 & \cdots & 1 \end{bmatrix} \begin{bmatrix} x_1 & 1 \\ x_2 & 1 \\ \vdots & \vdots \\ x_n & 1 \end{bmatrix} = \begin{bmatrix} \sum x_i^2 & \sum x_i \\ \sum x_i & n \end{bmatrix} \checkmark

验证 $X^T y$ （ $2 \times n$ 矩阵乘 $n \times 1$ 向量）：

X^T y = \begin{bmatrix} x_1 & x_2 & \cdots & x_n \\ 1 & 1 & \cdots & 1 \end{bmatrix} \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix} = \begin{bmatrix} \sum x_i y_i \\ \sum y_i \end{bmatrix} \checkmark

所以方程 (1')(2') 恰好就是：

\boxed{X^T X \begin{bmatrix} a \\ b \end{bmatrix} = X^T y}

这叫正规方程。两边左乘 $(X^T X)^{-1}$ （前提： $X^T X$ 可逆，见 § 六踩坑）：

\boxed{\beta = (X^T X)^{-1} X^T y, \quad \text{其中 } \beta = \begin{bmatrix} a \\ b \end{bmatrix}}

这就是 § 写在前面那行"超简洁"的公式——只是 4 个矩阵乘法的打包写法。

1.6 设计矩阵的几何意义

为什么 $X$ 要加一列 $1$ ？

$X$ 的第 $i$ 行 $[x_i, 1]$ 乘以 $\beta = [a, b]^T$ 得到预测值 $\hat y_i = a x_i + b$
这一列 $1$ 就是截距项——没它就只能拟合过原点的直线（ $y = ax$ ）

这个模式可以推广——只要模型是"参数的线性组合"就能用正规方程：

模型	设计矩阵 X 的行
$y = ax + b$	$[x_i, 1]$
$y = ax^2 + bx + c$	$[x_i^2, x_i, 1]$
$y = a_1 x_1 + a_2 x_2 + b$	$[x_{i,1}, x_{i,2}, 1]$
$k$ 阶多项式	$[x_i^k, x_i^{k-1}, \dots, x_i, 1]$

正规方程对所有"线性模型"都成立——这就是为什么它在机器学习里是"基础方程"。

二、Demo 1：图像去噪（行均值拟合）

应用场景

扫描的旧照片经常有这种系统性行噪点：

上面亮、下面暗（扫描灯光不均）
或某些行特别亮（扫描仪传感器问题）

算法

计算每行像素的平均亮度 → 得到 $n$ 个点 $(row_i, \text{avg\_brightness}_i)$
最小二乘拟合： $\text{avg\_brightness} = a \cdot \text{row} + b$
从原图每行减去拟合值 → 去除系统性行噪点

原图：                    处理后：
┌──────────────┐          ┌──────────────┐
│ ░░░ 亮       │          │              │
│ ░ 中等       │     →    │   一致亮度    │
│   ▓ 暗       │          │              │
└──────────────┘          └──────────────┘
   (上面亮/下面暗)            (统一了)

三、Demo 2：颜色校准（拍色卡拟合）

应用场景

显示器偏色（偏暖/偏冷）——用标准色卡（已知 RGB 值）拍照，得到实际显示的 RGB，拟合修正系数。

算法

准备 $n$ 组数据点 $(\text{target}_i, \text{actual}_i)$ ：每组是 3 维 RGB
对每通道独立拟合： $\text{actual} = m \cdot \text{target} + b$
对原图应用反向变换： $\text{corrected} = (\text{actual} - b) / m$

目标红色 (255, 0, 0)  →  实际显示 (260, 10, -5)
                                    ↓
                            拟合 m=1.02, b=10
                                    ↓
反向：corrected_r = (260 - 10) / 1.02 = 245

四、关键代码（用库调用）

Demo 1：行均值去噪

#[wasm_bindgen]
pub fn denoise_row_mean(pixels: &[u8], width: u32, height: u32) -> Vec<u8> {
    let w = width as usize;
    let h = height as usize;

    // 1. 每行平均亮度
    let mut row_means: Vec<f64> = Vec::with_capacity(h);
    for y in 0..h {
        let mut sum = 0.0;
        for x in 0..w {
            let idx = (y * w + x) * 4;
            sum += (pixels[idx] as f64 + pixels[idx + 1] as f64 + pixels[idx + 2] as f64) / 3.0;
        }
        row_means.push(sum / w as f64);
    }

    // 2. 最小二乘拟合 y = ax + b（手算公式）
    let n = h as f64;
    let x_mean = (h as f64 - 1.0) / 2.0;
    let y_mean = row_means.iter().sum::<f64>() / n;
    let mut num = 0.0;
    let mut den = 0.0;
    for (i, &y) in row_means.iter().enumerate() {
        num += (i as f64 - x_mean) * (y - y_mean);
        den += (i as f64 - x_mean).powi(2);
    }
    let a = if den != 0.0 { num / den } else { 0.0 };
    let b = y_mean - a * x_mean;

    // 3. 逐行减去拟合值
    let mut out = pixels.to_vec();
    for y in 0..h {
        let offset = a * y as f64 + b;
        for x in 0..w {
            let idx = (y * w + x) * 4;
            out[idx]     = (pixels[idx] as f64     - offset).clamp(0.0, 255.0) as u8;
            out[idx + 1] = (pixels[idx + 1] as f64 - offset).clamp(0.0, 255.0) as u8;
            out[idx + 2] = (pixels[idx + 2] as f64 - offset).clamp(0.0, 255.0) as u8;
        }
    }
    out
}

Demo 2：颜色校准（用库）

/// 单通道最小二乘：actual = m * target + b
fn solve_channel(target: &[f64], actual: &[f64]) -> (f64, f64) {
    let n = target.len();
    // 构造设计矩阵 X = [target, 1]（n × 2）
    let mut x_mat = DMatrix::<f64>::zeros(n, 2);
    for i in 0..n {
        x_mat[(i, 0)] = target[i];
        x_mat[(i, 1)] = 1.0;
    }
    let y_vec = DVector::from_vec(actual.to_vec());

    // 正规方程：β = (X^T X)^(-1) X^T y
    let xtx = x_mat.transpose() * &x_mat;
    let xty = x_mat.transpose() * y_vec;
    let beta = xtx.try_inverse().unwrap() * xty;

    (beta[0], beta[1])  // (m, b)
}

核心就 1 行库调用：β = (X^T X)⁻¹ X^T y——xtx.try_inverse() * xty。

五、前端效果展示

六、踩坑提醒

1. 协方差矩阵奇异时 `try_inverse` 返回 None

// ❌ 直接解引用会 panic
let inv = xtx.try_inverse().unwrap();

// ✅ 用 unwrap_or 兜底
let inv = xtx.try_inverse().unwrap_or(DMatrix::identity(2, 2));

2. 校准点 < 2 组不能拟合

if n < 2 {
    return pixels.to_vec();  // 至少 2 组才能确定一条直线
}

3. 校准的反向变换有除零风险

if mr.abs() > 1e-6 {
    out[idx] = ((r - br) / mr).clamp(0.0, 255.0) as u8;
} else {
    out[idx] = r as u8;  // m 接近 0：直接复制
}

4. 行均值去噪对整行偏色有用，对局部噪点没用

✅ 适用：扫描灯光不均、传感器行漂移
❌ 不适用：椒盐噪点、随机噪点（这些应该用前面任务 25 的去噪）

七、接下来

任务 22：深挖 PCA：从最大化方差到最小化重建误差 从优化问题本身推导 PCA 的数学原理。

关键 insight："最大化方差" ≡ "最小化重建误差" ≡ 协方差矩阵的特征向量——三个看似不同的视角其实指向同一个数学对象。

这一篇的"求偏导 = 0 → 联立方程组 → 矩阵形式"推导套路和§ 1.4-1.5一脉相承——PCA 和最小二乘都是"凸优化 → 求偏导 → 矩阵方程"。

一句话总结

最小二乘 = 找一条直线最佳拟合数据点—— $\beta = (X^T X)^{-1} X^T y$ 。

代码 1 行库调用：xtx.try_inverse() * xty。

两个真实应用：图像去噪（行噪点）+ 颜色校准（显示器偏色）。

📦 项目地址：pixel-math-wasm 🦀 Rust + WebAssembly 实战系列

🏷️ 标签：#Rust #WebAssembly #图像处理 #最小二乘 #回归 #矩阵求逆 #nalgebra

Rust图像处理第21节-最小二乘回归：用矩阵求逆解"拟合"问题

📌 写在前面

🚀 TL;DR

📖 目录

一、核心公式：从最小二乘到正规方程

1.1 最小二乘是什么？——拟合的直觉

1.2 优化目标：残差平方和

1.3 怎么求最小值？——求导

1.4 手算：对 a 和 b 求偏导

1.5 写成矩阵形式——正规方程

1.6 设计矩阵的几何意义

二、Demo 1：图像去噪（行均值拟合）

应用场景

算法

三、Demo 2：颜色校准（拍色卡拟合）

应用场景

算法

四、关键代码（用库调用）

Demo 1：行均值去噪

Demo 2：颜色校准（用库）

五、前端效果展示

六、踩坑提醒

1. 协方差矩阵奇异时 try_inverse 返回 None

2. 校准点 < 2 组不能拟合

3. 校准的反向变换有除零风险

4. 行均值去噪对整行偏色有用，对局部噪点没用

七、接下来

一句话总结

1. 协方差矩阵奇异时 `try_inverse` 返回 None