Singular Value Decomposition: SVD, 奇异值分解

矩阵 \(A\) 的奇异值分解指的是把 \(A\) 分解为三个矩阵的乘积,\(A = UDV^\trsps\),使得 \(U\) 和 \(V\) 的列向量都是正交的,而 \(D\) 是对角矩阵且对角元素为正实数。

实践中常常出现 \(A\) 很接近于一个低秩矩阵的情形 (对于方阵来说这意味着近乎不可逆),此时为 \(A\) 找一个低秩近似就很有用。事实上,对于任意 \(k\),\(A\) 的奇异值分解给出对 \(A\) 的最佳 \(k\)-秩近似。

常见本征值分解只能用于方阵,并且要求这方阵满足一些条件以保证本征矢量的正交性,而奇异值分解可用于任何矩阵。事实上,奇异值分解给出的矩阵 \(V\) 的列总是组成一个正交集合,而对 \(A\) 无任何要求。\(V\) 的列矢量称为 \(A\) 的右奇异矢量。同样,\(U\) 的列称为 \(A\) 的左奇异矢量并且组成正交集合。如果 \(A\) 是可逆方阵,其逆是 \(VD^{-1}U^\trsps\)。

设 \(A\) 是个 \(n\) 行 \(d\) 列矩阵 (\(n\times d\))。把它想象为 \(d\) 维空间的 \(n\) 个点,然后我们需要找到能“最好地”代表这些点的 \(k\) 维子空间。这里“最好”的意思是,要让所有点到这个子空间的垂直距离的平方和最小。我们从子空间是一维的情况开始,也就是一根穿过坐标原点的直线。之后我们将会看到,最佳拟合的 \(k\) 维子空间可以通过 \(k\) 次使用最佳直线拟合算法得到。这个问题被称为最佳最小二乘拟合。

考虑一个点 \(\V{a}_i = (a_{i1}, a_{i2}, \ldots, a_{id})\)。把它投影到一条穿过原点的直线时,下面的等式成立 \begin{equation} a_{i1}^2 + a_{i2}^2 + \cdots + a_{id}^2 = (\text{点在线上的投影长度})^2 + (\text{点到线的距离})^2. \end{equation} 考虑到左边的数都是给定的,于是,最小化所有点到线的距离平方和的问题,等价于最大化所有点在线上的投影长度平方和。在一般的最佳拟合子空间问题里也是这样。

关于“最佳拟合”,也不是非要通过最小化距离平方和来定义。比如也可以通过最小化距离的和来定义。但是前者的好处是平方比绝对值常常有更方便的数学性质,比如上面把最小化距离平方和转变成最大化投影长度平方和就是一个例子。

奇异矢量

现在定义 \(n\times d\) 矩阵 \(A\) 的奇异矢量。像前面说的那样,把 \(A\) 的行设想为 \(d\) 维空间的 \(n\) 个点,然后考虑用经过原点的直线来做最佳拟合。令 \(\V{v}\) 是沿着这条直线的单位矢量,那么 \(A\) 的第 \(i\) 行,\(\V{a}_i\) 在 \(\V{v}\) 方向上的投影是 \(\norm{\V{a}_i\cdot\V{v}}\)。于是可以看出,投影长度的平方和是 \(\norm{A\V{v}}^2\),而最佳拟合直线让它最大。

现在来定义 \(A\) 的第一奇异矢量 \(\V{v}_1\)。它是个单位列向量,使得 \(A\) 的所有行向量在它的方向上投影长度平方和最大,也就是 \begin{equation} \V{v}_1 = \argmax_{\norm{\V{v}}=1} \norm{A\V{v}}. \end{equation} 技术上讲,\(\argmax\) 给出的结果不一定唯一,这种情况下我们随便选出一个作为第一奇异矢量就好了。

定义 \(A\) 的第一奇异值为 \(\sigma_1(A) = \norm{A\V{v}_1}\)。注意到 \[\sigma_1^2 = \sum_{i=1}^n(\V{a}_i\cdot\V{v})^2,\] 也就是所有点在第一奇异矢量方向上投影的平方和。

现在考虑寻求 \(A\) 的最佳 2 维子空间拟合。一种贪婪算法是,采用 \(\V{v}_1\) 作为这个子空间的第一个基矢量,然后需要找到第二个基矢量。关于最佳拟合的定义中的距离平方和有助于下面的推演。对于每个包含 \(\V{v}_1\) 的 2 维子空间,矩阵 \(A\) 往这个子空间的投影距离的平方和等于往 \(\V{v}_1\) 方向上的投影长度的平方和加上往这个子空间内与 \(\V{v}_1\) 垂直的方向上的投影长度的平方和。因此,与其说找包含 \(\V{v}_1\) 的最佳子空间,我们来找垂直于 \(\V{v}_1\) 的单位矢量 \(\V{v}_2\),使得 \(\norm{A\V{v}}\) 最大。使用相同的贪婪策略可找出最佳三维、四维,以及更高维的最佳子空间,这也定义出 \(\V{v}_3\),\(\V{v}_4\), \(\ldots\)。下面会给出正式的定义。并没有先验的理由能保证贪婪算法给出最佳拟合,但事实上可以证明按照前述操作得到的一系列子空间的确就是最佳拟合。

第二奇异矢量定义为垂直于 \(\V{v}_1\) 的最佳拟合矢量 \begin{equation} \V{v}_2 = \argmax_{\substack{\V{v}\bot\V{v}_1\\ \norm{\V{v}}=1}} \norm{A\V{v}}. \end{equation} \(\sigma_2(A) = \norm{A\V{v}_2}\) 被称为 \(A\) 的第二奇异值。第三奇异矢量和奇异值定义为 \begin{equation} \V{v}_3 = \argmax_{\substack{\V{v}\bot\V{v}_1,\ \V{v}_2\\ \norm{\V{v}}=1}} \norm{A\V{v}}, \end{equation} \begin{equation} \sigma_3 = \norm{A\V{v}_3}. \end{equation} 如此继续下去,直到当找到 \(r\) 个奇异矢量和奇异值之后,出现 \begin{equation} \max_{\substack{\V{v}\bot\V{v}_1,\ \V{v}_2,\ \ldots, \V{v}_r\\ \norm{\V{v}}=1}} \norm{A\V{v}} = 0 \end{equation} 的情况。

这也意味着,\(\V{v}_1,\ \V{v}_2,\ \ldots,\ \V{v}_r\) 张成 \(A\) 的所有行向量组成的空间。否则,假设某一行 \(\V{a}_j\) 不能表示成这些奇异矢量的线性组合。考虑 \(\V{b} = \V{a}_j - \left[(\V{a}_j\cdot\V{v}_1)\V{v}_1 + (\V{a}_j\cdot\V{v}_2)\V{v}_2 + \cdots + (\V{a}_j\cdot\V{v}_r)\V{v}_r\right]\),\(\V{b}\) 不是零矢量,且 \(\V{b}\) 垂直于 \(\V{v}_1,\ \V{v}_2,\ \ldots,\ \V{v}_r\),与前提矛盾。于是我们也知道 \(A\) 的秩是 \(r\)。

上面的做法是,先找出最佳一维子空间,然后找出包含此一维子空间的最佳二维子空间。如果我们一开始就直接找最佳二维子空间,得到的会不会更好呢?答案是否定的。下面给出贪婪算法对每个维数都给出最佳子空间的证明。

\(\theorem{4.1}\) \(A\) 是个 \(n\times d\) 的矩阵,奇异矢量为 \(\V{v}_1,\ \V{v}_2,\ \ldots,\ \V{v}_r\)。对于 \(k\le r\),令 \(V_k\) 表示 \(\V{v}_1,\ \V{v}_2,\ \ldots,\ \V{v}_k\) 张成的子空间。那么,对每个 \(k\),\(V_k\) 是 \(A\) 的最佳拟合 \(k\) 维子空间。

方法是用数学归纳法,要点是,当按前述方法找出最佳 \(k-1\) 维子空间 \(V_{k-1}\) 后,对于最佳 \(k\) 维子空间 \(W\), 需要找出 \(W\) 内的一个矢量 \(\V{w}_k\) 垂直于 \(V_{k-1}\),并且把 \(\V{w}_k\) 作为 \(W\) 的一个基。把 \(V_{k-1}\) 的 \(k-1\) 个正交基投影到 \(W\) 中,然后在 \(W\) 中找出一个与这 \(k-1\) 个投影都垂直的矢量就可以了。这意味着求解 \(k-1\) 个齐次线性方程。矢量的维数 \(d\) 必定不小于 \(k\),也就是 \(d > k-1\),所以方程一定有非零解。然后,按照之前描述的构造法则,\(\norm{A\V{w}_k}\le \norm{A\V{v}_k}\)。于是,\(A\) 在 \(W\) 的其它基矢量上的投影和加上在 \(\V{w}_k\) 上的投影和必定不超过 \(\norm{A\V{v}_1}^2 + \norm{A\V{v}_2}^2 + \cdots + \norm{A\V{v}_k}^2\) (这里用到了归纳假设),这就证明了 \(V_k\) 是最优的。

\(A\V{v}_i\) 作为一个 \(n\) 维矢量,其分量是 \(A\) 的行在 \(\V{v}_i\) 方向上的投影长度 (有正负号)。把 \(\norm{A\V{v}_i}\) 想象为矩阵 \(A\) 在 \(\V{v}_i\) 方向上的“分量”。这么想是有道理的,因为 \(A\) 在所有这些奇异矢量上的分量的平方和等于 \(A\) 的所有元素的平方和。对每个行向量 \(\V{a}_j\),\(\sum_{i=1}^{r}(\V{a}_j\cdot\V{v}_i)^2 = \norm{\V{a}_j}^2\)。然后对所有行向量求和就得到结论。

所以我们可以说,所有奇异值的平方和等于 \(A\) 的“全部内容”。根据这个事实我们可以定义出 \(A\) 的 Frobenius 范数 \begin{equation} \lVert A\rVert_\text{F} = \sqrt{\sum_{j,k} a_{jk}^2}. \end{equation} 按定义,Frobenius 范数的平方等于所有奇异值的平方和。

\(\V{v}_1,\ \V{v}_2,\ \ldots,\ \V{v}_r\) 定义的是右奇异矢量。 左奇异矢量定义为 \begin{equation} \V{u}_i = \frac{1}{\sigma_i(A)} A \V{v}_i, \label{eqLeftSingularVectorDef} \end{equation} 这里 \(i\) 取值为 1 到 \(r\)。需要证明左奇异矢量的两个性质:1. 它们也组成正交基;2. \(\V{u}_1,\ \V{u}_2,\ \ldots,\ \V{u}_r\) 像 \(\V{v}_1,\ \V{v}_2,\ \ldots,\ \V{v}_r\) 那样依次最大化 \(\norm{\V{u}^\trsps A}\)。

先来证明 \(\V{u}_1,\ \V{u}_2,\ \ldots,\ \V{u}_r\) 组成正交基。方法是数学归纳法,归纳变量是 \(A\) 的秩。当 \(A\) 的秩为 1,只有一个左或者右奇异矢量,当然成立。现在进行归纳。考虑 \begin{equation} B = A - \sigma_1 \V{u}_1\V{v}_1^\trsps. \end{equation} 然后会发现,对 \(B\) 的奇异值分解步骤,相当于 \(A\) 的从第二个奇异值开始的奇异值分解步骤。为了看出这点,首先观察到,按定义,\(B\V{v}_1 = 0\)。于是 \(B\) 的第一个右奇异矢量——记为 \(\V{z}\)——必定垂直于 \(\V{v}_1\)。否则,如果有个 \(\V{v}_1\) 方向的分量——记为 \(\V{z}_1\),会有 \(\norm{B \frac{\V{z}-\V{z}_1}{\norm{\V{z}-\V{z}_1}}} = \frac{\norm{B\V{z}}}{\norm{\V{z}-\V{z}_1}} > \norm{B\V{z}}\)。注意到,“分量”的定义是 \(\V{z}_1 = (\V{z}\cdot\V{v}_1)\V{v}_1\),于是 \(\norm{\V{z}-\V{z}_1}^2 = 1 - 2(\V{z}\cdot\V{v}_1)^2 + (\V{z}\cdot\V{v}_1)^2 = 1 - (\V{z}\cdot\V{v}_1)^2 \le 1\)。 这就与 \(\V{z}\) 定义里的 \(\argmax\) 矛盾。所以 \(\V{z}\bot\V{v}_1\)。可是这样一来,\(B\V{z} = A\V{z}\),于是 \(B\) 的第一个奇异矢量是 \(A\) 的第二个奇异矢量;可如此顺次推下去。

这样,跑一遍算法得到 \(B\) 的右奇异矢量 \(\V{v}_2,\ \V{v}_3,\ \ldots,\ \V{v}_r\) 和相应的左奇异矢量 \(\V{u}_2,\ \V{u}_3,\ \ldots,\ \V{u}_r\)。按照归纳假设,\(\V{u}_2,\ \V{u}_3,\ \ldots,\ \V{u}_r\) 是正交的。

剩下的就是证明 \(\V{u}_1\) 与 \(\V{u}_2,\ \V{u}_3,\ \ldots,\ \V{u}_r\) 正交。否则的话,对某个 \(i\ge 2\),\(\V{u}_1^\trsps \V{u}_i\ne 0\)。不妨假定 \(\V{u}_1^\trsps \V{u}_i > 0\)。对于无穷小的 \(\varepsilon > 0\),下面这个矢量 \begin{equation} A\left(\frac{\V{v}_1 + \varepsilon\V{v}_i}{\norm{\V{v}_1 + \varepsilon\V{v}_i}}\right) = \frac{\sigma_1\V{u}_1 + \varepsilon\sigma_i\V{u}_i}{\sqrt{1+\varepsilon^2}} \end{equation} 的长度不会短于其沿着 \(\V{u}_1\) 方向分量的长度 \begin{equation} \begin{split} &\V{u}_1^\trsps\cdot\frac{\sigma_1\V{u}_1 + \varepsilon\sigma_i\V{u}_i}{\sqrt{1+\varepsilon^2}} = \left(\sigma_1 + \varepsilon\sigma_i\V{u}_1^\trsps\V{u}_i\right) \left(1-\frac{\varepsilon^2}{2}+O(\varepsilon^4\right) \\ =& \sigma_1 + \varepsilon\sigma_i\V{u}_1^\trsps\V{u}_i - O(\varepsilon^2) > \sigma_1, \end{split} \end{equation} 这是个矛盾,因为 \(\sigma_1\) 本来应该是单位矢量里 \(\norm{A\V{v}}\) 的最大值。这就证明了 \(\V{u}_1,\ \V{u}_2,\ \ldots,\ \V{u}_r\) 是正交的。

奇异值分解 (SVD)

\(\sigma_i\V{u}_i\) 是一个矢量,其分量是 \(A\) 的每行在 \(\V{v}_i\) 上的投影。\(\sigma_i\V{u}_i\V{v}_i^\trsps\) 是一个秩为 1 的矩阵,每一列相当于 \(\sigma_i\V{u}_i\) 乘以相应的 \(\V{v}_i\) 的坐标。

我们要证明,\(A\) 可以按下面的方式分解为一些秩为 1 的矩阵的和 \begin{equation} A = \sum_{i=1}^r \sigma_i\V{u}_i\V{v}_i^\trsps. \end{equation} 先证明一个引理:\(A\V{v} = B\V{v}\) 对所有矢量 \(\V{v}\) 都成立,当且仅当 \(A = B\)。只要把 \(\V{v}\) 的元素每次取成只有一个为 1 其余都为 0 就容易看出。

现在,任何矢量 \(\V{v}\) 可以写成奇异矢量的线性组合加上一个与所有奇异矢量垂直的矢量,\(\V{v} = \sum_{i=1}^r c_i\V{v}_i + \V{v}_\bot\)。而按照奇异矢量的定义,\(A\V{v}_\bot = \V{0}\)。于是 \(A\V{v} = \sum_{i=1}^rc_i A\V{v}_i = \sum_{i=1}^r c_i\sigma_i\V{u}_i\),这与 \(\left(\sum_{i=1}^r \sigma_i\V{u}_i\V{v}_i^\trsps\right) \V{v}\) 相等,于是得证。

\(A\) 的这种分解叫做奇异值分解。用矩阵记号可写为 \(A = U D V^\trsps\),这里 \(U\) 和 \(V\) 的列对应于左右奇异矢量,而 \(D\) 是个对角矩阵,对角元素是 \(A\) 的奇异值。\(U\) 是 \(n\) 行 \(r\) 列,\(D\) 是 \(r\) 行 \(r\) 列,而 \(V\) 是 \(r\) 行 \(d\) 列。\(A_{st} = \sum_{i=1}^r\sigma_i \V{u}_{i,s}\V{v}_{i,t} = \sum_{i=1}^r\sum_{j=1}^r \V{u}_{j,s}D_{j,i} \V{v}_{i,t}\),所以 \(U_{sj} = \V{u}_{j,s},\ V_{it}=\V{v}_{i,t}\),或者可以直观地想象为 \(U\) 里面 \(\V{u}_i\) 是竖着放的,而 \(V\) 里面 \(\V{v}_i\) 是横着放的。

矩阵的奇异值序列是唯一的。如果没有任何两个奇异值相等,那么奇异矢量序列也是唯一的。但是,如果一些奇异值相等,则相应的奇异矢量组成一个子空间。张成这个子空间的任何正交基都可作为奇异矢量。

根据奇异值分解,我们发现 \(\V{u}_i^\trsps A = \sigma_i\V{v}_i^\trsps\),所以 \(\norm{\V{u}_i^\trsps A} = \sigma_i\)。\(\V{u}_1\) 必定是所有单位矢量里让 \(\norm{\V{u}^\trsps A}\) 最大的。这可归结为求 \[ \argmax_{\substack{\{c_i,\ i=1,\ \ldots,\ r\}\\\sum_{i=1}^r c_i^2=1}}\sum_{i=1}^r c_i^2\sigma_i^2 \] 的问题,其答案是显然的。类似地,\(\V{u}_2\) 必定是所有与 \(\V{u}_1\) 垂直的单位矢量里让 \(\norm{\V{u}^\trsps A}\) 最大的。可如此推下去,于是我们证明了 \(\V{u}_i\) 序列与 \(\V{v}_i\) 序列有类似的性质。

最佳 \(k\) 秩近似

有两种重要的矩阵范数,一种是前面定义的 Frobenius 范数 \(\lVert A\rVert_\text{F}\),一种是 2-范数 \(\lVert A\rVert_2\)。后者定义为 \begin{equation} \max_{\norm{\V{v}}=1} \norm{A\V{v}}, \end{equation} 也就是等于最大的奇异值。回顾前面讲的,Frobenius 范数是矩阵代表的所有点到原点距离的平方和。而 2-范数是往某个方向的投影距离的平方和的最大值,也就是沿着第一奇异矢量的投影距离的平方和的最大值。

对于 \(k\in\{1,\ 2,\ \ldots,\ r\}\),令 \begin{equation} A_k = \sum_{i=1}^k \sigma_i\V{u}_i\V{v}_i^\trsps, \end{equation} 也就是 \(A\) 的奇异值分解的前 \(k\) 项。\(A_k\) 的秩为 \(k\)。\(A_1\) 相当于一个列向量乘以一个常数后横向重复了 \(d\) 次,因此一个列向量乘以一个行向量得到的矩阵的秩总是 1 (除非其中一个是零矢量,这样得到的矩阵为零矩阵)。一个矩阵的秩定义为这个矩阵的线性独立列向量的个数 (等于线性独立行向量的个数),等于这个矩阵的列数减去以这个矩阵为系数的齐次微分方程组的解空间的维数。比如当一个方阵可逆,则只有零解,解空间维数为 0,秩为其行数或列数。\(A_k\) 的列数为 \(d\)。以 \(A_k\) 为系数的齐次线性方程组的解空间由 \(\V{v}_{k+1},\ \V{v}_{k+2},\ \ldots,\ \V{v}_r\),以及 \(d-r\) 个垂直于 \(\V{v}_1,\ \V{v}_2,\ \ldots,\ \V{v}_r\) 的基矢量张成,所以解空间的维数是 \(r-k+d-r=d-k\),因此 \(A_k\) 的秩为 \(d-(d-k)=k\)。 下面我们证明,\(A_k\) 是 \(A\) 的最佳 \(k\) 秩近似,这里近似误差通过 Frobenius 范或者 2-范数定义。

\(A_k\) 的行可看成是由 \(A\) 的行投影到前 \(k\) 个奇异矢量张成的子空间得到。这个从 \(A_k\) 的定义即可得到。

\(\theorem{4.7}\) 对于任何秩不超过 \(k\) 的矩阵 \(B\),下式成立 \begin{equation} \froben{A-A_k} \le \froben{A-B}. \end{equation} 证明:令 \(B\) 是秩不超过 \(k\) 的让 \(\froben{A-B}^2\) 最小的矩阵。令 \(V\) 表示 \(B\) 的行张成的矩阵。\(V\) 的维数顶多为 \(k\)。\(B\) 的每行必定等于相应的 \(A\) 的行在 \(V\) 上的投影;否则,把 \(B\) 的一行换成相应的 \(A\) 的一行在 \(V\) 上的投影可以减小 \(\froben{A-B}^2\),给定一个向量和一个子空间,求子空间里的一个向量使得它与给定向量的差的长度最短。答案就是,当这个向量等于给定矢量在子空间里的投影时最短;通过把待求矢量表示为子空间的正交基的线性叠加可以证明这个。同时不改变 \(V\) (也就不会改变 \(B\) 的秩)。这样一来,\(\froben{A-B}^2\) 就是 \(A\) 的行向量与 \(V\) 的距离的平方和。之前已经证明,\(A\) 的最佳拟合 \(k\) 维子空间是 \(V_k\),也就是说,当 \(V=V_k\),\(A\) 的行向量与 \(V\) 的距离的平方和最小。换句话说,此时,\(B\) 的行向量张成的子空间必定等于 \(V_k\);再考虑到 \(A_k\) 的确等于 \(A\) 在 \(V_k\) 上的投影,于是必须有 \(B=A_k\),得证。

接下来考虑 2-范数。

引理: \begin{equation} \twonorm{A-A_k}^2 = \sigma_{k+1}^2. \end{equation} 证明:\(A-A_k = \sum_{i=k+1}^r \sigma_i\V{u}_i\V{v}_i^\trsps\)。我们需要求一个单位矢量 \(\V{v} = \sum_{i=1}^r \alpha_i\V{v}_i\) 使得 \(\norm{(A-A_k) \V{v}}\) 最大。 \begin{equation} \begin{split} &\norm{(A-A_k)\V{v}} = \norm{\sum_{i=k+1}^r \alpha_i\sigma_i\V{u}_i\V{v}_i^\trsps\V{v}_i}\\ =& \norm{\sum_{i=k+1}^r\alpha_i\sigma_i\V{u}_i} = \sqrt{\sum_{i=k+1}^r\alpha_i^2\sigma_i^2}, \end{split} \end{equation} 其中 \(\sum_{i=1}^r\alpha_i^2 = 1\)。显然,当 \(\alpha_{k+1}=1\) 而别的都为零时结果最大。

\(\theorem{4.9}\) 对于任何秩不超过 \(k\) 的矩阵 \(B\),下式成立 \begin{equation} \twonorm{A-A_k} \le \twonorm{A-B}. \end{equation} 证明:如果 \(A\) 的秩不超过 \(k\),则上式左边为 0,结论显然。所以假定 \(A\) 的秩大于 \(k\)。假定存在秩不超过 \(k\) 的矩阵 \(B\) 使得 \(\sigma_{k+1} > \twonorm{A-B}\)。\(B\) 的零空间,可以把一个矩阵的零空间想象为垂直于这个矩阵的所有行向量的向量张成的子空间。也就是以 \(B\) 为系数的齐次线性方程组的解张成的空间,其维数至少是 \(d-k\)。令 \(\V{v}_1,\ \V{v}_2,\ \ldots,\ \V{v}_{k+1}\) 表示 \(A\) 的前 \(k+1\) 个奇异矢量。基于维数的考虑,存在非零矢量 \(\V{z}\) 属于下述集合 \[ \text{Null}(B) \cap \text{Span}\left\{\V{v}_1,\ \V{v}_2,\ \ldots,\ \V{v}_{k+1}\right\}, \] 也就是这两个集合的交集不为空。令 \(\V{z}_1,\ \V{z}_2,\ \ldots,\ \V{z}_{d-k}\) 为 \(\text{Null}(B)\) 里的 \(d-k\) 个线性独立矢量。现在,\(\V{z}_1,\ \V{z}_2,\ \ldots,\ \V{z}_{d-k},\ \V{v}_1,\ \V{v}_2,\ \ldots,\ \V{v}_{k+1}\) 是 \(d\) 维空间里的 \(d+1\) 个矢量,因此线性相关。于是可以找到 \(\alpha_1,\ \alpha_2,\ \ldots,\ \alpha_{d-k}\) 以及 \(\beta_1,\ \beta_2,\ \ldots,\ \beta_{k+1}\),使得 \(\sum_{i=1}^{d-k}\alpha_i\V{z}_i = \sum_{j=1}^{k+1}\beta_j\V{v}_j\)。原书此处把上标弄错了。令 \(\V{z}=\sum_{i=1}^{d-k}\alpha_i\V{z}_i\) 且把 \(\V{z}\) 归一化。我们要证明,\(\norm{(A-B)\V{z}}\ge\sigma_{k+1}\)。首先,按照 2-范数的定义 \begin{equation} \twonorm{A-B}^2 \ge \norm{(A-B)\V{z}}^2. \end{equation} 由于 \(\V{z}\) 在 \(B\) 的零空间内,\(B\V{z}=\V{0}\)。于是 \begin{equation} \twonorm{A-B}^2 \ge \norm{A\V{z}}^2. \end{equation} 同时,由于 \(\V{z} = \sum_{j=1}^{k+1}\beta_j\V{v}_j\),把 \(A\) 的奇异值分解代入即可看出上式不会小于 \(\sigma_{k+1}^2\)。得证。

\(\theorem{4.10}\) (本征值和本征矢量的类比) \[A\V{v}_i = \sigma_i \V{u}_i\ \text{and}\ A^\trsps \V{u}_i = \sigma_i \V{v}_i.\] 第一个等式从左奇异矢量的定义得到 (方程(\ref{eqLeftSingularVectorDef}))。对于第二个等式,注意到从SVD得到\[A^\trsps \V{u}_i = \sum_j\sigma_j\V{v}_j\V{u}_j^\trsps\V{u}_i,\] 结合 \(\V{u}_j\)的正交归一性即得。

计算奇异值分解的幂方法

由 \(A\) 的奇异值分解可得 \begin{equation} B \equiv A^\trsps A = \sum_i \sigma_i^2 \V{v}_i\V{v}_i^\trsps. \end{equation} \(B\) 是一个对称方阵,于是左右奇异矢量相同。继续求平方 \begin{equation} B^2 = \sum_i \sigma_i^4 \V{v}_i\V{v}_i^\trsps. \end{equation} 一般情形 \begin{equation} B^k = \sum_i \sigma_i^{2k} \V{v}_i\V{v}_i^\trsps. \end{equation} 如果 \(\sigma_1 > \sigma_2\),则 \begin{equation} \frac{1}{\sigma_1^{2k}} B^k \to \V{v}_1\V{v}_1^\trsps. \end{equation} 注意到 \(\froben{B^k}\to \sigma_1^{2k}\),所以,可以先计算 \(B^k\),然后除以 \(\froben{B^k}\),再把得到的矩阵的第一个列向量归一化,就得到对 \(\V{v}_1\) 的近似。

问题在于,当 \(A\) 是个巨大 (但稀疏) 的矩阵时,比如大小为 \(10^8\times10^8\) 但仅包含 \(10^9\) 个非零元素,得到的 \(B\) 会是一个 \(10^8\times10^8\) 且稠密的矩阵。这种情况下 \(B\) 的存储都成问题,何况计算它的幂。即使 \(A\) 没那么大,计算矩阵的幂仍然是个慢速过程。因此需要更有效率的方法。

给定一个随机矢量 \(\V{x} = \sum c_i\V{v}_i\),我们有 \begin{equation} B^k\V{x} \to \sigma_1^{2k}\V{v}_1\V{v}_1^\trsps\;\sum c_i\V{v}_i = \sigma_1^{2k}c_1\V{v}_1. \end{equation} 归一化上面这个矢量就得到 \(\V{v}_1\)。注意到 \(B^k\V{x}\) 的计算只是 \(2k\) 次稀疏矩阵乘以矢量的计算,计算量小多了。

为了计算\(k\)个奇异矢量,选取随机矢量\(\V{r}\),然后为矢量空间\(\V{r},\ A\V{r},\ \ldots,\ A^{k-1}\V{r}\)选取正交归一基。然后计算\(A\)与每个基矢量的乘积,然后找出新的矢量空间的正交基。直观地说,\(A\)被应用到一个子空间上,而不是单个矢量上。不断地把\(A\)应用到子空间,然后正交归一化,很快就会收敛到前\(k\)个奇异矢量。

当第一个和第二个奇异值差别不大时,上面的方法收敛会比较慢。但下面会证明在这种情况幂方法会收敛到“接近最大”的那些奇异值对应的奇异矢量。

引理:令 \((x_1,\ x_2,\ \ldots,\ x_d)\) 是一个 \(d\) 维随机单位矢量。那么 \(\norm{x_1} \ge \frac{1}{20\sqrt{d}}\) 的几率至少是 \(9/10\)。

\(\theorem{4.11}\) 令\(A\)是\(n\times d\)矩阵,\(\V{x}\)是\(\mathbb{R}^d\)空间中的单位长度矢量,且\(\norm{\V{x}^\trsps \V{v}_1}\ge\delta\),这里\(\delta > 0\)。令\(V\)是\(A\)的对应的奇异值大于\((1-\varepsilon)\sigma_1\)的左奇异矢量张成的空间。令\(\V{w}\)是幂方法应用\(k=\frac{\ln(1/\varepsilon\delta)}{2\varepsilon}\)次后的单位矢量,即 \[\V{w} = \frac{\left(A^\trsps A\right)^k\V{x}}{\norm{\left(A^\trsps A\right)^k\V{x}}},\] 则\(\V{w}\)的垂直于\(V\)的分量的大小最多为\(\varepsilon\)。这就是说,如果一开始跟\(\V{v}_1\)不是完全垂直,则迭代之后跟对应于较大奇异值的那些奇异矢量张成的空间也不会很垂直。

\(\theorem{4.12}\) 令\(\V{y}\in\mathbb{R}^n\)是服从球对称单位方差高斯分布的随机矢量。将\(\V{y}\)归一化,即令\(\V{x} = \V{y}/\norm{\V{y}}\)。对于任意的单位矢量\(\V{v}\),有 \[\text{Prob}\left(\norm{\V{x}^\trsps\V{v}}\le \frac{1}{20\sqrt{d}}\right) \le \frac{1}{10} + 3 e^{-d/96}.\] 也就是说,随机选取的两个矢量几乎互相垂直的概率并不高。

奇异值分解的应用

主成分分析

奇异值分解(SVD)的传统应用是主成分分析(PCA)。考虑推荐电影的问题。有\(n\)个顾客,\(d\)个电影,矩阵\(A\)的元素\(a_{ij}\)表示顾客\(i\)对电影\(j\)的喜爱程度。假定有\(k\)个决定顾客对电影的喜欢程度的基本因素,并且\(k\)远小于\(d\)和\(n\)。比如,这些因素可以包括:喜剧、戏剧性、动作、新奇性的多少。每个电影可以通过一个\(k\)维矢量来描述其包含的每个因子的多少,而每个顾客可以通过一个\(k\)维矢量来描述其心目中不同因子的重要性,这两个矢量的点乘就表示一个顾客对一个电影的喜爱程度。也就是说,\(n\times d\)的矩阵\(A\)可以表示为描述顾客的\(n\times k\)矩阵\(U\)与描述电影的\(k\times d\)矩阵\(V\)的乘积。通过SVD找出最佳的\(k\)秩近似\(A_k\)给出\(U\)和\(V\)。\(A\)未必会精确等于\(U V\),这种情况下\(A - U V\)被当成噪声。另一个问题是,SVD会给出包含负数的矩阵。非负矩阵分解(NMF)在某些情况下是更恰当的选择。

在上面的描述中,\(A\)是可完整获得的,我们希望通过\(U\)和\(V\)来找出基础矢量。但在现实情形,比如在电影推荐的情形,每个顾客可能只看了少量电影,也就是说更自然的假设是我们只知道\(A\)的一部分元素,然后希望估计完整的\(A\)。如果\(A\)是个一般的\(n\times d\)矩阵,则包含的信息量是\(\Omega(n d)\)的量级,因此不可能用很少的数来估计。如果假定\(A\)是一个低秩矩阵叠加上了一部分噪音,并且其元素是从某个分布随机取样的,则有可能可以用SVD来估计\(A\)的总体。这被称为协同过滤,这被用于推荐系统。

球形高斯分布混合的聚类

谱分解

奇异矢量和文件分级

在一个离散优化问题里的应用

奇异矢量与本征矢量