矩阵分解及其Eigen实现
主要是用来记录自己的学习过程，内容也主要来自于网上的各种资料，然后自己总结而来，参考的资料都以注明，感谢这些作者的分享。如果内容有误，请大家指点。
LU分解1
理论
定义
将矩阵等价为两个矩阵 L L L和 U U U的乘积 ，其中 L L L和 U U U分别是单位下三角矩阵和上三角矩阵，当 A A A的前rank(A)阶顺序主子式都不为0时，矩阵 A A A可以分解为
A = L U A=LU A=LU
其中 L L L是下三角矩阵， U U U是上三角矩阵。
LU分解求线性方程组
对于求解线性方程，将 A A A进行 L U LU LU分解，得到 L ( U x ) = b L(Ux)=b L(Ux)=b，令 y = U x y = Ux y=Ux，则有 L y = b Ly=b Ly=b，由于 L L L是下三角矩阵，可以方便求出 y y y，求出 y y y之后，在求解 U x = y Ux = y Ux=y，由于 U U U是上三角矩阵，这样就可以方便求出 x x x。
LU分解数值稳定值
考虑一个矩阵 A = [ 0 1 1 1 ] A=\left[\begin{array}{ll}0 & 1 \\ 1 & 1\end{array}\right] A=[0111] ，虽然A非奇异，且条件数很小 κ ( A ) ≈ 2.62 \kappa(A) \approx 2.62 κ(A)≈2.62 ，但是LU分解算法第一次执行就失败了，因为 A 0 , 0 = 0 A_{0,0}=0 A0,0=0 在分母上。
考虑另一个矩阵 A = [ 1 0 − 20 1 1 1 ] A=\left[\begin{array}{cc}10^{-20} & 1 \\ 1 & 1\end{array}\right] A=[10−20111] ，给 A 0 , 0 A_{0,0} A0,0 一个误差。可以计算得到 L = [ 1 0 1 0 20 1 ] L=\left[\begin{array}{cc}1 & 0 \\ 10^{20} & 1\end{array}\right] L=[1102001] ， U = [ 1 0 − 20 1 0 1 − 1 0 20 ] U=\left[\begin{array}{cc}10^{-20} & 1 \\ 0 & 1-10^{20}\end{array}\right] U=[10−20011−1020] ，如果这里用的是浮点数的话，那么 U ≈ [ 1 0 − 20 1 0 − 1 0 20 ] U \approx\left[\begin{array}{cc}10^{-20} & 1 \\ 0 & -10^{20}\end{array}\right] U≈[10−2001−1020] ， 把LU相乘得到 L U = [ 1 0 − 20 1 1 0 ] ≠ [ 1 0 − 20 1 1 1 ] = A L U=\left[\begin{array}{cc}10^{-20} & 1 \\ 1 & 0\end{array}\right] \neq\left[\begin{array}{cc}10^{-20} & 1 \\ 1 & 1\end{array}\right]=A LU=[10−20110]=[10−20111]=A 。如果用这组LU去计算问题，会出现条件数很小，但相对误差很大的情况。
注: 因为浮点数的固有问题, 在处理数值计算问题时候，一定要尽可能避免大数和小数相加减的问题。
LUP分解1
理论
在高斯消元的时候，做一些行置换（pivoting）。
LUP分解求线性方程组
- step1：LUP分解
- step2：求解 L y = P b Ly = Pb Ly=Pb
- step3：求解 U x = y Ux = y Ux=y
LUP分解数值稳定性
LUP算法比LU算法稳定
LU分解的Eigen实现
Eigen提供了两种LU分解的方法，分别为Eigen::FullPivLU和Eigen::PartialPivLU。其中，Eigen::PartialPivLU类提供了可逆方阵的LU分解，具有部分置换(partial pivoting)，将矩阵 A A A分解为 A = P L U A = PLU A=PLU， L L L为单位下三角矩阵， U U U为上三角矩阵， P P P为置换矩阵(permutation matrix)。但是，Eigen::FullPivLU类可以对任意矩阵进行LU分解，并具有完全置换(complete pivoting)，矩阵 A A A分解为 A = P − 1 L U Q − 1 A=P^{-1}LUQ^{-1} A=P−1LUQ−1，其中， L L L为单位下三角矩阵， U U U为上三角矩阵， P , Q P, Q 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;
}
结果分析：
初始矩阵A = 0.680375 0.59688 -0.329554
-0.211234 0.823295 0.5364590.566198 -0.604897 -0.444451
b = 0.10794
-0.04520590.257742
L = 1 0 00.72499 1 0
-0.734727 0.493089 1
U = 0.823295 -0.211234 0.5364590 0.833518 -0.7184820 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.5364590.566198 -0.604897 -0.444451
x = 0.60876
-0.2312820.51038
方法二：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;
}
结果分析：
A = 0.680375 0.59688 -0.329554
-0.211234 0.823295 0.5364590.566198 -0.604897 -0.444451
0.208598
P =
1 0 0
0 0 1
0 1 0
L = 1 0 00.832185 1 0
-0.310467 -0.915573 1
U = 0.680375 0.59688 -0.3295540 -1.10161 -0.17020 0 0.278313
重构后的系统矩阵A: 0.680375 0.59688 -0.329554
-0.211234 0.823295 0.5364590.566198 -0.604897 -0.444451
Cholesky分解1
理论
定义
是LU分解的一种特殊情况，当系统矩阵 A A A为对称正定矩阵时，此时可以得到 U = L T U = L^T U=LT，即
A = L L T A = LL^T A=LLT
Choleskyfen分解数值稳定性
不需要想普通的矩阵一样还需要置换操作，因此相比于LU，LUP分解数值更加稳定。
Cholesky分解Eigen实现
同样地，EIgen也提供了两种方法，分别为Eigen::LDLT以及Eigen::LLT。Eigen::LLT提供了对称正定矩阵的分解方法，使得 A = L L T A=LL^T A=LLT，其中 L L L为下三角矩阵。Eigen::LDLT是EIgen提供的更加鲁棒的Cholesky分解方法，只需要矩阵为半正定或者半负定即可，将 A A A分解为 A = P T L D L ∗ P A = P^TLDL^*P A=PTLDL∗P，其中， P P P为置换矩阵， L L L 是具有单位对角线的下三角矩阵， D D 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;
}
结果分析：
The matrix A is4 -1 2
-1 6 02 0 5
L = 2 0 0-0.5 2.39792 01 0.208514 1.9891
重构后的矩阵A: 4 -1 2
-1 6 02 0 5
方法二：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);
}
结果分析：
The matrix A is4 -1 2
-1 6 02 0 5
b = 0.680375
-0.2112340.566198
L = 1 0 00 1 0
-0.166667 0.4 1
D = 65
3.03333
x = 0.13803
-0.01220070.0580278
QR分解2,3
理论
定义
一个矩阵 A ∈ R m × n , m ≥ n A \in \mathbb{R}^{m \times n}, m \geq n A∈Rm×n,m≥n 可以被分解成 A = Q R A=Q R A=QR, 其中:
- Q ∈ R m × m Q \in \mathbb{R}^{m \times m} Q∈Rm×m 是正交矩阵
R ≡ [ R ^ 0 ] ∈ R m × n R \equiv\left[\begin{array}{c}\hat{R} \\ 0\end{array}\right] \in \mathbb{R}^{m \times n} R≡[R^0]∈Rm×n - R ^ ∈ R n × n \hat{R} \in \mathbb{R}^{n \times n} R^∈Rn×n 是上三角矩阵
正交矩阵的性质
- Q T Q = Q Q T = I Q^T Q=Q Q^T=I QTQ=QQT=I
- 左乘一个正交矩阵对欧式范数的结果不影响
∥ Q v ∥ 2 2 = v T Q T Q v = v T v = ∥ v ∥ 2 2 \|Q v\|_2^2=v^T Q^T Q v=v^T v=\|v\|_2^2 ∥Qv∥22=vTQTQv=vTv=∥v∥22
QR分解求解线性最小二乘
对于超定线性最小二乘问题 A x ≈ b Ax \approx b Ax≈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 ∈ R m × m Q \in \mathbb{R}^{m \times m} Q∈Rm×m， Q b ∈ R m × 1 Qb \in \mathbb{R}^{m \times 1} Qb∈Rm×1， [ R ^ 0 ] ∈ R m × n \left[\begin{array}{c} \hat{R} \\ 0 \end{array}\right] \in \mathbb{R}^{m \times n} [R^0]∈Rm×n， [ R ^ 0 ] x ∈ R m × 1 \left[\begin{array}{c} \hat{R} \\ 0 \end{array}\right]x \in \mathbb{R}^{m \times 1} [R^0]x∈Rm×1。如果把 Q T b Q^Tb QTb拆分成上下两部分，即 Q T b = [ c 1 c 2 ] Q^Tb=\left[ \begin{array}{c} c_1 \\ c_2 \end{array}\right] QTb=[c1c2]，其中 c 1 ∈ R n c_1 \in \mathbb{R}^n c1∈Rn， c 2 ∈ R m − n c_2 \in \mathbb{R}^{m-n} c2∈Rm−n。那么目标函数可转换为：
∥ r ( x ) ∥ 2 2 = ∥ c 1 − R ^ x ∥ 2 2 + ∥ c 2 ∥ 2 2 \left\| r(x) \right\|_2^2=\left\| c_1 - \hat{R}x\right\|_2^2 + \left\| c_2 \right\|_2^2 ∥r(x)∥22= c1−R^x 22+∥c2∥22
显然， ∥ c 2 ∥ 2 2 \left\| c_2 \right\|_2^2 ∥c2∥22与 x x x无关，能优化的只有前半部分，使得 ∥ c 1 − R ^ x ∥ 2 2 \left\| c_1 - \hat{R}x\right\|_2^2 c1−R^x 22等于0，即 R ^ x = c 1 \hat{R}x=c_1 R^x=c1，\left| r(x) \right|_2^2的最小值为 ∥ c 2 ∥ 2 2 \left\| c_2 \right\|_2^2 ∥c2∥22。这样处理之后就可以避免正规方程中的 ( A T A ) − 1 (A^TA)^{-1} (ATA)−1，避免了条件数变成 c o n d ( A T A ) = c a n d ( A ) 2 cond(A^TA)=cand(A)^2 cond(ATA)=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 | ColPivHouseholderQR | FullPivHouseholderQR | CompleteOrthogonalDecomposition |
|---|---|---|---|
| A = Q R A=QR A=QR | A P = Q R AP=QR AP=QR | P A P ′ = Q R PAP^{'}=QR PAP′=QR | A P = Q [ T 0 0 0 ] Z AP = Q[\begin{array}{cc} T & 0\\ 0 & 0 \end{array}]Z AP=Q[T000]Z |
| Q Q Q为正交矩阵， R R R为上三角矩阵 | Q Q Q为正交矩阵， R R R为上三角矩阵， P P P为置换矩阵 | Q Q Q为正交矩阵， R R R为上三角矩阵， P , P ′ P,P^{'} P,P′为置换矩阵 | Q , Z Q,Z Q,Z为正交矩阵， P P P为置换矩阵， T T T为上三角矩阵，大小为A的秩 |
| 最快，但是最不稳定 | 稳定性较好，比HouseholderQR慢，比FullPivHouseholderQR快 | 数值最稳定，但是最慢 |
方法一：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;
}
结果分析：
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.05177030.58236 -0.461852 0.555351 0.176282 -0.328724-0.173814 -0.746785 -0.00906669 -0.348021 0.539351x = -0.262592
-0.0749101-0.593202
方法二：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;
}
结果分析：
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.7262860.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.150080 0.860025 -0.01190550 0 0.4383410 0 00 0 0
x = -0.262592
-0.0749101-0.593202
方法三：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;
}
结果分析：
Q = 0.124977 -0.919134 0.334273 0.162086 -0.03954220.467087 0.131775 -0.232567 0.826756 -0.1638640.493559 -0.0702728 -0.202056 -0.160453 0.82758-0.629484 -0.274961 -0.623708 0.274145 0.2529410.35547 -0.239345 -0.635875 -0.435088 -0.471929
x = -0.262592
-0.0749099-0.593202
方法四：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;
}
结果分析：
= -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.7262860.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.150080 0.860025 -0.01190550 0 0.4383410 0 00 0 0
Z =
1 0 0
0 1 0
0 0 1
x = -0.262592
-0.0749101-0.593202
SVD分解4,5
定义
对于任意一个矩阵 A ∈ C m × n A \in \mathbb{C}^{m \times n} A∈Cm×n，可以分解为：
A = U Σ V H A = U \Sigma V^H A=UΣVH
其中， U U U是 m × m m \times m m×m酉矩阵， V V V是 N × N N \times N N×N酉矩阵， Σ \Sigma Σ是 m × n m \times n m×n对角矩阵，对角线上的元素称为矩阵 A A A的奇异值( A T A A^TA ATA或者 A T A A^TA ATA特征值开根号)，矩阵 U U U的列向量称为左奇异向量(left singular vector)，矩阵 V V V的列向量称为右奇异向量(right singular vector)。 A A A的左奇异向量是 A A T AA^T AAT 的特征向量， A A A的左奇异向量是 A T A A^TA ATA 的特征向量。
SVD的核心就是找到一组A的行空间中正交基 v 1 , v 2 , … v r \mathbf{v}_1, \mathbf{v}_2, \ldots \mathbf{v}_r v1,v2,…vr ，使得:
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[v1v2⋯vr]=[σ1u1σ2u2⋯σrur]=[u1u2⋯ur] σ11σ22⋱σrr
其中 u 1 , u 2 , … u r \mathbf{u}_1, \mathbf{u}_2, \ldots \mathbf{u}_r u1,u2,…ur 是一组在A列空间中的正交基。写成矩阵形式： A V = U Σ A V=U \Sigma AV=UΣ 。其中 V V V和 U U U分别是行空间和列空间中的基，通常情况下： U ≠ V U \neq V U=V 。但是当A是正定矩阵的时候， V V V和 U U U是相同的。也就是说一组正交基在行空间和列空间中是一样的。所以，特征值分解或者叫对角分解实际上是SVD的一种特殊形式。
SVD求解线性最小二乘6
用SVD求解：设 A A 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)(Σ10)VH
记 V H x = y , U 1 H b = b 1 , ( U 1 为前 n 列构成的矩阵 ) , U 2 H b = b 2 V^H \boldsymbol{x}=\boldsymbol{y}, U_1^H \boldsymbol{b}=\boldsymbol{b}_1,(U_1为前n列构成的矩阵) ,U_2^H \boldsymbol{b}=\boldsymbol{b}_2 VHx=y,U1Hb=b1,(U1为前n列构成的矩阵),U2Hb=b2 则
∥ A x − b ∥ 2 2 = ∥ U Σ y − b ∥ 2 2 = ∥ Σ y − U H b ∥ 2 2 = ∥ ( Σ 1 0 ) y − ( b 1 b 2 ) ∥ 2 2 = ∥ Σ 1 y − b 1 ∥ 2 2 + ∥ b 2 ∥ 2 2 \begin{aligned} & \|A \boldsymbol{x}-\boldsymbol{b}\|_2^2=\|U \Sigma \boldsymbol{y}-\boldsymbol{b}\|_2^2=\left\|\Sigma \boldsymbol{y}-U^H \boldsymbol{b}\right\|_2^2=\left\|\left(\begin{array}{c} \Sigma_1 \\ 0 \end{array}\right) \boldsymbol{y}-\left(\begin{array}{l} \boldsymbol{b}_1 \\ \boldsymbol{b}_2 \end{array}\right)\right\|_2^2 \\ & =\left\|\Sigma_1 \boldsymbol{y}-\boldsymbol{b}_1\right\|_2^2+\left\|\boldsymbol{b}_2\right\|_2^2 \end{aligned} ∥Ax−b∥22=∥UΣy−b∥22= Σy−UHb 22= (Σ10)y−(b1b2) 22=∥Σ1y−b1∥22+∥b2∥22
当 A , b A, \boldsymbol{b} A,b 给定，则 ∥ b 2 ∥ 2 2 = ∥ U 2 H b ∥ 2 2 \left\|\boldsymbol{b}_2\right\|_2^2=\left\|U_2^H \boldsymbol{b}\right\|_2^2 ∥b2∥22= U2Hb 22 是一个常数，为使 ∥ A x − b ∥ 2 2 = min \|A \boldsymbol{x}-\boldsymbol{b}\|_2^2=\min ∥Ax−b∥22=min 只要取 Σ 1 y = b 1 ⇔ y = Σ 1 − 1 b 1 = Σ 1 − 1 U 1 H b \Sigma_1 \boldsymbol{y}=\boldsymbol{b}_1 \Leftrightarrow \boldsymbol{y}=\Sigma_1^{-1} \boldsymbol{b}_1=\Sigma_1^{-1} U_1^H \boldsymbol{b} Σ1y=b1⇔y=Σ1−1b1=Σ1−1U1Hb
这就得到线性方程组的最小二乘解
x = V Σ 1 − 1 U 1 H b = ∑ k = 1 n u k H b σ k v k \boldsymbol{x}=V \Sigma_1^{-1} U_1^H \boldsymbol{b}=\sum_{k=1}^n \frac{\boldsymbol{u}_k^H \boldsymbol{b}}{\sigma_k} \boldsymbol{v}_k x=VΣ1−1U1Hb=k=1∑nσkukHbvk
SVD分解的数值稳定性7
稳定性: 在实际使用计算机计算时，SVD分解在数值计算上是很稳定的，这是因为SVD分解是对 A A A 进行一系列正交变换(旋转和平移就是正交变换)，而正交变换在实际计算中是很稳定的(这是因为正交矩阵的条件数(酉不变范数下)都是 1 )，例如：对于一个真实数据 A A A 和扰动后的矩阵 A ^ = A + E \hat{A}=A+E A^=A+E ，其中 E E E 是误差(存储数据和对数据进行计算时都会使真实数据产生误差, 这一误差 可以在某种程度上减弱，但很难彻底消除，则有下面的式子成立:
∣ σ i − σ ^ i ∣ ≤ ∥ A ^ − A ∥ 2 = ∥ E ∥ 2 \left|\sigma_i-\hat{\sigma}_i\right| \leq\|\hat{A}-A\|_2=\|E\|_2 ∣σi−σ^i∣≤∥A^−A∥2=∥E∥2
其中， σ ^ i \hat{\sigma}_i σ^i 是 A ^ \hat{A} A^ 奇异值，也就是说，只要误差足够小我们算出的奇异值就会足够靠近真实矩阵的奇异值，而特征值一般是不具有的这么良好的稳定性，也就是说，即使 ∥ E ∥ 2 \|E\|_2 ∥E∥2 很小， A ^ \hat{A} A^ 的特征值 有可能和 A 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=(10500001),E=(010−600),A^=(110−6500001)
不难看出 A A A 的两个特征值都是 1 ，而我们只在 A A A 的左下角的分量添加了一个很小的扰动 1 0 − 6 10^{-6} 10−6 ， 但是 A ^ \hat{A} A^ 的特征值是 1.223606797749979 1.223606797749979 1.223606797749979 和 0.776393202250021 0.776393202250021 0.776393202250021 ，与 1 相差太大，根本达不到 一般情况下我们所要求的精度。但是
A A A 两个奇异值是: 50000.00002 50000.00002 50000.00002 和 0.000020 0.000020 0.000020
A ^ \hat{A} A^ 两个奇异值是: 50000.00002 50000.00002 50000.00002 和 0.000019 0.000019 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;
}
结果分析：
Here is the matrix m:0.68 0.597
-0.211 0.8230.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.8656770.711313 -0.0636487-0.585856 0.496541
Its right singular vectors are the columns of the thin V matrix:
-0.182601 0.9831870.983187 0.182601
Now consider this rhs vector:
1
0
0
A least-squares solution of m*x = rhs is:
0.888048
0.496366
方法二：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;
}
结果分析：
Here is the matrix m:0.68 0.597
-0.211 0.8230.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.8656770.711313 -0.0636487-0.585856 0.496541
Its right singular vectors are the columns of the thin V matrix:
-0.182601 0.9831870.983187 0.182601
Now consider this rhs vector:
1
0
0
A least-squares solution of m*x = rhs is:
0.888048
0.496366
矩阵分解 ↩︎ ↩︎ ↩︎
QR分解 ↩︎
常见矩阵分解 ↩︎
SVD分解 ↩︎
MIT线性代数笔记3.5(SVD分解) ↩︎
svd求解线性最小二乘 ↩︎
SVD分解数值稳定性 ↩︎
矩阵分解及其Eigen实现相关推荐
- Eigen密集矩阵求解 1 - 线性代数及矩阵分解
简介 这里介绍线性系统的解析,如何进行各种分解计算,如LU,QR,SVD,特征值分解等. 简单线性求解 在一个线性系统,常如下表示,其中A,b分别是一个矩阵,需要求x: Ax=bAx \:= \: b ...
- tf 如何进行svd_Tensorflow快餐教程(6) - 矩阵分解
摘要: 特征分解,奇异值分解,Moore-Penrose广义逆 矩阵分解 特征向量和特征值 我们在<线性代数>课学过方阵的特征向量和特征值. 定义:设A∈Fn×n是n阶方阵.如果存在非零向 ...
- Tensorflow快餐教程(6) - 矩阵分解
摘要: 特征分解,奇异值分解,Moore-Penrose广义逆 矩阵分解 特征向量和特征值 我们在<线性代数>课学过方阵的特征向量和特征值. 定义:设A∈Fn×nA∈Fn×n是n阶方阵.如 ...
- 视觉SLAM中的数学——解方程AX=b与矩阵分解：奇异值分解（SVD分解） 特征值分解 QR分解 三角分解 LLT分解
前言 本博客主要介绍在SLAM问题中常常出现的一些线性代数相关的知识,重点是如何采用矩阵分解的方法,求解线性方程组AX=B.主要参考了<计算机视觉--算法与应用>附录A以及Eigen库的方 ...
- java 矩阵分解_计算方法（三）矩阵分解1-正交分解(QR分解)
正交分解 矩阵的正交分解又称为QR分解,是将矩阵分解为一个正交矩阵Q和一个上三角矩阵的乘积的形式. 任意实数方阵A,都能被分解为 .这里的Q为正交单位阵,即 R是一个上三角矩阵.这种分解被称为QR分解 ...
- 【机器学习的数学基础】（六）矩阵分解(Matrix Decomposition)(上)
文章目录 4 矩阵分解(Matrix Decompositions)(上) 4.1 行列式与迹 4.2 特征值和特征向量 4 矩阵分解(Matrix Decompositions)(上) 在第2章和第 ...
- 特征值分解（Eigen Value Decomposition，EVD）、奇异值分解（Singular Value Decomposition，SVD）原理、公式推导及应用
1 正交矩阵&正交变换 正交变换是保持图形形状和大小不变的几何变换,包含旋转.平移.轴对称及这些变换的复合形式,正交变换可以保持向量的长度和向量之间的角度不变.特别的,标准正交基经正交变换后仍 ...
- 计算方法（三）矩阵分解1-正交分解(QR分解)
为什么80%的码农都做不了架构师?>>> 正交分解 矩阵的正交分解又称为QR分解,是将矩阵分解为一个正交矩阵Q和一个上三角矩阵的乘积的形式. 任意实数方阵A,都能被分解为 .这 ...
- 矩阵奇异值分解特征值分解_推荐系统中的奇异值分解与矩阵分解
矩阵奇异值分解特征值分解 Recently, after watching the Recommender Systems class of Prof. Andrew Ng's Machine Lea ...
最新文章
- JDBC编程的事务处理
- Ajax的get、post和ajax提交
- 几个常见的 slice 错误
- 腾讯高性能分布式路由技术，亮相亚太网络研讨会APNet
- 使用Azure云原生构建博客是怎样一种体验？（上篇）
- Android环境的安装遇到的问题
- 阶段项目:学生信息管理系统数据库设计
- 产品经理学技术之数据结构
- 构建高性能ASP.NET站点
- QtCreator格式化代码---Beautifier插件使用方式
- 利用图片的 onerror 事件载入默认图片
- OC---Math公式
- 银行卡四要素API 方便好用
- 乐队设备--反馈抑制器学习笔记
- 【观察】西部数据：再定义分层存储架构，赋能数据中心新基建
- 哈尔滨工业大学计算机系统大作业2022春
- 前端-获取treegrid的选中数据
- 克隆好的CentOS6虚拟机如何联网，解决报错Device eth0 does not seem to be present, delaying initialization
- Google Earth Engine（GEE）——
- 7.动态绘制(Jig)