矩阵分解及其Eigen实现_eigen lu分解

主要是用来记录自己的学习过程，内容也主要来自于网上的各种资料，然后自己总结而来，参考的资料都以注明，感谢这些作者的分享。如果内容有误，请大家指点。

LU分解¹

理论

定义

       将矩阵等价为两个矩阵$L$和$U$的乘积 ，其中$L$和$U$分别是单位下三角矩阵和上三角矩阵，当$A$的前rank(A)阶顺序主子式都不为0时，矩阵$A$可以分解为
$$A=LU$$
其中$L$是下三角矩阵，$U$是上三角矩阵。

LU分解求线性方程组

       对于求解线性方程，将$A$进行$LU$分解，得到$L(Ux)=b$，令$y=Ux$，则有$Ly=b$，由于$L$是下三角矩阵，可以方便求出$y$，求出$y$之后，在求解$Ux=y$，由于$U$是上三角矩阵，这样就可以方便求出$x$。

LU分解数值稳定值

       考虑一个矩阵 A = [ 0 1 1 1 ] A=\left[

$$\begin{array}{ll}0 & 1 \\ 1 & 1\end{array}$$ \right] A=[01​11​] ，虽然A非奇异，且条件数很小 $\kappa (A)\approx 2.62$ ，但是LU分解算法第一次执行就失败了，因为 $A_{0,0}=0$ 在分母上。
考虑另一个矩阵 $A=\left[\begin{array}{cc}10^{-20}& 1\\ 1& 1\end{array}\right]$ ，给 $A_{0,0}$ 一个误差。可以计算得到 $L=\left[\begin{array}{cc}1& 0\\ 10^{20}& 1\end{array}\right]$ ， $U=\left[\begin{array}{cc}10^{-20}& 1\\ 0& 1-10^{20}\end{array}\right]$ ，如果这里用的是浮点数的话，那么 $U\approx \left[\begin{array}{cc}10^{-20}& 1\\ 0& -10^{20}\end{array}\right]$ ， 把LU相乘得到 $LU=\left[\begin{array}{cc}10^{-20}& 1\\ 1& 0\end{array}\right]\ne \left[\begin{array}{cc}10^{-20}& 1\\ 1& 1\end{array}\right]=A$ 。如果用这组LU去计算问题，会出现条件数很小，但相对误差很大的情况。
注: 因为浮点数的固有问题, 在处理数值计算问题时候，一定要尽可能避免大数和小数相加减的问题。

LUP分解¹

理论

       在高斯消元的时候，做一些行置换（pivoting）。

LUP分解求线性方程组

• step1：LUP分解
• step2：求解$Ly=Pb$
• step3：求解$Ux=y$

LUP分解数值稳定性

LUP算法比LU算法稳定

LU分解的Eigen实现

       Eigen提供了两种LU分解的方法，分别为Eigen::FullPivLU和Eigen::PartialPivLU。其中，Eigen::PartialPivLU类提供了可逆方阵的LU分解，具有部分置换(partial pivoting)，将矩阵$A$分解为$A=PLU$，$L$为单位下三角矩阵，$U$为上三角矩阵，$P$为置换矩阵(permutation matrix)。但是，Eigen::FullPivLU类可以对任意矩阵进行LU分解，并具有完全置换(complete pivoting)，矩阵$A$分解为$A=P^{-1}LU{Q}^{-1}$，其中，$L$为单位下三角矩阵，$U$为上三角矩阵，$P,Q$为置换矩阵。

需要注意的是：Eigen::PartialPivLU类会断言矩阵是方阵，但是不会判断矩阵是否可逆，需要在使用时我们自己去判断。Eigen::FullPivLU类非常稳定，并且在大型矩阵上经过了很好的测试。但是，在某些用例中，SVD 分解本质上更稳定和/或更灵活。例如，在计算矩阵的核时，使用 SVD 可以选择矩阵的最小奇异值，这是 LU 分解看不到的。

       接下里，具体看看这两种算法在Eigen中是如何实现。

方法一：Eigen::FullPivLU

int main(int argc, char *argv[])
{
     Matrix3d A = Matrix3d::Random();
     cout << "初始矩阵A = " << '\n'
          << A << endl;
     Vector3d b = Vector3d::Random();
     cout << "b = " << '\n'
          << b << endl;
     Eigen::FullPivLU<Matrix3d> lu(A);
     Matrix3d U = lu.matrixLU().triangularView<Upper>();
     Matrix3d L = lu.matrixLU().triangularView<UnitLower>();
     Matrix3d P = lu.permutationP();
     Matrix3d Q = lu.permutationQ();
     cout << "L = \n"
          << L << endl;
     cout << "U = \n"
          << U << endl;
     cout << "P = \n"
          << P << endl;
     cout << "Q = \n"
          << Q << endl;
     cout << "重构原始矩阵后的结果为：\n"
          << P.inverse() * L * U * Q.inverse() << endl;
     cout << "x = \n"
          << lu.solve(b) << endl;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

结果分析：

初始矩阵A = 
 0.680375   0.59688 -0.329554
-0.211234  0.823295  0.536459
 0.566198 -0.604897 -0.444451
b = 
   0.10794
-0.0452059
  0.257742
L = 
        1         0         0
  0.72499         1         0
-0.734727  0.493089         1
U = 
 0.823295 -0.211234  0.536459
        0  0.833518 -0.718482
        0         0  0.303977
P = 
0 1 0
1 0 0
0 0 1
Q = 
0 1 0
1 0 0
0 0 1
重构原始矩阵后的结果为：
 0.680375   0.59688 -0.329554
-0.211234  0.823295  0.536459
 0.566198 -0.604897 -0.444451
x = 
  0.60876
-0.231282
  0.51038
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32

方法二：Eigen::PartialPivLU

int main(int argc, char *argv[])
{
     Matrix3d A = Matrix3d::Random();
     cout << "A = \n"
          << A << endl;
     cout << A.determinant() << endl;
     Eigen::PartialPivLU<Eigen::Matrix3d> plu(A);
     Matrix3d P = plu.permutationP();
     cout << "P = \n"
          << P << endl;
     Matrix3d L = plu.matrixLU().triangularView<UnitLower>();
     Matrix3d U = plu.matrixLU().triangularView<Upper>();
     cout << "L = \n"
          << L << endl;
     cout << "U = \n"
          << U << endl;
     cout << "重构后的系统矩阵A: \n"
          << P * L * U << endl;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

结果分析：

A = 
 0.680375   0.59688 -0.329554
-0.211234  0.823295  0.536459
 0.566198 -0.604897 -0.444451
0.208598
P = 
1 0 0
0 0 1
0 1 0
L = 
        1         0         0
 0.832185         1         0
-0.310467 -0.915573         1
U = 
 0.680375   0.59688 -0.329554
        0  -1.10161   -0.1702
        0         0  0.278313
重构后的系统矩阵A: 
 0.680375   0.59688 -0.329554
-0.211234  0.823295  0.536459
 0.566198 -0.604897 -0.444451
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

Cholesky分解¹

理论

定义

是LU分解的一种特殊情况，当系统矩阵$A$为对称正定矩阵时，此时可以得到$U=L^T$，即
$$A=L{L}^T$$

Choleskyfen分解数值稳定性

不需要想普通的矩阵一样还需要置换操作，因此相比于LU，LUP分解数值更加稳定。

Cholesky分解Eigen实现

       同样地，EIgen也提供了两种方法，分别为Eigen::LDLT以及Eigen::LLT。Eigen::LLT提供了对称正定矩阵的分解方法，使得$A=L{L}^T$，其中$L$为下三角矩阵。Eigen::LDLT是EIgen提供的更加鲁棒的Cholesky分解方法，只需要矩阵为半正定或者半负定即可，将$A$分解为$A=P^{T}LD{L}^{\ast }P$，其中，$P$为置换矩阵，$L$ 是具有单位对角线的下三角矩阵，$D$ 是对角矩阵，Eigen::LDLT这个类支持就地分解机制。推荐使用Eigen::LDLT。

方法一：Eigen::LLT

int main(int argc, char *argv[])
{
     MatrixXd A(3, 3);
     A << 4, -1, 2, -1, 6, 0, 2, 0, 5;
     cout << "The matrix A is" << endl
          << A << endl;
     LLT<MatrixXd> llt(A);
     MatrixXd L = llt.matrixL();
     cout << "L = \n"
          << L << endl;
     cout << "重构后的矩阵A: \n"
          << L * L.transpose() << endl;
}
1
2
3
4
5
6
7
8
9
10
11
12
13

结果分析：

The matrix A is
 4 -1  2
-1  6  0
 2  0  5
L = 
       2        0        0
    -0.5  2.39792        0
       1 0.208514   1.9891
重构后的矩阵A: 
 4 -1  2
-1  6  0
 2  0  5
1
2
3
4
5
6
7
8
9
10
11
12

方法二：Eigen::LDLT

int main(int argc, char *argv[])
{
     MatrixXd A(3, 3);
     A << 4, -1, 2, -1, 6, 0, 2, 0, 5;
     Vector3d b = Vector3d::Random();
     cout << "The matrix A is" << endl
          << A << endl;
     cout << "b = \n"
          << b << endl;
     LDLT<MatrixXd> ldlt(A);
     MatrixXd L = ldlt.matrixL();
     auto D = ldlt.vectorD();
     cout << "L = \n"
          << L << endl;
     cout << "D = \n"
          << D << endl;
     cout << "x = \n"
          << ldlt.solve(b);
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

结果分析：

The matrix A is
 4 -1  2
-1  6  0
 2  0  5
b = 
 0.680375
-0.211234
 0.566198
L = 
        1         0         0
        0         1         0
-0.166667       0.4         1
D = 
      6
      5
3.03333
x = 
   0.13803
-0.0122007
 0.0580278
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

QR分解²,³

理论

定义

       一个矩阵 $A\in {\mathbb{R}}^{m\times n},m\ge n$ 可以被分解成 $A=QR$, 其中:

1. $Q\in {\mathbb{R}}^{m\times m}$ 是正交矩阵
   $R\equiv \left[\begin{array}{c}\hat{R}\\ 0\end{array}\right]\in {\mathbb{R}}^{m\times n}$
2. $\hat{R}\in {\mathbb{R}}^{n\times n}$ 是上三角矩阵

正交矩阵的性质

3. $Q^{T}Q=Q{Q}^T=I$
4. 左乘一个正交矩阵对欧式范数的结果不影响
   $$\Vert Qv{\Vert }_2^2=v^T{Q}^{T}Qv=v^{T}v=\Vert v{\Vert }_2^2$$

QR分解求解线性最小二乘

       对于超定线性最小二乘问题$Ax\approx b$，目标函数为：
ϕ ( x ) = ∥ r ( x ) ∥ 2 2 = ∥ b − A x ∥ 2 2 = ∥ b − Q [ R ^ 0 ] x ∥ 2 2 = ∥ Q T ( b − Q [ R ^ 0 ] x ) ∥ 2 2 = ∥ Q T b − [ R ^ 0 ] x ∥ 2 2

$$\begin{aligned} \phi(x)=\|r(x)\|_2^2 & =\|b-A x\|_2^2=\left\|b-Q\left[\begin{array}{c} \hat{R} \\ 0 \end{array}\right] x\right\|_2^2 \\ & =\left\|Q^T\left(b-Q\left[\begin{array}{c} \hat{R} \\ 0 \end{array}\right] x\right)\right\|_2^2 \\ & =\left\|Q^T b-\left[\begin{array}{c} \hat{R} \\ 0 \end{array}\right] x\right\|_2^2 \end{aligned}$$ ϕ(x)=∥r(x)∥22​​=∥b−Ax∥22​= ​b−Q[R^0​]x ​22​= ​QT(b−Q[R^0​]x) ​22​= ​QTb−[R^0​]x ​22​​
其中，$Q\in {\mathbb{R}}^{m\times m}$，$Qb\in {\mathbb{R}}^{m\times 1}$， $\left[\begin{array}{c}\hat{R}\\ 0\end{array}\right]\in {\mathbb{R}}^{m\times n}$，$\left[\begin{array}{c}\hat{R}\\ 0\end{array}\right]x\in {\mathbb{R}}^{m\times 1}$。如果把$Q^{T}b$拆分成上下两部分，即$Q^{T}b=\left[\begin{array}{c}{c}_1\\ c_2\end{array}\right]$，其中$c_1\in {\mathbb{R}}^n$，$c_2\in {\mathbb{R}}^{m-n}$。那么目标函数可转换为：
$${\Vert r(x)\Vert }_2^2={\Vert c_1-\hat{R}x\Vert }_2^2+{\Vert c_2\Vert }_2^2$$
       显然， ${\Vert c_2\Vert }_2^2$与$x$无关，能优化的只有前半部分，使得${\Vert c_1-\hat{R}x\Vert }_2^2$等于0，即$\hat{R}x=c_1$，\left| r(x) \right|_2^2的最小值为${\Vert c_2\Vert }_2^2$。这样处理之后就可以避免正规方程中的$(A^{T}A{)}^{-1}$，避免了条件数变成$cond(A^{T}A)=cand(A{)}^2$，所有QR分解比正规方程求解最小二乘问题更加稳定。

计算QR分解的方法

1 Gram–Schmidt Orthogonalization
2 Householder Triangularization
3 Givens Rotations
当A是稠密矩阵，Givens Rotations并没有比另外两种算法更高效，但如果A是稀疏矩阵，那么Givens Rotations大小为0的元素可以直接被忽略。另一个优势是，Givens Rotations更容易并行化，因为Givens Rotations只对两个元素进行操作，处理不同列的时候可以完全的独立。

QR分解的Eigen实现

       Eigen提供了四种计算QR分解的方法，分别为：Eigen::ColPivHouseholderQR，Eigen::CompleteOrthogonalDecomposition，Eigen::FullPivHouseholderQR以及Eigen::HouseholderQR。

方法一：HouseholderQR

typedef Matrix<float, 5, 1> Vector5f;
int main(int argc, char *argv[])
{
     MatrixXf A(MatrixXf::Random(5, 3)), Q;
     A.setRandom();
     Vector5f b = MatrixXf::Random(5, 1);
     HouseholderQR<MatrixXf> qr(A);
     Q = qr.householderQ();
     std::cout << "The complete unitary matrix Q is:\n"
               << Q << "\n\n";
     cout << "x = \n"
          << qr.solve(b) << endl;
}
1
2
3
4
5
6
7
8
9
10
11
12
13

结果分析：

The complete unitary matrix Q is:
  -0.676275   0.0792999    0.713112  -0.0788488   -0.147031
  -0.220518   -0.322306   -0.370288   -0.365561   -0.759436
  -0.353086   -0.344724   -0.214153      0.8414  -0.0517703
    0.58236   -0.461852    0.555351    0.176282   -0.328724
  -0.173814   -0.746785 -0.00906669   -0.348021    0.539351

x = 
 -0.262592
-0.0749101
 -0.593202
1
2
3
4
5
6
7
8
9
10
11

方法二：Eigen::ColPivHouseholderQR

typedef Matrix<float, 5, 1> Vector5f;
int main(int argc, char *argv[])
{
     MatrixXf A(MatrixXf::Random(5, 3)), Q;
     A.setRandom();
     Vector5f b = MatrixXf::Random(5, 1);
     ColPivHouseholderQR<MatrixXf> cqr(A);
     Q = cqr.householderQ();
     MatrixXf R = cqr.matrixR().triangularView<Upper>();
     cout << "Q = \n"
          << Q << endl;
     cout << "R = \n"
          << R << endl;
     cout << "x = \n"
          << cqr.solve(b) << endl;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

结果分析：

Q = 
 -0.603651   0.704296   0.334273   0.165978  -0.016934
 -0.320874  -0.364108  -0.232567   0.833958  -0.122029
  -0.45273  -0.208744  -0.202056  -0.427952  -0.726286
  0.379609   0.572496  -0.623708   0.173782  -0.330052
  -0.42846 0.00820044  -0.635875  -0.252233   0.590251
R = 
   1.60258     1.3316   -1.15008
         0   0.860025 -0.0119055
         0          0   0.438341
         0          0          0
         0          0          0
x = 
 -0.262592
-0.0749101
 -0.593202
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

方法三：Eigen::FullPivHouseholderQR

typedef Matrix<float, 5, 1> Vector5f;
int main(int argc, char *argv[])
{
     MatrixXf A(MatrixXf::Random(5, 3)), Q;
     A.setRandom();
     Vector5f b = MatrixXf::Random(5, 1);
     FullPivHouseholderQR<MatrixXf> fqr(A);
     Q = fqr.matrixQ();
     cout << "Q = \n"
          << Q << endl;
     cout << "x = \n"
          << fqr.solve(b) << endl;
}
1
2
3
4
5
6
7
8
9
10
11
12
13

结果分析：

Q = 
  0.124977  -0.919134   0.334273   0.162086 -0.0395422
  0.467087   0.131775  -0.232567   0.826756  -0.163864
  0.493559 -0.0702728  -0.202056  -0.160453    0.82758
 -0.629484  -0.274961  -0.623708   0.274145   0.252941
   0.35547  -0.239345  -0.635875  -0.435088  -0.471929
x = 
 -0.262592
-0.0749099
 -0.593202
1
2
3
4
5
6
7
8
9
10

方法四：CompleteOrthogonalDecomposition

typedef Matrix<float, 5, 1> Vector5f;
int main(int argc, char *argv[])
{
     MatrixXf A(MatrixXf::Random(5, 3)), Q;
     A.setRandom();
     Vector5f b = MatrixXf::Random(5, 1);
     CompleteOrthogonalDecomposition<MatrixXf> cod(A);
     Q = cod.householderQ();
     MatrixXf T = cod.matrixT().triangularView<Upper>();
     MatrixXf Z = cod.matrixZ();
     cout << "Q = \n"
          << Q << endl;
     cout << "T = \n"
          << T << endl;
     cout << "Z = \n"
          << Z << endl;
     cout << "x = \n"
          << cod.solve(b) << endl;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

结果分析：

 = 
 -0.603651   0.704296   0.334273   0.165978  -0.016934
 -0.320874  -0.364108  -0.232567   0.833958  -0.122029
  -0.45273  -0.208744  -0.202056  -0.427952  -0.726286
  0.379609   0.572496  -0.623708   0.173782  -0.330052
  -0.42846 0.00820044  -0.635875  -0.252233   0.590251
T = 
   1.60258     1.3316   -1.15008
         0   0.860025 -0.0119055
         0          0   0.438341
         0          0          0
         0          0          0
Z = 
1 0 0
0 1 0
0 0 1
x = 
 -0.262592
-0.0749101
 -0.593202
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

SVD分解⁴,⁵

定义

       对于任意一个矩阵$A\in {\mathbb{C}}^{m\times n}$，可以分解为：
$$A=U\Sigma V^H$$
其中，$U$是$m\times m$酉矩阵，$V$是$N\times N$酉矩阵，$\Sigma$是$m\times n$对角矩阵，对角线上的元素称为矩阵$A$的奇异值($A^{T}A$或者$A^{T}A$特征值开根号)，矩阵$U$的列向量称为左奇异向量(left singular vector)，矩阵$V$的列向量称为右奇异向量(right singular vector)。$A$的左奇异向量是$A{A}^T$ 的特征向量，$A$的左奇异向量是$A^{T}A$ 的特征向量。
       SVD的核心就是找到一组A的行空间中正交基 ${\mathbf{v}}_1,{\mathbf{v}}_2,\dots {\mathbf{v}}_r$ ，使得:
A [ v 1    v 2 ⋯ v r ] = [ σ 1 u 1    σ 2 u 2 ⋯ σ r u r ] = [ u 1    u 2 ⋯ u r ] [ σ 11 σ 22 ⋱ σ r r ]

$$\begin{aligned} A[\mathbf{v}_1 \; \mathbf{v}_2 \cdots \mathbf{v}_r] &= [\sigma_1\mathbf{u}_1 \; \sigma_2\mathbf{u}_2 \cdots \sigma_r\mathbf{u}_r] \\ & = [\mathbf{u}_1 \; \mathbf{u}_2 \cdots \mathbf{u}_r] \begin{bmatrix} \sigma_{11} & & & \\ & \sigma_{22} & & \\ & & \ddots &\\ & & &\sigma_{rr} \end{bmatrix} \end{aligned}$$ A[v1​v2​⋯vr​]​=[σ1​u1​σ2​u2​⋯σr​ur​]=[u1​u2​⋯ur​] ​σ11​​σ22​​⋱​σrr​​ ​​
其中 ${\mathbf{u}}_1,{\mathbf{u}}_2,\dots {\mathbf{u}}_r$ 是一组在A列空间中的正交基。写成矩阵形式： $AV=U\Sigma$ 。其中$V$和$U$分别是行空间和列空间中的基，通常情况下： $U\ne V$ 。但是当A是正定矩阵的时候，$V$和$U$是相同的。也就是说一组正交基在行空间和列空间中是一样的。所以，特征值分解或者叫对角分解实际上是SVD的一种特殊形式。

SVD求解线性最小二乘⁶

       用SVD求解：设 $A$ 的SVD为
A = U Σ V H = ( U 1 , U 2 ) ( Σ 1 0 ) V H A=U \Sigma V^H=\left(U_1, U_2\right)\left(

$$\begin{array}{c} \Sigma_1 \\ 0 \end{array}$$ \right) V^H A=UΣVH=(U1​,U2​)(Σ1​0​)VH
记 $V^H\mathbf{x}=\mathbf{y},U_1^H\mathbf{b}={\mathbf{b}}_1,(U_1为前n列构成的矩阵),U_2^H\mathbf{b}={\mathbf{b}}_2$ 则
$$\begin{array}{rl}& \Vert A\mathbf{x}-\mathbf{b}{\Vert }_2^2=\Vert U\Sigma \mathbf{y}-\mathbf{b}{\Vert }_2^2={\Vert \Sigma \mathbf{y}-U^H\mathbf{b}\Vert }_2^2={\Vert \left(\begin{array}{c}{\Sigma }_1\\ 0\end{array}\right)\mathbf{y}-\left(\begin{array}{l}{\mathbf{b}}_1\\ {\mathbf{b}}_2\end{array}\right)\Vert }_2^2\\ & ={\Vert {\Sigma }_1\mathbf{y}-{\mathbf{b}}_1\Vert }_2^2+{\Vert {\mathbf{b}}_2\Vert }_2^2\end{array}$$
当 $A,\mathbf{b}$ 给定，则 ${\Vert {\mathbf{b}}_2\Vert }_2^2={\Vert U_2^H\mathbf{b}\Vert }_2^2$ 是一个常数，为使 $\Vert A\mathbf{x}-\mathbf{b}{\Vert }_2^2=\min$ 只要取 ${\Sigma }_1\mathbf{y}={\mathbf{b}}_1\iff \mathbf{y}={\Sigma }_1^{-1}{\mathbf{b}}_1={\Sigma }_1^{-1}{U}_1^H\mathbf{b}$
这就得到线性方程组的最小二乘解
$$\mathbf{x}=V{\Sigma }_1^{-1}{U}_1^H\mathbf{b}=\sum _{k=1}^n\frac{{\mathbf{u}}_k^H\mathbf{b}}{{\sigma }_k}{\mathbf{v}}_k$$

SVD分解的数值稳定性⁷

       稳定性: 在实际使用计算机计算时，SVD分解在数值计算上是很稳定的，这是因为SVD分解是对 $A$ 进行一系列正交变换(旋转和平移就是正交变换)，而正交变换在实际计算中是很稳定的(这是因为正交矩阵的条件数(酉不变范数下)都是 1 )，例如：对于一个真实数据 $A$ 和扰动后的矩阵 $\hat{A}=A+E$ ，其中 $E$ 是误差(存储数据和对数据进行计算时都会使真实数据产生误差, 这一误差 可以在某种程度上减弱，但很难彻底消除，则有下面的式子成立:
$$\left|{\sigma }_i-{\hat{\sigma }}_i\right|\le \Vert \hat{A}-A{\Vert }_2=\Vert E{\Vert }_2$$
其中， ${\hat{\sigma }}_i$ 是 $\hat{A}$ 奇异值，也就是说，只要误差足够小我们算出的奇异值就会足够靠近真实矩阵的奇异值，而特征值一般是不具有的这么良好的稳定性，也就是说，即使 $\Vert E{\Vert }_2$ 很小， $\hat{A}$ 的特征值 有可能和 $A$ 的特征值相差很大，例如
A = ( 1 50000 0 1 ) , E = ( 0 0 1 0 − 6 0 ) , A ^ = ( 1 50000 1 0 − 6 1 ) A=\left(

$$\begin{array}{cc} 1 & 50000 \\ 0 & 1 \end{array}$$ \right), E=\left( $$\begin{array}{cc} 0 & 0 \\ 10^{-6} & 0 \end{array}$$ \right), \hat{A}=\left( $$\begin{array}{cc} 1 & 50000 \\ 10^{-6} & 1 \end{array}$$ \right) A=(10​500001​),E=(010−6​00​),A^=(110−6​500001​)
       不难看出 $A$ 的两个特征值都是 1 ，而我们只在 $A$ 的左下角的分量添加了一个很小的扰动 $10^{-6}$ ， 但是 $\hat{A}$ 的特征值是 $1.223606797749979$ 和 $0.776393202250021$ ，与 1 相差太大，根本达不到 一般情况下我们所要求的精度。但是
       $A$ 两个奇异值是: $50000.00002$ 和 $0.000020$
       $\hat{A}$ 两个奇异值是: $50000.00002$ 和 $0.000019$
几乎相差无几。也正是因为SVD分解的数值稳定性，才使得其有很多广泛的应用。

SVD分解的Eigen实现

       Eigen提供了两种计算SVD分解的方法，分别为Eigen::BDCSVD和Eigen::JacobiSVD。对于小矩阵(<16)，最好使用Eigen::JacobiSVD，但是对于较大的矩阵强烈建议使用BDCSVD，可以快几个数量级。

方法一：JacobiSVD

int main(int argc, char *argv[])
{
     MatrixXf m(3, 2);
     m << 0.68, 0.597, -0.211, 0.823, 0.566, -0.605;
     cout << "Here is the matrix m:" << endl
          << m << endl;
     JacobiSVD<MatrixXf> svd;
     svd.compute(m, ComputeThinU | ComputeThinV);
     cout << "Its singular values are:" << endl
          << svd.singularValues() << endl;
     cout << "Its left singular vectors are the columns of the thin U matrix:" << endl
          << svd.matrixU() << endl;
     cout << "Its right singular vectors are the columns of the thin V matrix:" << endl
          << svd.matrixV() << endl;
     Vector3f rhs(1, 0, 0);
     cout << "Now consider this rhs vector:" << endl
          << rhs << endl;
     cout << "A least-squares solution of m*x = rhs is:" << endl
          << svd.solve(rhs) << endl;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

结果分析：

Here is the matrix m:
  0.68  0.597
-0.211  0.823
 0.566 -0.605
Its singular values are:
 1.19173
0.898234
Its left singular vectors are the columns of the thin U matrix:
  0.388338   0.865677
  0.711313 -0.0636487
 -0.585856   0.496541
Its right singular vectors are the columns of the thin V matrix:
-0.182601  0.983187
 0.983187  0.182601
Now consider this rhs vector:
1
0
0
A least-squares solution of m*x = rhs is:
0.888048
0.496366
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

方法二：BDCSVD

int main(int argc, char *argv[])
{
     MatrixXf m(3, 2);
     m << 0.68, 0.597, -0.211, 0.823, 0.566, -0.605;
     cout << "Here is the matrix m:" << endl
          << m << endl;
     BDCSVD<MatrixXf> svd;
     svd.compute(m, ComputeThinU | ComputeThinV);
     cout << "Its singular values are:" << endl
          << svd.singularValues() << endl;
     cout << "Its left singular vectors are the columns of the thin U matrix:" << endl
          << svd.matrixU() << endl;
     cout << "Its right singular vectors are the columns of the thin V matrix:" << endl
          << svd.matrixV() << endl;
     Vector3f rhs(1, 0, 0);
     cout << "Now consider this rhs vector:" << endl
          << rhs << endl;
     cout << "A least-squares solution of m*x = rhs is:" << endl
          << svd.solve(rhs) << endl;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

结果分析：

Here is the matrix m:
  0.68  0.597
-0.211  0.823
 0.566 -0.605
Its singular values are:
 1.19173
0.898234
Its left singular vectors are the columns of the thin U matrix:
  0.388338   0.865677
  0.711313 -0.0636487
 -0.585856   0.496541
Its right singular vectors are the columns of the thin V matrix:
-0.182601  0.983187
 0.983187  0.182601
Now consider this rhs vector:
1
0
0
A least-squares solution of m*x = rhs is:
0.888048
0.496366
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

---

1. 矩阵分解 ↩︎ ↩︎ ↩︎

2. QR分解 ↩︎

3. 常见矩阵分解 ↩︎

4. SVD分解 ↩︎

5. MIT线性代数笔记3.5(SVD分解) ↩︎

6. svd求解线性最小二乘 ↩︎

7. SVD分解数值稳定性 ↩︎