主要是用来记录自己的学习过程,内容也主要来自于网上的各种资料,然后自己总结而来,参考的资料都以注明,感谢这些作者的分享。如果内容有误,请大家指点。
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.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=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 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 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.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}
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.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
-
矩阵分解 ↩︎ ↩︎ ↩︎
-
QR分解 ↩︎
-
常见矩阵分解 ↩︎
-
SVD分解 ↩︎
-
MIT线性代数笔记3.5(SVD分解) ↩︎
-
svd求解线性最小二乘 ↩︎文章来源:https://www.toymoban.com/news/detail-437571.html
-
SVD分解数值稳定性 ↩︎文章来源地址https://www.toymoban.com/news/detail-437571.html
到了这里,关于矩阵分解及其Eigen实现的文章就介绍完了。如果您还想了解更多内容,请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章,希望大家以后多多支持TOY模板网!