Part I · 第 1 章线性代数 · 1.12

奇异值分解 SVD

任何矩阵，不管看起来多乱，都只做三件事：先转一下，再沿坐标轴拉伸，最后再转一下。所以单位圆永远会变成一个椭圆，椭圆的半轴长就是奇异值。从奇异值能读出机械臂往哪个方向动得快、往哪个方向动不了，也能读出一个问题离“无解”还有多远。

约 35 分钟前置：1.3 矩阵即变换 · 1.10 正交矩阵与旋转 · 1.11 特征值 7 个交互

1为什么机器人需要它

一台机械臂的每个关节都能转，但末端执行器在不同方向上的“灵活程度”并不一样。关节速度 $\dot{\mathbf q}$ 经过雅可比矩阵 (Jacobian) $J$ 变成末端速度 $\dot{\mathbf x} = J\dot{\mathbf q}$。给所有关节同样的速度预算（一个单位圆），末端能达到的速度就铺成一个椭圆，叫可操作度椭圆 (manipulability ellipse)。椭圆在某个方向很扁，机械臂在那个方向就使不上劲；椭圆压成一条线段，机械臂就进入了奇异位形 (singularity)，在某个方向上彻底动不了。

这个椭圆的长短轴和朝向，由 $J$ 的奇异值分解（SVD）完整描述。同一个工具在机器人软件里随处可见：

最小二乘与伪逆：数值逆运动学、相机和手眼标定，都要解“方程比未知数多（或少）”的线性方程组（1.13 节）。
点云 PCA：用一小片点估计桌面的法向量，或者估计一根杆子的主轴方向。
秩与条件数：我的雅可比离奇异还有多远？这个标定问题的数据够不够？
压缩与降维：用很少的数近似一张图像或一个大矩阵。

2直觉：圆变成椭圆

1.3 节里，我们看矩阵把方格网变成平行四边形网。这一节换一个更好用的“探针”：单位圆。圆上每个点都是一个长度为 1 的方向 $\mathbf x$，经过矩阵 $A$ 变成 $A\mathbf x$。所有 $A\mathbf x$ 合在一起，是什么形状？

答案永远是一个椭圆（极端情况下退化成一条线段，甚至一个点）。椭圆有最长的半轴，也有最短的半轴。于是问题变成：哪一个输入方向被拉得最长？哪一个被压得最短？

在下图左边拖动圆上的点 $\mathbf x$。旁边的曲线图画出 $|A\mathbf x|$ 随方向角的变化，也可以直接在曲线上拖动。

FIG 1拉一拉单位圆

输入 x 与输出 Ax

|Ax| 随输入方向变化

虚线圆是所有输入方向，青色椭圆是它们的像。黄色箭头 $A\mathbf x$ 的长度就是曲线图的高度。注意曲线在一圈里有两个峰、两个谷：$\mathbf x$ 和 $-\mathbf x$ 被拉伸的倍数相同。

找到被拉伸得最多的输入方向（让 $|A\mathbf x|$ 达到最大）。
再找到被拉伸得最少的方向。它和上一个方向是什么关系？

换几个矩阵试试，你会发现两件总是成立的事：

拉伸最多的输入方向 $\cv{\mathbf v_1}$ 和拉伸最少的输入方向 $\cw{\mathbf v_2}$ 互相垂直。
它们的像 $A\cv{\mathbf v_1}$ 和 $A\cw{\mathbf v_2}$ 正好落在椭圆的长轴和短轴上，所以也互相垂直。

核心洞见

每个矩阵都有一对“特别的”正交输入方向 $\cv{\mathbf v_1}, \cw{\mathbf v_2}$，它们被送到一对正交输出方向 $\cv{\mathbf u_1}, \cw{\mathbf u_2}$ 上，长度只乘了 $\sigma_1, \sigma_2$ 倍：$A\mathbf v_i = \sigma_i \mathbf u_i$。数 $\sigma_1 \ge \sigma_2 \ge 0$ 叫奇异值 (singular values)，也就是椭圆的两个半轴长。这就是奇异值分解的全部内容，剩下的只是把它写成矩阵。

3动手玩：转、拉、转

既然 $A$ 把正交的 $\cv{\mathbf v_1}, \cw{\mathbf v_2}$ 送到正交的 $\cv{\mathbf u_1}, \cw{\mathbf u_2}$ 方向，它的动作就可以拆成三步：

$V^\top$：旋转。把 $\cv{\mathbf v_1}, \cw{\mathbf v_2}$ 转到 x 轴和 y 轴上。
$\Sigma$：沿坐标轴拉伸。x 方向乘 $\sigma_1$，y 方向乘 $\sigma_2$。圆在这一步变成椭圆。
$U$：再旋转。把两根坐标轴转到 $\cv{\mathbf u_1}, \cw{\mathbf u_2}$ 方向（如果 $A$ 会让平面翻面，这一步还带一次镜像）。

拖动阶段滑块，或者点按钮，看单位圆一步一步变成 $A$ 的椭圆。圆上的彩色小点和字母 F 帮你追踪每个点的去向；紫色虚线椭圆是终点。

FIG 2转–拉–转实验室

矩阵 A（拖动数字或输入）

预设

黄色箭头从 $\cv{\mathbf v_1}$ 出发，在第 ① 步结束时正好躺在 x 轴上（图中标作 $\mathbf e_1$，即 x 方向的单位向量；$\mathbf e_2$ 是 y 方向），第 ② 步被拉长到 $\sigma_1$，第 ③ 步转到 $\sigma_1\cv{\mathbf u_1}$；紫色箭头 $\cw{\mathbf v_2}$ 同理。红绿箭头是 $\ihat, \jhat$ 被带去的位置，最后正好是 $A$ 的两列。下方公式里，$U$ 的两列和 $V^\top$ 的两行按同样的颜色标出。

试一试

改矩阵，让最后的椭圆变成一个比单位圆大（或小）的正圆。这时第 ② 步在做什么？这类矩阵叫什么？
选一个会让平面翻面的矩阵（$\det A < 0$），把阶段拖过第 ③ 步，看字母 F 怎样被镜像。
找一个非对角矩阵，让第 ① 步和第 ③ 步恰好互相抵消（$U = V$）。提示：让两个非对角元相等，再让对角元足够大。

把三步写成符号，就是本节的主角：

$$ A \;=\; \underbrace{\hl{U}{U}}_{(3)}\;\underbrace{\hl{S}{\Sigma}}_{(2)}\;\underbrace{\hl{Vt}{V^\top}}_{(1)}, \qquad A\,\hl{v1}{\cv{\mathbf v_1}} = \hl{s1}{\sigma_1}\,\hl{u1}{\cv{\mathbf u_1}}, \qquad A\,\hl{v2}{\cw{\mathbf v_2}} = \hl{s2}{\sigma_2}\,\hl{u2}{\cw{\mathbf u_2}} $$

提示：把鼠标悬停在公式中的符号上，FIG 2 会画出对应阶段的形状，或加粗对应的向量。注意矩阵从右往左读：最右边的 $V^\top$ 最先作用。

4奇异值是矩阵的体检报告

奇异值不关心 $A$ 往哪个方向转，只关心它拉伸了多少。几个数就能说明很多事：

$\sigma_1$：$A$ 能把单位向量拉到的最大长度，也叫 $A$ 的 2-范数 $\|A\|_2$。
$\sigma_2$（一般地，最小的 $\sigma_n$）：最弱方向上还剩多少。$\sigma_n = 0$ 说明有一个方向被压成了零，信息丢了，$A$ 不可逆。
秩 (rank) 等于非零奇异值的个数，也就是椭圆“真正有几维”。
条件数 $\kappa = \sigma_1/\sigma_n$：椭圆有多扁。解 $A\mathbf x = \mathbf b$ 时，$\mathbf b$ 里的相对误差最多会被放大 $\kappa$ 倍。

拖动两根列向量 Aî、Aĵ，看右边的柱状图和读数怎么变。

FIG 3秩与条件数

黄色和紫色线段是椭圆的两条半轴 $\sigma_1\cv{\mathbf u_1}$、$\sigma_2\cw{\mathbf u_2}$。当两列平行时，整个平面被压到一条直线上：椭圆塌成线段，$\sigma_2 = 0$，秩降为 1。数值计算里不会有精确的 0，所以“秩”要配一个容差，这里取 0.02。

让矩阵变成秩 1：$\sigma_2 \to 0$，椭圆塌成一条线段。这时条件数是多少？

真实系统里

numpy.linalg.matrix_rank 和 Eigen 的秩判断，内部都是先做 SVD，再数有多少个奇异值大于某个容差。标定程序常把整组奇异值打印出来：如果最小的那个比最大的小十几个数量级，说明某个参数根本没被数据约束住，比如手眼标定时机械臂只绕一个轴转过。

5第二种视角：一层一层叠起来

把 $A = U\Sigma V^\top$ 按列和行展开，同一个分解可以读成一个求和：

$$ A \;=\; \sigma_1\,\cv{\mathbf u_1}\cv{\mathbf v_1}^{\!\top} + \sigma_2\,\cw{\mathbf u_2}\cw{\mathbf v_2}^{\!\top} + \cdots + \sigma_r\,\mathbf u_r\mathbf v_r^{\top}. $$

每一项 $\mathbf u_i\mathbf v_i^\top$ 都是“一列乘一行”，是秩为 1 的矩阵，画出来像一块格子布。奇异值从大到小排好，前面几层承载了大部分内容。只保留前 $k$ 层，就得到 $A$ 的秩 $k$ 近似 $A_k$。Eckart–Young 定理说：在所有秩不超过 $k$ 的矩阵里，$A_k$ 离 $A$ 最近。

一张灰度图就是一个矩阵。下面这张 80×80 的图有 6400 个数，拖动 $k$ 看看要叠多少层才认得出来。

FIG 4低秩图像近似

奇异值谱（对数坐标，可点击）

左：原图。中：前 $k$ 层之和 $A_k$。右：第 $k$ 层 $\sigma_k\mathbf u_k\mathbf v_k^\top$ 本身（灰色为 0，亮为正，暗为负），每一层都是横竖条纹的“格子布”。存储 $A_k$ 只需要 $k$ 个 $\sigma$、$k$ 列 $\mathbf u$ 和 $k$ 行 $\mathbf v^\top$，共 $k(80+80+1)$ 个数。

找到最小的 $k$，使相对误差 $\|A-A_k\|_F/\|A\|_F$ 低于 5%。它需要存多少个数？

同一个分解，两种读法

几何读法：$A = U\Sigma V^\top$ 是“转、拉、转”三个动作的复合，适合理解单个变换（雅可比、坐标变换）。
代数读法：$A = \sum_i \sigma_i\mathbf u_i\mathbf v_i^\top$ 是一摞按重要性排好序的秩 1 层，适合处理数据（图像、点云、传感器记录）。两种读法说的是同一组 $\sigma, \mathbf u, \mathbf v$。

6再看数学

定义

任意一个 $m\times n$ 的实矩阵 $A$（不必是方阵，不必可逆）都可以写成

$$ A = \hl{U}{U}\,\hl{S}{\Sigma}\,\hl{Vt}{V^\top}, \qquad U \in \R^{m\times m},\;\; V\in\R^{n\times n},\qquad \Sigma = \operatorname{diag}(\sigma_1,\dots,\sigma_p) \in \R^{m\times n}, $$ $$ \sigma_1 \ge \sigma_2 \ge \dots \ge \sigma_p \ge 0, \qquad p = \min(m, n). $$

其中 $U$ 和 $V$ 是 1.10 节的正交矩阵：列向量两两垂直、长度为 1，作用是旋转（可能带镜像）。$\Sigma$ 只在对角线上有数。把 $AV = U\Sigma$ 一列一列地写出来，就回到了 FIG 1 的发现：

$$ A\,\hl{v1}{\cv{\mathbf v_1}} = \hl{s1}{\sigma_1}\,\hl{u1}{\cv{\mathbf u_1}}, \quad A\,\hl{v2}{\cw{\mathbf v_2}} = \hl{s2}{\sigma_2}\,\hl{u2}{\cw{\mathbf u_2}}, \quad \dots,\quad A\mathbf v_i = \sigma_i\mathbf u_i. $$

它们从哪里来：$A^\top A$ 的特征分解

代入 $A = U\Sigma V^\top$，利用 $U^\top U = I$：

$$ A^\top A = V\,\Sigma^\top\Sigma\,V^\top = V \operatorname{diag}(\sigma_1^2,\dots,\sigma_n^2)\, V^\top, \qquad A A^\top = U\,\Sigma\Sigma^\top U^\top. $$

这正是 1.11 节的特征分解：$V$ 的列是 $A^\top A$ 的特征向量，$\sigma_i^2$ 是它的特征值；$U$ 的列是 $AA^\top$ 的特征向量。$A^\top A$ 是对称矩阵，对称矩阵的特征向量天然正交，这就解释了为什么 $\cv{\mathbf v_1} \perp \cw{\mathbf v_2}$。它还是半正定的（$\mathbf x^\top A^\top A\mathbf x = |A\mathbf x|^2 \ge 0$），所以特征值非负，开方得到的 $\sigma_i$ 是实数。求出 $\mathbf v_i$ 以后，$\mathbf u_i = A\mathbf v_i/\sigma_i$。

也可以这样理解：$|A\mathbf x|^2 = \mathbf x^\top (A^\top A)\mathbf x$。FIG 1 里那条曲线的平方，就是对称矩阵 $A^\top A$ 的二次型，它的最大值、最小值落在特征向量上。

能直接读出来的量

$$ \|A\|_2 = \sigma_1, \qquad \operatorname{rank} A = \#\{\,i : \sigma_i > 0\,\}, \qquad \kappa(A) = \frac{\sigma_1}{\sigma_p}, \qquad |\det A| = \sigma_1\sigma_2\cdots\sigma_n \;\;(m = n). $$

最后一条很直观：旋转不改变面积，只有第 ② 步的拉伸改变面积，倍数是 $\sigma_1\sigma_2$。$\det A$ 的符号则藏在 $U$、$V$ 是否带镜像里。

最佳低秩近似（Eckart–Young）

$$ A_k = \sum_{i=1}^{k} \sigma_i\,\mathbf u_i\mathbf v_i^\top, \qquad \|A - A_k\|_2 = \sigma_{k+1}, \qquad \|A - A_k\|_F = \sqrt{\sigma_{k+1}^2 + \dots + \sigma_p^2}. $$

而且没有任何秩 $\le k$ 的矩阵能做得比这更好。直观地说：扔掉最薄的那几层，误差就是被扔掉的那几层的厚度。FIG 4 里的相对误差就是用第二个式子算的。

伪逆预告

求逆就是把三步倒着走：先用 $U^\top$ 转回来，再把每个方向除以 $\sigma_i$，最后用 $V$ 转回去。遇到 $\sigma_i = 0$ 的方向没法除，那就干脆放弃这个方向：

$$ A^{+} = V\,\Sigma^{+}\,U^\top, \qquad \Sigma^{+} = \operatorname{diag}\!\Big(\tfrac{1}{\sigma_1}, \dots, \tfrac{1}{\sigma_r}, 0, \dots, 0\Big). $$

这就是 Moore–Penrose 伪逆，1.13 节会看到它给出最小二乘解。注意 $1/\sigma_i$：奇异值越小，逆在那个方向上放大得越厉害。

三维：球变成椭球

一切照搬到三维：$3\times3$ 矩阵把单位球变成椭球，三个半轴长是 $\sigma_1\ge\sigma_2\ge\sigma_3$。秩 2 的矩阵把球压成一张椭圆盘，秩 1 压成一根针。选一个预设，或者编辑矩阵，再拖动阶段滑块；拖动画面可以转动视角。

FIG 5单位球 → 椭球

矩阵 A（3×3）

预设

黄、紫、青三根箭头从 $\mathbf v_1, \mathbf v_2, \mathbf v_3$ 出发，最后变成椭球的三根半轴 $\sigma_i\mathbf u_i$。细的红绿蓝箭头是 $\ihat,\jhat,\khat$ 的去向。这里的 3×3 SVD 用单边 Jacobi 法数值求出：反复对 $A$ 的两列做平面旋转，直到所有列两两正交，列长就是 $\sigma_i$。

实际怎么算

教科书上“先求 $A^\top A$ 的特征值”的做法适合理解，不适合编程：构造 $A^\top A$ 会把条件数平方，小奇异值的精度大量丢失。LAPACK（NumPy、Eigen 背后的库）先把 $A$ 化为双对角形，再迭代求解（Golub–Kahan 方法）；小矩阵常用 Jacobi 法。直接调用 np.linalg.svd 就好。

7机器人中的它：可操作度椭圆与点云法向量

平面两连杆机械臂

两连杆机械臂的末端位置是 $x = \ell_1\cos\theta_1 + \ell_2\cos(\theta_1+\theta_2)$，$y = \ell_1\sin\theta_1 + \ell_2\sin(\theta_1+\theta_2)$。对两个关节角求偏导，得到 $2\times2$ 的雅可比：

$$ \dot{\mathbf x} = J\dot{\mathbf q}, \qquad J = \begin{bmatrix} \ci{-\ell_1 s_1 - \ell_2 s_{12}} & \cj{-\ell_2 s_{12}} \\ \ci{\ell_1 c_1 + \ell_2 c_{12}} & \cj{\ell_2 c_{12}} \end{bmatrix}, \qquad |\det J| = \sigma_1\sigma_2 = \ell_1\ell_2\,|\sin\theta_2|. $$

（$s_1 = \sin\theta_1$，$c_{12} = \cos(\theta_1+\theta_2)$，依此类推。）第 1 列是只转关节 1 时末端的速度，第 2 列是只转关节 2 时的速度。所有满足 $|\dot{\mathbf q}| \le 1$ 的关节速度，被 $J$ 映射成末端的速度椭圆：半轴为 $\sigma_i\mathbf u_i$。拖动肘部或末端，或者用滑块调关节角。

FIG 6两连杆机械臂的可操作度椭圆

连杆长 $\ell_1 = 1.2$，$\ell_2 = 1.2/\sqrt2 \approx 0.85$；椭圆按 0.5 倍缩放画出。红、绿短箭头是 $J$ 的两列。手臂伸直（$\theta_2 = 0$）或完全折叠（$\theta_2 = \pm180°$）时，两列平行，椭圆塌成线段：末端无法沿手臂方向运动。淡色圆环是工作空间的内外边界，奇异位形恰好都落在边界上。

把机械臂带到接近奇异的位形（$\sigma_2/\sigma_1 < 0.05$）。椭圆的短轴指向哪里？
找到一个各向同性的姿态：椭圆几乎是正圆（$\kappa < 1.05$）。这时 $\theta_2$ 是多少？

力椭圆是速度椭圆的“对偶”。末端施加的力 $\mathbf F$ 需要的关节力矩是 $\boldsymbol\tau = J^\top\mathbf F$（3.4 节会推导）。在力矩预算 $|\boldsymbol\tau|\le1$ 下，末端能施加的力构成另一个椭圆：主轴方向同样是 $\mathbf u_i$，但半轴长变成 $1/\sigma_i$。速度上最“懒”的方向，恰好是力上最“壮”的方向，就像自行车的低速挡：转得慢，但蹬得动。切换上图的“力椭圆”看看。手臂伸直时你能单手撑起很重的东西，就是这个道理。

真实系统里

Yoshikawa 可操作度 $w = \sqrt{\det(JJ^\top)} = \sigma_1\sigma_2\cdots$ 常被用作运动规划的代价项，让机械臂远离奇异位形。实时逆运动学求解器会监控 $\sigma_{\min}$：它一变小，$J^{-1}$ 就会要求巨大的关节速度，所以要改用阻尼最小二乘，把 $1/\sigma$ 换成 $\sigma/(\sigma^2+\lambda^2)$（1.13 节和 3.4 节）。

点云 PCA：法向量和主轴

深度相机拍到桌面上的一小片点。想知道桌面的法向量，可以把这 $N$ 个点减去均值，按行堆成一个 $N\times3$ 的数据矩阵 $X$，然后做 SVD：$X = U\Sigma V^\top$。$V$ 的三列是点云的三个主方向，$\sigma_i/\sqrt N$ 是点云沿该方向的标准差。点云铺在一个平面上，所以沿法向几乎没有厚度：最小奇异值对应的 $\mathbf v_3$ 就是法向量。反过来，一根杆子形状的点云，最大奇异值对应的 $\mathbf v_1$ 就是它的主轴。

FIG 7用 SVD 估计法向量和主轴

灰点是 160 个采样点，三根箭头是 $\mathbf v_i$，长度为 2 倍标准差 $2\sigma_i/\sqrt N$（太短的按最小长度画）。黄色箭头是我们要的那个方向（平面模式下是 $\mathbf v_3$，杆模式下是 $\mathbf v_1$），绿色虚线是真实方向。拖动画面可以转动视角。

在平面模式下加大噪声，直到估计的法向量偏离真值超过 20°。这时 $\sigma_2$ 和 $\sigma_3$ 是什么关系？

PCL（Point Cloud Library）和 Open3D 的法向量估计做的正是这件事：对每个点取 k 个邻居，求协方差矩阵 $X^\top X/N$（其特征值就是 $\sigma_i^2/N$），取最小特征值的特征向量。

8常见误区

误区 1：奇异值就是特征值

不是。特征向量是被 $A$ 送回自己方向的向量（$A\mathbf v = \lambda\mathbf v$）；奇异向量是被送到另一个正交方向的向量（$A\mathbf v_i = \sigma_i\mathbf u_i$）。旋转 90° 没有实特征向量，但两个奇异值都是 1；剪切 $\begin{bmatrix}1&1\\0&1\end{bmatrix}$ 的特征值是 1 和 1，奇异值却是 1.618 和 0.618。两者只在对称半正定矩阵上一致（此时 $U = V$，正是 FIG 2 的第三个任务）。

误区 2：奇异向量是唯一的

把 $\mathbf v_i$ 和 $\mathbf u_i$ 同时取反，$A\mathbf v_i = \sigma_i\mathbf u_i$ 依然成立，所以不同的库、甚至同一个库在相邻两帧，可能给出符号相反的结果。如果两个奇异值相等（椭圆是圆），对应平面内任意一对正交方向都行。不要逐元素比较两次 SVD 的结果；点云法向量要用“朝向相机”之类的规则统一符号。

误区 3：小奇异值无害

求逆要除以 $\sigma_i$。$\sigma_2 = 10^{-6}$ 意味着数据里 $10^{-6}$ 量级的噪声会被放大成 1 量级的解。常用对策：把小于容差的 $\sigma$ 当作 0（截断 SVD），或者用阻尼 $\sigma/(\sigma^2+\lambda^2)$ 代替 $1/\sigma$。

误区 4：U 和 V 是同一组基

$V$ 生活在输入空间，$U$ 生活在输出空间。对 7 自由度机械臂的 $6\times7$ 雅可比，$V$ 是 $7\times7$（关节空间方向），$U$ 是 $6\times6$（末端运动方向），连维数都不一样。“机械臂在哪个方向动不了”要看 $U$ 的最后一列；“哪种关节运动不起作用”要看 $V$ 的最后几列（零空间）。

9小练习

1. $A = \begin{bmatrix}3&0\\0&-2\end{bmatrix}$ 的奇异值是？

奇异值是长度，永远非负。负号是一次镜像，归到 $U$ 里：$A = \begin{bmatrix}1&0\\0&-1\end{bmatrix}\begin{bmatrix}3&0\\0&2\end{bmatrix}I$。9 和 4 是 $A^\top A$ 的特征值，即 $\sigma_i^2$。

2. 某个姿态下，两连杆机械臂雅可比的奇异值是 $\sigma_1 = 1.2,\ \sigma_2 = 0.004$。下面哪句话对？

$J\mathbf v_2 = \sigma_2\mathbf u_2$，所以 $\dot{\mathbf q} = \mathbf v_2/\sigma_2$ 才能产生单位速度 $\mathbf u_2$，其大小为 $1/0.004 = 250$。机械臂已经非常接近奇异位形，$\kappa = 300$。

3. 对哪一类矩阵，SVD 和特征分解完全是同一回事？

对称矩阵有正交的特征向量 $A = Q\Lambda Q^\top$；如果特征值还都 $\ge 0$，这就已经是一个 SVD（$U = V = Q$，$\Sigma = \Lambda$）。旋转矩阵的奇异值全是 1，特征值却是复数。

4. 桌面上一片去均值的点云，数据矩阵的奇异值为 $(9.1,\ 6.3,\ 0.05)$。桌面的法向量是？

点沿法向几乎没有散布，所以法向是“拉伸最少”的输入方向 $\mathbf v_3$。$\mathbf u_i$ 是 $N$ 维向量（每个点一个分量），不是三维空间里的方向。

10延伸阅读

G. Strang, Introduction to Linear Algebra, Ch. 7 — 奇异值分解
G. Strang, MIT 18.06 Linear Algebra (OCW) — 第 29 讲：奇异值分解
G. Strang, Linear Algebra and Learning from Data, Part I — 从数据角度讲 SVD、PCA 与低秩近似
L. N. Trefethen & D. Bau, Numerical Linear Algebra, Lectures 4–5 — SVD 的几何意义与数值性质
K. Lynch & F. Park, Modern Robotics, §5.4 — 可操作度椭球与力椭球
T. Yoshikawa, “Manipulability of Robotic Mechanisms”, Int. J. Robotics Research 4(2), 1985 — 可操作度指标的原始论文