PVNet 关键点投票异常：从求解器到最小二乘条件数

学习教程更新于 2026-10-03编辑本文

1. 保留的故障记录与证据边界

原记录发表于 2024-08-08：在 NVIDIA Orin Nano 上，将基于 clean-pvnet 的自定义数据集流程改为 TensorRT/Python 后，mask 看起来正常，但由向量场投票得到的关键点严重偏离预期，最终位姿不准确。当时预计坐标为数百像素，日志却出现约两万的值，例如 (20114.7383, -663.1294)。

原始 GPU 张量日志：多个关键点坐标达到一万至两万量级，第一点为 20114.7383、-663.1294

当时记录的处理是“将最后一步求解改为 QR 后恢复”。这个现象值得保留，但原稿没有保存失败样本的矩阵、修改前后代码、库版本、精度和评估结果，无法据此证明 LU 本身不适用于这类问题，也不能断言换 QR 就能修复所有异常。本文补充数学模型与可运行的局部检查；没有重新运行那份 TensorRT 引擎或原始数据集。

2. 投票后的精化为什么是最小二乘

对一个关键点，设前景像素位置为 $p_i=(x_i,y_i)^\mathsf{T}$ ，预测单位方向为 $d_i=(d_{ix},d_{iy})^\mathsf{T}$ 。关键点 $k$ 应落在从像素沿该方向延伸的直线上。取垂直法向量：

n_i=\begin{bmatrix}d_{iy}\\-d_{ix}\end{bmatrix}, \qquad n_i^\mathsf{T}(k-p_i)=0.

RANSAC 先找候选交点及其内点；只用选中的内点建立：

A=\begin{bmatrix}n_1^\mathsf{T}\\n_2^\mathsf{T}\\\vdots\\n_m^\mathsf{T}\end{bmatrix}, \qquad b=\begin{bmatrix}n_1^\mathsf{T}p_1\\n_2^\mathsf{T}p_2\\\vdots\\n_m^\mathsf{T}p_m\end{bmatrix}, \qquad \hat{k}=\underset{k}{\operatorname{argmin}}\,\lVert Ak-b\rVert_2^2.

实际方向含噪声，直线通常不会精确相交。若有置信权重 $w_i\geq0$ ，应把第 $i$ 行的 $A$ 和 $b$ 同时乘以 $\sqrt{w_i}$ ；直接乘以 $w_i$ 会改变预期权重。0/1 内点筛选是特殊情况。

本次检查 clean-pvnet 的 ransac_voting_gpu.py：代码把像素索引从行列顺序换成 (x,y)，对内点法向量构造 ATA 与 ATb；不同函数分支使用矩阵求逆或旧 torch.solve 辅助函数。它说明实际移植需要逐步对齐什么，并不说明历史故障使用了完全相同的提交。应固定仓库提交再比较，不能只记录“基于 clean-pvnet”。上游实现

条件数比求解器名字更能解释失真

法向量方向太接近时， $A$ 的两列接近线性相关；很多几乎平行的投票线也不能稳定确定二维交点。对于满列秩的实矩阵，以 2 范数计算：

\kappa_2(A)=\frac{\sigma_{\max}(A)}{\sigma_{\min}(A)}, \qquad \kappa_2(A^\mathsf{T}A)=\kappa_2(A)^2.

因此先形成正规方程 $A^\mathsf{T}Ak=A^\mathsf{T}b$ 会放大病态性，并在有限精度下损失有效信息。对**原始 $A$ **做 QR/SVD，可避免额外构造正规方程；对已经形成的 ATA 再做 QR，则不能恢复此前舍入丢失的信息。几何退化、错误坐标和外点也不会因为更换分解方法而消失。

3. LU 与 QR 各自在做什么

LU：方程求解与主元

无主元交换、 $L$ 对角为 1 的 $3\times3$ 分解写作：

A=LU, \quad L=\begin{bmatrix}1&0&0\\l_{21}&1&0\\l_{31}&l_{32}&1\end{bmatrix}, \quad U=\begin{bmatrix}u_{11}&u_{12}&u_{13}\\0&u_{22}&u_{23}\\0&0&u_{33}\end{bmatrix}.

从 $A=LU$ 逐项比较可得：

\begin{aligned} u_{11}&=a_{11},&u_{12}&=a_{12},&u_{13}&=a_{13},\\ l_{21}&=a_{21}/u_{11},&l_{31}&=a_{31}/u_{11},\\ u_{22}&=a_{22}-l_{21}u_{12},&u_{23}&=a_{23}-l_{21}u_{13},\\ l_{32}&=(a_{32}-l_{31}u_{12})/u_{22},\\ u_{33}&=a_{33}-l_{31}u_{13}-l_{32}u_{23}. \end{aligned}

这些公式要求相关主元非零。非奇异方阵的各阶顺序主子式非零，是这种无交换分解存在的条件；数值实现通常采用部分主元交换，写作 $PA=LU$ ，有些库使用等价但排列矩阵定义相反的 $A=PLU$ 。必须以接口文档为准。LU 分解本身也可推广到矩形矩阵；“只适用于方阵”过于绝对，常见方阵线性方程求解才是这里的比较对象。SciPy LU 约定

对于良态、满秩的方阵系统，带主元的 LU 可以得到正确结果。这里更应避免显式求逆后再乘右端项，并先确认求解器究竟接收 A，还是 ATA。

QR：列正交与最小二乘

当 $A\in\mathbb{R}^{m\times n}$ 且 $m\geq n$ ，简约 QR 为：

A=QR,\quad Q\in\mathbb{R}^{m\times n},\quad Q^\mathsf{T}Q=I_n, \quad R\in\mathbb{R}^{n\times n}.

这里 $Q$ 是列正交归一矩阵；矩形 $Q$ 一般不满足 $QQ^\mathsf{T}=I_m$ ，不应不加区分地称为方形正交矩阵。满列秩时 $R$ 非奇异，最小二乘解通过三角方程求出：

R\hat{k}=Q^\mathsf{T}b.

用 Gram–Schmidt 理解分解过程时，令 $a_k$ 为 $A$ 的第 $k$ 列：

\begin{aligned} v_1&=a_1,&r_{11}&=\lVert v_1\rVert_2,&q_1&=v_1/r_{11},\\ r_{jk}&=q_j^\mathsf{T}a_k\quad(j<k),\\ v_k&=a_k-\sum_{j=1}^{k-1}q_jr_{jk},&r_{kk}&=\lVert v_k\rVert_2,&q_k&=v_k/r_{kk}. \end{aligned}

投影也可写成 $\operatorname{proj}_{q_j}(a_k)=(q_j^\mathsf{T}a_k)q_j$ ，因为 $q_j$ 已归一化。这修复了原稿把下标写成星号的问题，并用 $v_k$ 区分归一化之前的向量。若 $r_{kk}=0$ ，不能继续相除；经典 Gram–Schmidt 在浮点病态问题中也可能丢失正交性。工程上使用库的 Householder QR、带列主元的 QR 或 SVD，并根据秩选择处理策略。NumPy QR、LAPACK 最小二乘驱动

4. 可运行的 CPU 参考：只验证内点精化

下面接收一个关键点、已经筛选好的内点，统一为 (x,y) 像素坐标和朝向关键点的方向。它不实现 mask、网络推理、RANSAC 或 PnP。使用 float64 建立调试基准；实际 GPU 流程是否能用相同精度、吞吐和 API 需另行评估。

import numpy as np

def refine_keypoint(points, directions, max_condition=1e8):
    points = np.asarray(points, dtype=np.float64)
    directions = np.asarray(directions, dtype=np.float64)
    if points.ndim != 2 or points.shape[1] != 2:
        raise ValueError("points must have shape (N, 2)")
    if directions.shape != points.shape or len(points) < 2:
        raise ValueError("need matching directions and at least two points")
    if not np.isfinite(points).all() or not np.isfinite(directions).all():
        raise ValueError("non-finite input")
    if not np.isfinite(max_condition) or max_condition < 1:
        raise ValueError("invalid condition limit")
    lengths = np.linalg.norm(directions, axis=1)
    if np.any(lengths == 0) or not np.isfinite(lengths).all():
        raise ValueError("invalid direction length")
    directions = directions / lengths[:, None]
    A = np.column_stack((directions[:, 1], -directions[:, 0]))
    b = np.sum(A * points, axis=1)
    if not np.isfinite(b).all():
        raise ValueError("non-finite linear-system data")
    keypoint, _, rank, singular = np.linalg.lstsq(A, b, rcond=None)
    if not np.isfinite(keypoint).all():
        raise ValueError("non-finite solution")
    if rank < 2:
        raise ValueError("parallel voting lines: keypoint is not unique")
    condition = singular[0] / singular[-1]
    if condition > max_condition:
        raise ValueError("ill-conditioned voting geometry")
    residual = A @ keypoint - b
    return keypoint, {
        "rank": int(rank),
        "condition": float(condition),
        "rms_distance": float(np.sqrt(np.mean(residual ** 2))),
    }

points = np.array([[100., 100.], [500., 100.],
                   [100., 400.], [500., 400.]])
expected = np.array([320., 240.])
directions = expected - points
keypoint, diagnostics = refine_keypoint(points, directions)
assert np.allclose(keypoint, expected, atol=1e-9)
print(keypoint, diagnostics)

max_condition=1e8 是这个 float64 教学基准的示例门限，不是网络精度或像素误差的保证。应依据允许误差、输入噪声和实际数值精度设置条件数、内点数量及残差门限。这里的残差是归一化法向量下的点到直线距离；它小也不保证关键点准确，因为近乎平行的线可能在很远处相交。

lstsq 返回的 residuals 在某些形状或秩不足时是空数组，所以例子主动计算残差。秩不足时，库仍可能给出一个最小范数解；这不等于数据已经唯一确定二维关键点，本例明确拒绝这种情况。NumPy lstsq

5. 在真实移植中逐层比较

固定一张输入图、一份权重和随机种子，保存模型提交、TensorRT/CUDA/PyTorch 版本、输入尺寸与精度模式。
比较原框架和引擎的 mask 及向量场，确认通道顺序、(x,y) 与 (row,column)、缩放/裁剪/填充逆变换。mask 正常不能证明向量头正确。
保存同一关键点的 points、directions、内点索引、原始 A,b 和求解结果。检查有限值、方向长度和批次维度，不能在 reshape/transpose 时混淆关键点。
用相同内点分别比较 float64 CPU 参考与生产求解路径，记录秩、最小奇异值、条件数、残差和关键点误差。再逐项比较 FP32/FP16，避免同时更改数据和求解器。
若几何退化，返回明确的失败/置信状态，并按应用规则保留 RANSAC 候选或跳过位姿更新；不要静默用零点或单位矩阵假装求解成功。关键点允许落在图像外时，也不能简单裁剪坐标掩盖问题。
关键点恢复后，继续检查相机内参、图像缩放和 PnP 的 2D/3D 点对应关系，用重投影误差确认位姿。

本页补充验证运行于 NumPy 1.26.4：良态四条投票线恢复 (320,240)；退化平行线、近乎平行线和非法输入被拒绝；同一良态问题的 QR 与正规方程求解相符；另外的人工病态矩阵展示了 FP32 正规方程可能丢失秩。它们验证数学说明与参考代码，不代替原故障复现。