矩阵分解及其Eigen实现

这篇具有很好参考价值的文章主要介绍了矩阵分解及其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=[1020111] ,给 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=[10200111020] ,如果这里用的是浮点数的话,那么 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[1020011020] , 把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=[1020110]=[1020111]=A如果用这组LU去计算问题,会出现条件数很小,但相对误差很大的情况。
注: 因为浮点数的固有问题, 在处理数值计算问题时候,一定要尽可能避免大数和小数相加减的问题。

LUP分解1

理论

       在高斯消元的时候,做一些行置换(pivoting)。
矩阵分解及其Eigen实现

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=P1LUQ1,其中, 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.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

方法二: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.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

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=PTLDLP,其中, 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 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

方法二: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 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

QR分解2,3

理论

定义

       一个矩阵 A ∈ R m × n , m ≥ n A \in \mathbb{R}^{m \times n}, m \geq n ARm×n,mn 可以被分解成 A = Q R A=Q R A=QR, 其中:

  1. Q ∈ R m × m Q \in \mathbb{R}^{m \times m} QRm×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
  2. R ^ ∈ R n × n \hat{R} \in \mathbb{R}^{n \times n} R^Rn×n 是上三角矩阵

正交矩阵的性质

  1. Q T Q = Q Q T = I Q^T Q=Q Q^T=I QTQ=QQT=I
  2. 左乘一个正交矩阵对欧式范数的结果不影响
    ∥ 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 Qv22=vTQTQv=vTv=v22

QR分解求解线性最小二乘

       对于超定线性最小二乘问题 A x ≈ b Ax \approx b Axb,目标函数为:
ϕ ( 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=bAx22= bQ[R^0]x 22= QT(bQ[R^0]x) 22= QTb[R^0]x 22
其中, Q ∈ R m × m Q \in \mathbb{R}^{m \times m} QRm×m Q b ∈ R m × 1 Qb \in \mathbb{R}^{m \times 1} QbRm×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]xRm×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 c1Rn c 2 ∈ R m − n c_2 \in \mathbb{R}^{m-n} c2Rmn。那么目标函数可转换为:
∥ 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= c1R^x 22+c222
       显然, ∥ c 2 ∥ 2 2 \left\| c_2 \right\|_2^2 c222 x x x无关,能优化的只有前半部分,使得 ∥ c 1 − R ^ x ∥ 2 2 \left\| c_1 - \hat{R}x\right\|_2^2 c1R^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 c222。这样处理之后就可以避免正规方程中的 ( 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.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

方法二: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.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

方法三: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.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

方法四: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.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

SVD分解4,5

定义

       对于任意一个矩阵 A ∈ C m × n A \in \mathbb{C}^{m \times n} ACm×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[v1v2vr]=[σ1u1σ2u2σrur]=[u1u2ur] σ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的一种特殊形式。
矩阵分解及其Eigen实现

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} Axb22=UΣyb22= ΣyUHb 22= (Σ10)y(b1b2) 22=Σ1yb122+b222
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 b222= U2Hb 22 是一个常数,为使 ∥ A x − b ∥ 2 2 = min ⁡ \|A \boldsymbol{x}-\boldsymbol{b}\|_2^2=\min Axb22=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=b1y=Σ11b1=Σ11U1Hb
这就得到线性方程组的最小二乘解
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Σ11U1Hb=k=1nσ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σ^iA^A2=E2
其中, σ ^ i \hat{\sigma}_i σ^i A ^ \hat{A} A^ 奇异值,也就是说,只要误差足够小我们算出的奇异值就会足够靠近真实矩阵的奇异值,而特征值一般是不具有的这么良好的稳定性,也就是说,即使 ∥ E ∥ 2 \|E\|_2 E2 很小, 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=(010600),A^=(1106500001)
       不难看出 A A A 的两个特征值都是 1 ,而我们只在 A A A 的左下角的分量添加了一个很小的扰动 1 0 − 6 10^{-6} 106 , 但是 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.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

方法二: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.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. QR分解 ↩︎

  3. 常见矩阵分解 ↩︎

  4. SVD分解 ↩︎

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

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

  7. SVD分解数值稳定性 ↩︎文章来源地址https://www.toymoban.com/news/detail-437571.html

到了这里,关于矩阵分解及其Eigen实现的文章就介绍完了。如果您还想了解更多内容,请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章,希望大家以后多多支持TOY模板网!

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处: 如若内容造成侵权/违法违规/事实不符,请点击违法举报进行投诉反馈,一经查实,立即删除!

领支付宝红包 赞助服务器费用

相关文章

  • 小波分解及其Python实现方法

    小波分解及其Python实现方法 摘要: 小波分解是一种信号处理方法,可以将信号分解为不同频率和时间分辨率的小波系数。在本文中,我们将详细介绍小波分解的原理和算法,并使用Python编程语言演示实现。我们将介绍小波函数的选择、离散小波变换(DWT)和连续小波变换(

    2024年02月12日
    浏览(31)
  • Matlab实现矩阵的QR分解和奇异值分解

    1.安装并运行matlab软件; 2.在命令窗口行输入需要进行QR分解的矩阵,并输入求秩及进行QR分解的函数,如下图; 3.点击回车键,则可得Q及R矩阵; 4.若要查看之前所输入的矩阵及所获得的相关变量,可从右侧工作区窗口查看; 5.单击需要查看的变量名,则相关变量会被显示在主窗口

    2024年02月16日
    浏览(67)
  • 车辆姿态表达:旋转矩阵、欧拉角、四元数的转换以及eigen、matlab、pathon方法实现

    旋转矩阵、欧拉角、四元数主要用于表示坐标系中的旋转关系,通过三者之间的转换可以减小一些算法的复杂度。 本文主要概述旋转矩阵、欧拉角、四元数的基本理论、三者之间的转换关系以及三者转换在eigen、matlab和pathon上的方法实现。 对于两个三维点 p1 、 p2 : p 1 ( x

    2023年04月11日
    浏览(44)
  • Python实现矩阵奇异值分解(SVD)

    Python实现矩阵奇异值分解(SVD) 矩阵奇异值分解(Singular Value Decomposition, SVD)是一种重要的矩阵分解方法,可以将一个矩阵分解成三个矩阵的乘积,即 A = U Σ V T A=USigma V^{T} A = U Σ

    2024年02月10日
    浏览(42)
  • 奇异值分解与矩阵逆:数值实现与优化

    奇异值分解(Singular Value Decomposition, SVD)和矩阵逆(Matrix Inverse)是线性代数和数值分析领域中非常重要的概念和方法。这两者在现实生活中的应用非常广泛,例如图像处理、信号处理、数据挖掘、机器学习等领域。在这篇文章中,我们将从以下几个方面进行深入的讨论: 背景介绍

    2024年04月12日
    浏览(31)
  • 【SIMULINK】simulink实现信号矩阵整合、求逆、转置、分解、向量矩阵相乘(非matlab)

    simulink实现信号矩阵,并实现分解 simulink实现信号矩阵求逆 simulink实现信号矩阵转置 simulink矩阵向量相乘

    2024年02月11日
    浏览(41)
  • 推荐系统 | 基础推荐模型 | 矩阵分解模型 | 隐语义模型 | PyTorch实现

    基础推荐模型——传送门 : 推荐系统 | 基础推荐模型 | 协同过滤 | UserCF与ItemCF的Python实现及优化 推荐系统 | 基础推荐模型 | 矩阵分解模型 | 隐语义模型 | PyTorch实现 推荐系统 | 基础推荐模型 | 逻辑回归模型 | LS-PLM | PyTorch实现 推荐系统 | 基础推荐模型 | 特征交叉 | FM | FFM |

    2023年04月09日
    浏览(50)
  • LightFM:一款开源推荐系统框架,可以轻松实现大规模矩阵分解,快速、高效地处理大型矩阵

    作者:禅与计算机程序设计艺术 LightFM 是由 Yelp 开发的一款开源推荐系统框架,可以轻松实现大规模矩阵分解。该项目基于 TensorFlow 和 Keras 框架,可以快速、高效地处理大型矩阵。它具有以下特点: 提供了一种简单的方法来训练矩阵分解模型,即通过定义项间的交互矩阵和用

    2024年02月10日
    浏览(45)
  • 理解非负矩阵和张量分解:快速算法的Matlab实现与优化实践

    第一部分:非负矩阵分解(Non-negative Matrix Factorization,NMF)的基本原理 非负矩阵分解(NMF)是一种广泛应用的线性代数技术,特别适用于大规模的数据集分析。其基本思想是将一个非负矩阵分解为两个低秩的非负矩阵的乘积,使得矩阵的内在结构得以暴露并利于进一步分析。

    2024年02月16日
    浏览(54)
  • 机器学习实战教程(四):从特征分解到协方差矩阵:详细剖析和实现PCA算法

    方差和标准差的原理和实例演示,请参考 方差 方差(Variance)是度量一组数据的分散程度。方差是各个样本与样本均值的差的平方和的均值: 标准差 标准差是数值分散的测量。 标准差的符号是 σ (希腊语字母 西格马,英语 sigma) 公式很简单:方差的平方根。 协方差 通俗

    2024年02月02日
    浏览(50)

觉得文章有用就打赏一下文章作者

支付宝扫一扫打赏

博客赞助

微信扫一扫打赏

请作者喝杯咖啡吧~博客赞助

支付宝扫一扫领取红包,优惠每天领

二维码1

领取红包

二维码2

领红包