奇异值分解（SVD）在机器学习、信号处理、统计学以及金融等领域中具有广泛应用。

但你有没有想过它是怎么来的呢？我猜十有八九没看到过吧，没关系，本文很快就告诉你。

在 1870 年代由意大利数学家贝尔特拉米（Beltrami）和法国数学家若尔当（Jordan）引入。英国数学家西尔维斯特（Sylvester）于 1889 年独立于贝尔特拉米和若尔当提出了对矩阵进行了奇异值分解。

在刚提出的那个时候，都是针对实数方阵而言的。1902 年，由 Autonne 引入了复数矩阵，并在 1939 年由 Eckhart 和 Young 引入了一般矩阵（即实数/复数和方阵/非方阵）。

而本文主要看一看在早期提出 SVD 的时候，有哪些数学家分别做了什么工作。目的是加深对 SVD 的了解，以便能更好地应用这个线代中的强大工具。

1引言

矩阵理论中最重要的思想之一就是矩阵分解。矩阵分解的理论实用性早已受到数学家们的肯定。而随着计算机的发明，它们逐渐成为数值线性代数的支柱，它们已经成为可以解决各种问题的计算平台。

而矩阵分解的发展离不开所谓的规范型。为了理清这一发展思路，让我们从简要介绍相关的历史背景开始。

大多数经典的矩阵分解方法的提出时间是要早于矩阵概念得到广泛使用的年代。矩阵分解法是随着行列式、线性方程组，尤其是双线性形式和二次型等问题的研究而提出来的。

高斯可以说是做了一些早期工作。他在 1823 年撰写的论文中使用他早在 1801 年就使用的消元法来实现如下任务，

具体而言，可以将函数（关于的一个二次型）简化为以下形式：

\frac{u^{0} u^{0}}{\mathcal{A}^{0}}+\frac{u^{\prime} u^{\prime}}{\mathcal{B}^{\prime}}+\frac{u^{\prime \prime} u^{\prime \prime}}{\mathcal{C}^{\prime \prime}}+\frac{u^{\prime \prime \prime} u^{\prime \prime \prime}}{\mathcal{D}^{\prime \prime \prime}}+\mathrm{etc.}+M

其中，除数等是常数，而等是等的线性函数。但是，第二个函数独立于；第三个独立于和；第四个独立于和，依此类推。最后一个函数仅取决于最后一个未知数。此外，系数在中分别乘以。

从中我们可以看到高斯这里的方法将二次型的矩阵分解为乘积，其中是对角矩阵，而是与具有相同对角元素的上三角矩阵。

高斯这里的等函数是向量对应的元素。

另外，高斯通过消元法还能有效地获得矩阵的逆，从而将方程组转换为逆线性系统。高斯处理二次型和线性方程组的技巧使他对最小二乘法的理论和应用的一般化处理成为可能。

随后，针对不同的问题出现了一系列进展，

1829 年，柯西（Cauchy）通过考虑二次型和相应的齐次方程组，建立了对称矩阵的特征值和特征向量的特性。
1846 年，雅可比（Jacobi）给出了著名的对称矩阵对角化算法，在 1857 年的遗作中，他通过以高斯形式分解双线性形式获得了 LU 分解。
1868 年，魏尔斯特拉斯（Weierstrass）为一对双线性函数建立了规范型，该问题现在被称为广义特征值问题。因此，从一定角度看来，1873 年提出的奇异值分解被视为规范型上取得的众多成果之一。

2预备知识

双线性形式

向量空间上的双线性形式是域上的一个双线性映射。换句话说，双线性形式是一个函数，关于它的每个参数都是线性的，即有

B(\mathbf{u}+\mathbf{v}, \mathbf{w})=B(\mathbf{u}, \mathbf{w})+B(\mathbf{v}, \mathbf{w}) \quad \text {和} \quad B(\lambda \mathbf{u}, \mathbf{v})=\lambda B(\mathbf{u}, \mathbf{v})

B(\mathbf{u}, \mathbf{v}+\mathbf{w})=B(\mathbf{u}, \mathbf{v})+B(\mathbf{u}, \mathbf{w}) \quad \text {和} \quad B(\mathbf{u}, \lambda \mathbf{v})=\lambda B(\mathbf{u}, \mathbf{v})

坐标表示

令为一个 -维向量空间，它的一组基为。

由定义的矩阵称为双线性形式在基上的矩阵。

矩阵表示向量在这组基上的坐标，并且类似地，表示另一个向量，则，

B(\mathbf{v}, \mathbf{w})=\mathbf{x}^{\top} \mathbf{A} \mathbf{y}=\sum_{i, j=1}^{n} x_{i} a_{i j} y_{j}.

同一个双线性形式在不同的基上对应不同的矩阵。但是，一个双线性形式在不同的基上的矩阵相互之间是合同的。换句话说，假设是的另一组基，则有

\mathbf{f}_{j}=\sum_{i=1}^{n} C_{i, j} \mathbf{e}_{i}

其中形成一个可逆矩阵。然后，该双线性形式在新基上的矩阵为。

3三位作者的独立工作

1、问题

我们将使用现代矩阵符号来描述有关奇异值分解的早期工作。与高斯对其分解的描述一样，大多数内容都很容易转化为使用矩阵术语。

你要知道的是，使用现代矩阵记法可能会让你觉得他们的成就貌似微不足道。但是完全以原始的标量形式来表示的话，对于现代的读者来说会不习惯。这里，为了便以阅读，我们将几位作者的标记符号等统一起来。

在本文中，我们将针对如下所示的奇异值分解问题，

\mathbf{A}=\mathbf{U} \mathbf{\Sigma} \mathbf{V}^{\mathrm{T}}

其中是阶实矩阵，

\mathbf{\Sigma}=\operatorname{diag}\left(\sigma_{1}, \sigma_{2}, \ldots, \sigma_{n}\right)

是对角矩阵，其对角线元素非负，并且按从大到下排列。而矩阵和

\mathbf{U}=\left(\mathbf{u}_{1} \mathbf{u}_{2} \cdots \mathbf{u}_{n}\right) \quad \text { 以及 } \quad \mathbf{V}=\left(\mathbf{v}_{1} \mathbf{v}_{2} \cdots \mathbf{v}_{n}\right)

都是正交矩阵。后面用到的函数表示如下定义的 Frobenius-范数，

\|\mathbf{A}\|^{2}=\sum_{i, j} a_{i j}^{2}=\sum_{i} \sigma_{i}^{2}.

要知道，三位作者的论文在内容上比这里提到的精华要丰富得多，如果想看作者们详尽的论述，建议读者参考原始论文资料。

2、Beltrami，1873 年

Beltrami 和 Jordan 被认为是奇异值分解的共同开创者，Beltrami 第一个发表相关论文，而 Jordan 稍晚一点发表，但后者的贡献在于问题处理的完整性和优雅性。Beltrami 的贡献发表在意大利大学学生使用的数学杂志上，其目的是鼓励学生熟悉双线性形式。

推导

Beltrami 从如下双线性形式开始，

f(\mathbf{x}, \mathbf{y})=\mathbf{x}^{\mathrm{T}} \mathbf{A} \mathbf{y}

其中，为阶实矩阵。如果引入如下变量代换，

\mathbf{x}=\mathbf{U} \boldsymbol{\xi} \quad \text {以及} \quad \mathbf{y}=\mathbf{V} \boldsymbol{\eta} \qquad (2.0)

则有，

f(\mathbf{x}, \mathbf{y})=\xi^{\mathrm{T}} \mathbf{S} \eta

其中，

\mathbf{S}=\mathbf{U}^{\mathrm{T}} \mathbf{A} \mathbf{V} \qquad (2.1)

Beltrami 观察到，如果要求和是正交的，这样的话在选择它们的元素时将有个自由度，他建议用这些自由度来消除中的非对角元素。

假设是对角矩阵，即

\mathbf{S}=\mathbf{\Sigma}=\operatorname{diag}\left(\sigma_{1}, \ldots, \sigma_{n}\right).

然后从式和的正交性可得出，

\mathbf{U}^{\top} \mathbf{A}=\mathbf{\Sigma} \mathbf{V}^{\top} \qquad (2.2)

同样的有，

\mathbf{A V}=\mathbf{U} \boldsymbol{\Sigma} \qquad (2.3)

从式可得到（其中是将矩阵对角线上的元素取倒数），将其代入式，可得，

\left(\mathbf{A} \mathbf{A}^{\top}\right)\mathbf{U}=\mathbf{U}\mathbf{\Sigma}^{2} \qquad (2.4)

类似地，我们可以得到，

\left(\mathbf{A}^{\mathrm{T}} \mathbf{A}\right) \mathbf{V}=\mathbf{V} \mathbf{\Sigma}^{2}

因此，是如下方程的根，

\operatorname{det}\left(\mathbf{A} \mathbf{A}^{\top}-\sigma^{2} \mathbf{I}\right)=0 \qquad (2.5)

和

\operatorname{det}\left(\mathbf{A}^{\top} \mathbf{A}-\sigma^{2} \mathbf{I}\right)=0 \qquad (2.6)

请注意，Beltrami 这里的推导假设是非奇异的，因此也是非奇异的。

Beltrami 现在认为两个函数和是相同的，因为它们都是次多项式，并且在的值，以及在的值相同。

Beltrami 接下来声明，根据一个众所周知的定理，式的根是实数，并且他们是正的。为了说明这一点，他指出

0<\left\|\mathbf{x}^{\top} \mathbf{A}\right\|^{2}=\mathbf{x}^{\top} \mathbf{AA}^{\top} \mathbf{x}=\boldsymbol{\xi}^{\top} \boldsymbol{\Sigma}^{2} \boldsymbol{\xi} \qquad (2.7)

这里需要用到两个等式，由式可得，以及由式可得。另外，上面的不等式也意味着这里的是正的。

Beltrami 给出如下对角化的算法步骤，

1、计算式的根。
2、从式确定出。Beltrami 在此处指出，的列是可以决定的，只差一个因子，但实际上这是在所有不相同时才成立。他还默认生成的是正交的，这也要求是不同的。
3、根据式得出，这步要求非奇异。

讨论

从前面的内容可以明显看出，Beltrami 对于一个非奇异的、具有不同奇异值的实方阵推导出了奇异值分解。他的推导有一个小缺陷就是缺少了对于退化情况的处理。Beltrami 省略了这些额外情况，可能是在为他的学生简化过程，但也有可能他确实没有考虑到这些问题。

3、Jordan，1874 年

Jordan 可以称为奇异值分解的共同提出者。尽管他在 Beltrami 出版一年后才发表他的著作，但很显然，他的这项工作是由不同问题导出的独立工作。

推导

Jordan 从如下双线性形式开始，

P=\mathbf{x}^{\top} \mathbf{A} \mathbf{y}

并计算在如下约束下的最大值和最小值，

\|\mathbf{x}\|^{2}=\|\mathbf{y}\|^{2}=1 \qquad (3.1)

最大值由如下公式确定

0=d P=d \mathbf{x}^{\top} \mathbf{A} \mathbf{y}+\mathbf{x}^{\top} \mathbf{A} d \mathbf{y} \qquad (3.2)

而该式对于满足下式的所有和都成立，

d \mathbf{x}^{\top} \mathbf{x}=\mathbf{0} \quad \text {和} \quad d \mathbf{y}^{\top} \mathbf{y}=0 \qquad (3.3)

Jordan 然后断言: 方程将是方程的组合，从中可以得出

\mathbf{A} \mathbf{y}= \sigma \mathbf{x} \qquad (3.4)

以及

\mathbf{x}^{\mathrm{T}} \mathbf{A}=\tau \mathbf{y}^{\mathrm{T}} \qquad (3.5)

从式可得所求的最大值为，

\mathbf{x}^{\mathrm{T}}(\mathbf{A} \mathbf{y})=\sigma \mathbf{x}^{\mathrm{T}} \mathbf{x}=\sigma

同样地，最大值也为，因此有。

上面是 Jordan 的推导，但感觉到有些地方或许并不是很清晰。

为了大家更好地理解，本文改用拉格朗日乘子法。先构造如下拉格朗日函数，

\Phi = \mathbf{x}^{\top} \mathbf{A} \mathbf{y} - \frac{\sigma}{2} \left(\mathbf{x}^{\top} \mathbf{x} - 1\right) - \frac{\tau}{2} \left(\mathbf{y}^{\top} \mathbf{y}-1\right)

为了下文方便，上式两个乘子分别除以 2。

\begin{array}{c} \displaystyle \frac{\partial \Phi}{\partial \mathbf{x}} = \mathbf{Ay} - \sigma \mathbf{x} = \mathbf{0} \\[1em] \Downarrow \\[1em] \mathbf{Ay} = \sigma \mathbf{x} \end{array}

上面最后的式子两边左乘，因为是单位向量，因此得，

\mathbf{x}^{\top} \mathbf{A} \mathbf{y} = \sigma \mathbf{x}^{\top}\mathbf{x} = \sigma

可见，目标双线性形式的最大值为。

\begin{array}{c} \displaystyle \frac{\partial \Phi}{\partial \mathbf{y}} = \mathbf{x}^{\top}\mathbf{A} - \tau \mathbf{y}^{\top} = \mathbf{0}\\[1em] \Downarrow \\[1em] \mathbf{x}^{\top}\mathbf{A} = \tau \mathbf{y}^{\top} \end{array}

上面最后的式子两边右乘得，

\mathbf{x}^{\top} \mathbf{A} \mathbf{y} = \tau \mathbf{y}^{\top}\mathbf{y} = \tau

同样地，最大值也为，因此有。因此式可以改为，

\mathbf{A} \mathbf{y}= \sigma \mathbf{x} \qquad (3.4)

\mathbf{A}^{\top}\mathbf{x}=\sigma \mathbf{y} \qquad (3.5\star)

接下来看由式和式组成的线性方程组，

\left[\begin{array}{cc} -\sigma \mathbf{I} & \mathbf{A} \\ \mathbf{A}^{\mathbf{T}} & -\sigma \mathbf{I} \end{array}\right] \left[\begin{array}{cc} \mathbf{x} \\ \mathbf{y} \end{array}\right] = \mathbf{0}

若尔当观察到，是上面线性方程组对应的行列式决定的，即

D=\left|\begin{array}{cc} -\sigma \mathbf{I} & \mathbf{A} \\ \mathbf{A}^{\mathbf{T}} & -\sigma \mathbf{I} \end{array}\right|

要有非零解，该行列式必须为零。他显示该行列式仅包含的偶次幂。现在让为方程式的根，让和满足式和式，其中。（Jordan 指出，即使它不是唯一的，也可以找到这样的解。）

令

\hat{\mathbf{U}}=\left(\mathbf{u} \; \mathbf{U}_{\star}\right) \quad \text {和} \quad \hat{\mathbf{V}}=\left(\mathbf{v} \; \mathbf{V}_{\star}\right)

是正交的，并令

\mathbf{x}=\hat{\mathbf{U}} \hat{\mathbf{x}} \quad \text {和} \quad \mathbf{y}=\hat{\mathbf{V}} \hat{y}

通过这样的代换，得

P=\hat{\mathbf{x}}^{\top} \hat{\mathbf{A}} \hat{\mathbf{y}}

当时，其中

e_{1}=(1,0, \ldots, 0)^{\top},

取得最大值。而且，在取最大值处由、和的值容易验证得到如下结论，

\hat{\mathbf{A}} \hat{\mathbf{y}}=\sigma_{1} \hat{\mathbf{x}} \quad \text {以及} \quad \hat{\mathbf{x}}^{\top} \hat{\mathbf{A}}=\sigma_{1} \hat{\mathbf{y}}^{\top}

这意味着

\hat{\mathbf{A}}=\left[\begin{array}{cc} \sigma & 0 \\ 0 & {\mathbf{A}}_{1} \end{array}\right].

于是，令以及，因此，的形式为

\sigma_{1} \xi_{1} \eta_{1}+P_{1}

其中独立于和。接着，Jordan 用归纳法得出目标双线性形式的规范型为，

P=\xi^{\mathrm{T}} \mathbf{\Sigma} \boldsymbol{\eta}

最后，Jordan 指出，当特征方程的根方便计算时，可以直接从式、式和式计算出和的列。

讨论

Jordan 以简洁、优雅的手法来处理该问题。在数学上可以称为 deflation。他先解决一部分问题，从而成功将问题的规模缩小，以至于用归纳法解决整个问题，也避免了 Beltrami 方法中的的退化问题。

顺便说一句，舒尔（Schur）后来使用同样的策略建立了对一般矩阵的三角化。

另外，形成行列式的矩阵

\left[\begin{array}{cc} 0 & \mathbf{A} \\ \mathbf{A}^{\mathrm{T}} & \mathbf{0} \end{array}\right]

也得到了广泛使用。

4、Sylvester，1889 年

Sylvester 就奇异值分解的主题写了一个 Note 和两篇论文。Sylvester 在一篇论文中提出了一种将二次型简化为标准型的迭代算法，他指出，类似的迭代方法可以用于对角化双线性形式。在 1889 年的一篇论文中，Sylvester 提出了迭代算法和具体法则。

法则

Sylvester 从如下双线性形式开始

{B}=\mathbf{x}^{\top} \mathbf{A} \mathbf{y}.

并考虑二次型

M=\sum_{i}\left(\frac{d {B}}{d y_{i}}\right)^{2} = \mathbf{x}^{\top}\mathbf{A} \mathbf{A}^{\mathrm{T}} \mathbf{x}.

令为的规范型。如果具有规范型，则等价于，这意味着在一定排序下有。

为了找到代换，Sylvester 引入矩阵和，并断言对的代换是对角化的代换，对的代换是对角化的代换。一般来说，仅当具有不同的奇异值时，这点才成立。

无穷小迭代

Sylvester 提出此方法是为了对角化二次型，后来他将其推广到双线性形式。对于二次型，它已经足够复杂，这里将限于这种情况。

Sylvester 使用归纳法，即先假定已经解决阶为时的问题。因此，对于，可以假设矩阵的形式为，

\mathbf{A}=\left[\begin{array}{lll} \color{blue}{a} & \color{blue}{0} & \color{red}{f} \\ \color{blue}{0} & \color{blue}{b} & \color{red}{g} \\ \color{red}{f} & \color{red}{g} & \color{red}{c} \end{array}\right]

左上角蓝色的的矩阵是归纳法假设已经搞定了。接下来的问题是在保留先前引入的的情况下消除和。

Sylvester 建议进行以下形式的所谓无穷小正交变换

\left[\begin{array}{c} x_{1} \\ x_{2} \\ x_{3} \end{array}\right]=\left[\begin{array}{ccc} 1 & \epsilon & \eta \\ -\epsilon & 1 & \theta \\ -\eta & -\theta & 1 \end{array}\right]\left[\begin{array}{l} \xi_{1} \\ \xi_{2} \\ \xi_{3} \end{array}\right]

其中矩阵非对角线的元素非常小，因此比 1 次高的幂次项可以忽略。然后，转换矩阵

\left[\begin{array}{ccc} 1 & -\epsilon & -\eta \\ \epsilon & 1 & -\theta \\ \eta & \theta & 1 \end{array}\right] \left[\begin{array}{lll} \color{blue}{a} & \color{blue}{0} & \color{red}{f} \\ \color{blue}{0} & \color{blue}{b} & \color{red}{g} \\ \color{red}{f} & \color{red}{g} & \color{red}{c} \end{array}\right] \left[\begin{array}{ccc} 1 & \epsilon & \eta \\ -\epsilon & 1 & \theta \\ -\eta & -\theta & 1 \end{array}\right]

中第和两个元素为

(a-b) \epsilon-f \theta-g \eta \qquad (4.1)

而的变化由

\frac{1}{2} \delta\left(f^{2}+g^{2}\right)=(a-c) f \eta+(b-c) g \theta

如果或均不为零，则可以适当选择和的值来减小。如果非零，则可以选择，以使式等于 0，即保留先前引入的 0。对于诸如之类的特殊情况，Sylvester 也提出了如何通过显式缩小问题规模来处理的办法。

Sylvester 现在声称，这些无穷小转换的无穷序列将使或之一减为零，或者将问题减少为一种特殊情况。

讨论

Sylvester 在当时并不了解上面 Jordan 在 1874 年的工作以及雅可比（Jacobi）在 1846 年提出的针对二次型对角化问题的迭代算法（该算法在 1954 年被 Kogbetliantz 推广到了奇异值分解）。Sylvester 在没有证明的情况下将其描述出来，从而给读者留下了太多细节。

4小结

由前文我们知道，数学家从二次型的研究导出了矩阵特征值问题以及特征值分解。而从本文可以知道，从双线性形式的规范型问题导出了矩阵的奇异值分解这一强大的理论和应用工具。

本文主要就早期提出奇异值分解的几位数学家的工作做了简单介绍，目的是对它有更多的了解。但我们学习这一理论的主要目的还是在于应用，这也是这个主题接下来的重点，敬请关注这个系列。

相关阅读

矩阵前传 - 克莱姆没能证明的法则被他两行搞定

矩阵前传 - 矩阵之父 Sylvester 为什么提出 Matrix

矩阵前传 - 柯西-比内公式及其用初等矩阵的证明

二次型和矩阵合同原来是这么一回事

拉格朗日乘子法的来历与直观解释

矩阵特征值是这么来的，以及有趣的盖尔圆

万能的 SVD 分解是哪位牛人提出来的？

1引言

2预备知识

双线性形式

坐标表示

3三位作者的独立工作

1、问题

2、Beltrami，1873 年

推导

讨论

3、Jordan，1874 年

推导

讨论

4、Sylvester，1889 年

法则

无穷小迭代

讨论

4小结

矩阵和线性代数原来是这么来的

概率论原来可以这样优雅地入门

机器学习的数学基础之向量范数

机器学习的数学基础之矩阵范数

拉格朗日乘子法的来历与直观解释

万能的 SVD 分解是哪位牛人提出来的？

1引言

2预备知识

双线性形式

坐标表示

3三位作者的独立工作

1、 问题

2、Beltrami，1873 年

推导

讨论

3、Jordan，1874 年

推导

讨论

4、Sylvester，1889 年

法则

无穷小迭代

讨论

4小结

矩阵和线性代数原来是这么来的

概率论原来可以这样优雅地入门

机器学习的数学基础 之 向量范数

机器学习的数学基础 之 矩阵范数

拉格朗日乘子法的来历与直观解释

1、问题

机器学习的数学基础之向量范数

机器学习的数学基础之矩阵范数