SVD方法求解线性方程组
本文详细推导如果利用奇异值分解得到
A
X
=
0
AX=0
A
X
=
0
及
A
X
=
b
AX=b
A
X
=
b
的解,推导如下。
齐次方程组求解AX=0
对于一个维度为
m
×
n
m\times n
m
×
n
的矩阵
A
A
A
来说，我们可以将其进行奇异值分解得到:
A
m
×
n
=
U
m
×
m
W
m
×
n
V
n
×
n
T
A_{m\times n} = U_{m \times m} W_{m\times n}V^T_{n\times n}
A
m
×
n
=
U
m
×
m
W
m
×
n
V
n
×
n
T
我们知道对于齐次线性方程组来说，如果位置量的维度为
n
n
n
，并且
r
a
n
k
(
A
)
=
m
>
n
rank(A)=m>n
r
ank
(
A
)
=
m
>
n
的话，约束多于位置量个数，一般这种情况我们会使用最小二乘法进行求解:
x
+
=
arg
min
∥
A
X
∥
2
x^{+}=\arg \min\Vert AX\Vert^2
x
+
=
ar
g
min
∥
A
X
∥
2
此时，
x
+
x^{+}
x
+
代表最小二乘法最优值。
基于拉格朗日法推导
我们将式(2)转换为一个带等式约束的最优问题：
x
+
=
arg
min
∥
A
X
∥
2
s
.
t
.
∥
X
∥
2
=
1
x^{+} = \arg\min\Vert AX \Vert^2 \\ s.t.~~ \Vert X \Vert^2 = 1
x
+
=
ar
g
min
∥
A
X
∥
2
s
.
t
.
∥
X
∥
2
=
1
我们知道
∀
k
\forall k
∀
k
那么
k
A
X
=
0
kAX=0
k
A
X
=
0
,所以我们可以这么限制。
针对式(3)使用拉格朗日乘子法:
L
(
X
,
λ
)
=
∥
A
X
∥
2
+
λ
(
1
−
∥
X
∥
2
)
=
A
X
X
T
A
T
+
λ
(
1
−
X
T
X
)
L(X, \lambda) = \Vert AX \Vert^2 + \lambda(1 – \Vert X \Vert^2)\\ = AXX^TA^T + \lambda(1 – X^TX)
L
(
X
,
λ
)
=
∥
A
X
∥
2
+
λ
(
1
−
∥
X
∥
2
)
=
A
X
X
T
A
T
+
λ
(
1
−
X
T
X
)
让朗格朗日函数
L
L
L
对
λ
\lambda
λ
与
x
x
x
求导得:
∂
L
∂
λ
=
1
−
X
T
X
=
0
∂
L
∂
X
=
2
A
T
A
X
−
2
λ
X
=
0
\frac{\partial L}{\partial \lambda} = 1-X^TX = 0\\ \frac{\partial L}{\partial X} = 2A^TAX-2\lambda X = 0
∂
λ
∂
L
=
1
−
X
T
X
=
0
∂
X
∂
L
=
2
A
T
A
X
−
2
λ
X
=
0
于是得到:
A
T
A
=
λ
A^TA = \lambda
A
T
A
=
λ
于是，式(2)可以化简为:
x
+
=
arg
min
∥
A
X
∥
2
=
arg
min
X
T
A
T
A
X
=
arg
min
X
T
λ
X
x^{+}=\arg \min\Vert AX\Vert^2\\ =\arg\min X^TA^TAX\\ =\arg\min X^T\lambda X
x
+
=
ar
g
min
∥
A
X
∥
2
=
ar
g
min
X
T
A
T
A
X
=
ar
g
min
X
T
λ
X
由式(7)，我们知道这里的
λ
\lambda
λ
其实是
A
T
A
A^TA
A
T
A
的特征值，也就是说这里的
λ
\lambda
λ
最小为
A
T
A
A^TA
A
T
A
最小特征值。带入式(8)，
x
+
x^{+}
x
+
对应着
A
T
A
A^TA
A
T
A
最小特征值对应的特征向量。对应着也就是矩阵
A
A
A
对应的右奇异矩阵最后一行特征向量。
基于SVD推导
我们直接将式(1)带入式(2)，得到:
x
+
=
arg
min
∥
U
W
V
T
X
∥
2
=
arg
min
∥
W
V
T
X
∥
2
x^{+}=\arg \min\Vert UWV^TX\Vert^2\\ =\arg\min\Vert WV^TX \Vert^2\\
x
+
=
ar
g
min
∥
U
W
V
T
X
∥
2
=
ar
g
min
∥
W
V
T
X
∥
2
我们令
Y
=
V
T
X
Y = V^TX
Y
=
V
T
X
，于是我们要最小化
∥
W
Y
∥
2
\Vert WY\Vert^2
∥
WY
∥
2
，其中
W
=
d
i
g
(
σ
1
,
σ
2
,
⋯
,
σ
n
)
W=dig(\sigma_1,\sigma_2,\cdots, \sigma_n)
W
=
d
i
g
(
σ
1
,
σ
2
,
⋯
,
σ
n
)
，于是上式等价于:
x
+
=
arg
min
∥
W
Y
∥
2
=
arg
min
Y
T
W
T
W
Y
x^{+}=\arg\min\Vert WY \Vert^2\\ =\arg\min Y^TW^TWY
x
+
=
ar
g
min
∥
WY
∥
2
=
ar
g
min
Y
T
W
T
WY
其中:
Y
T
W
T
W
Y
=
σ
1
2
y
1
2
+
σ
2
2
y
2
2
+
⋯
+
σ
n
2
y
n
2
Y^TW^TWY=\sigma_1^2y_1^2+\sigma_2^2y_2^2+\cdots+\sigma_n^2y_n^2
Y
T
W
T
WY
=
σ
1
2
y
1
2
+
σ
2
2
y
2
2
+
⋯
+
σ
n
2
y
n
2
，我们知道奇异值是按大小排序得，也就是最后那个奇异值最小，那么当
Y
=
(
0
,
0
,
⋯
,
1
)
Y = (0,0,\cdots,1)
Y
=
(
0
,
0
,
⋯
,
1
)
时取最小值。
我们反代回去得：
X
=
V
T
[
0
0
⋮
1
]
X = V^T\begin{bmatrix} 0\\0\\\vdots\\1 \end{bmatrix}
X
=
V
T
⎣
⎡
0
0
⋮
1
⎦
⎤
可知当
X
X
X
取右奇异值最后一行时，式(2)取最小值。
非齐次方程求解AX=b
这里使用SVD进行求解，与齐次方程类似，我们定义最小二乘法的代价函数:
x
+
=
arg
min
∥
A
X
−
b
∥
2
=
arg
min
∥
U
W
V
T
X
−
b
∥
2
=
arg
min
∥
W
V
T
X
−
U
T
b
∥
2
x^{+} = \arg\min \Vert AX -b\Vert^2\\ =\arg\min\Vert UWV^TX-b\Vert^2 \\ = \arg\min\Vert WV^TX-U^Tb\Vert ^2
x
+
=
ar
g
min
∥
A
X
−
b
∥
2
=
ar
g
min
∥
U
W
V
T
X
−
b
∥
2
=
ar
g
min
∥
W
V
T
X
−
U
T
b
∥
2
我们定义
Y
=
V
T
X
Y = V^TX
Y
=
V
T
X
以及
b
ˉ
=
U
T
b
\bar{b}=U^Tb
b
ˉ
=
U
T
b
，得：
x
+
=
arg
min
∥
W
Y
−
b
ˉ
∥
2
x^{+} = \arg \min \Vert WY-\bar{b}\Vert^2
x
+
=
ar
g
min
∥
WY
−
b
ˉ
∥
2
要使得上式最小：
W
m
×
n
Y
n
×
1
=
b
ˉ
m
×
1
W_{m\times n}Y_{n\times 1}= \bar{b}_{m\times 1}
W
m
×
n
Y
n
×
1
=
b
ˉ
m
×
1
对于一般非齐次线性方程
A
m
×
n
x
n
×
1
=
b
m
×
1
A_{m\times n}x_{n\times 1}=b_{m\times 1}
A
m
×
n
x
n
×
1
=
b
m
×
1
，一般我们遇到得问题通常
m
>
n
m > n
m
>
n
，于是我们可以得到最优解:
x
+
=
V
m
×
m
W
m
×
n
−
1
U
n
×
n
T
b
m
×
1
x^{+} = V_{m\times m} W_{m\times n}^{-1}U^T_{n\times n} b_{m\times 1}
x
+
=
V
m
×
m
W
m
×
n
−
1
U
n
×
n
T
b
m
×
1
另外，我们知道
x
+
x^{+}
x
+
得维度为
n
×
1
n\times 1
n
×
1
，如果
m
>
n
m >n
m
>
n
的话，我们只取到有效为的奇异值
σ
i
\sigma_i
σ
i
，其他的用1代替。即:
x
+
=
V
n
×
n
[
{
U
m
×
m
T
B
m
×
1
}
1
:
r
[
σ
1
σ
2
⋯
σ
r
⏞
r
]
r
×
1
1
…
1
⏞
n
−
r
]
n
×
1
T
x^{+} = V_{n\times n} \begin{bmatrix} \frac{ \{U_{m \times m}^TB_{m \times 1}\}_{1:r} } {\begin{bmatrix} \overbrace{ \sigma_1 ~~ \sigma_2 ~~ \cdots ~~ \sigma_r}^{r} \end{bmatrix}}_{r\times 1} &\overbrace{ 1 ~~\dots ~~ 1}^{n-r} \end{bmatrix}_{n\times 1}^T
x
+
=
V
n
×
n
⎣
⎡
[
σ
1
σ
2
⋯
σ
r
r
]
{
U
m
×
m
T
B
m
×
1
}
1
:
r
r
×
1
1
…
1
n
−
r
⎦
⎤
n
×
1
T
上述的
r
∈
[
1
,
n
]
r\in [1,n]
r
∈
[
1
,
n
]
。
MATLAB求解
% determine the effective rank r of A using singular values
r = 1;
while( r < size(A,2) & s(r+1) >= max(size(A))*eps*s(1) )
r = r+1;
end
d = U'*b;
x = V* ( [d(1:r)./s(1:r); zeros(n-r,1) ] );
c++代码求解
#include <bits/stdc++.h>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;
int main() {
srand(time(NULL));
Matrix3d A;
A << rand()%10, rand()%10, rand()%10,
rand()%10, rand()%10, rand()%10,
rand()%10, rand()%10, rand()%10;
cout << "A = " << endl << A << endl;
Vector3d B(1, 4, 6);
cout << "B = " << B.transpose() << endl;
JacobiSVD<MatrixXd> svd(A, ComputeThinU | ComputeThinV);
Matrix3d U = svd.matrixU();
Matrix3d V = svd.matrixV();
Matrix3d W = Matrix3d::Zero();
Vector3d w = svd.singularValues();
for(int i = 0; i < 3; ++i) W.row(i)[i] = w[i];
cout << "U = " << endl << U << endl;
cout << "W = " << endl << W << endl;
cout << "V = " << endl << V << endl;
Vector3d x = V*W.inverse()*U.transpose()*B;
cout << "x = " << x.transpose() << endl;
Vector3d b_ = A*x;
cout << "A*x = " << b_ << endl;
return 0;
}
输出:
A =
1 5 8
1 2 5
1 4 6
B = 1 4 6
U =
0.723064 0.26416 -0.638278
0.412049 -0.906546 0.0915979
0.554432 0.329233 0.764337
W =
13.1138 0 0
0 0.987561 0
0 0 0.231648
V =
0.128837 -0.317097 0.939601
0.507645 0.835026 0.212197
0.851879 -0.449645 -0.268555
x = 18 3 -4
A*x = 1
4
6