算法:矩阵 LUP 分解
本文着笔于矩阵 LUP 分解算法,以及利用矩阵的 LUP 分解来解线性方程组、求矩阵对应行列式的值、求逆矩阵。
对于矩阵的定义代码如下:
struct Matrix
{
double dat[MAX_N][MAX_N],det,LUdat[MAX_N][MAX_N],rev[MAX_N][MAX_N];
int dgr,pi[MAX_N];
bool LUP_Decomposition();
VectorColumn LUP_Solve(VectorColumn);
void GetDeterminant();
void GetReverse();
};
矩阵的 LUP 分解
矩阵 A 的 LUP 分解即为找到下三角矩阵 L、上三角矩阵 U、行置换矩阵 P 使得 PA = LU ,实现方法为高斯消元。
E = I = ( 1 0 ⋯ 0 0 1 ⋯ 0 ⋮ ⋮ ⋱ ⋮ 0 0 ⋯ 1 ) , e i j = { 1 i = j 0 i ≠ j \begin{matrix} E \end{matrix} = \begin{matrix} I \end{matrix} = \begin{pmatrix} 1 & 0 & \cdots & 0 \\ 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1 \end{pmatrix}, e_{ij} = \begin{cases} 1 & i = j \\0 & i \neq j \end{cases} E=I=⎝⎜⎜⎜⎛10⋮001⋮0⋯⋯⋱⋯00⋮1⎠⎟⎟⎟⎞,eij={10i=ji̸=j
初始矩阵 P = E。对于矩阵 A,我们可以将它分解成这样两个矩阵的乘积(这里直接借用了 Introduction to Algorithms 里的描述)
A
=
(
a
11
ω
T
ν
A
′
)
=
(
1
0
ν
/
a
11
I
n
−
1
)
(
a
11
ω
T
0
A
′
−
ν
ω
T
/
a
11
)
\begin{matrix} A \end{matrix} = \begin{pmatrix} a_{11} & \omega^T \\ \nu & A' \end{pmatrix} = \begin{pmatrix} 1 & 0 \\ \nu / a_{11} & I_{n-1} \end{pmatrix} \begin{pmatrix} a_{11} & \omega^T \\ 0 & A' - \nu \omega^T / a_{11} \end{pmatrix}
A=(a11νωTA′)=(1ν/a110In−1)(a110ωTA′−νωT/a11)
这样,将这个规模为 n 的问题分解为一个规模为 n - 1 的子问题,迭代之即可求得矩阵 PA 的 LU 分解。最后得到的矩阵 U 相当于高斯消元后的结果,因为求矩阵 U 的实质是用当前规模矩阵的第一行的相应倍数减其它行,将其它行首列元素消为 0。
为了避免除数为零或过小引起精度问题,在迭代到矩阵第 i 行时,我们可以找到 j 使得 ∣ A j , i ∣ \vert A_{j, i} \vert ∣Aj,i∣ 尽可能大,其中 i ≤ j ≤ n i \leq j \leq n i≤j≤n。当 i 不等于 j 时,交换矩阵 PA 的第 i 行与第 j 行,相当于矩阵 A 不变,矩阵 P 的第 i 行与第 j 行交换。
对于矩阵的 LUP 分解代码如下:
#define ABS(x) \
({ \
typeof(x) __tmp_x=(x); \
__tmp_x<0?-__tmp_x:__tmp_x; \
})
#define EXCHANGE(a,b) \
{ \
typeof(a) __tmp_a=(a); \
a=b,b=__tmp_a; \
}
bool Matrix::LUP_Decomposition()
{
int p;
for(int i=1;i<=dgr;i++)
{
pi[i]=i;
for(int j=1;j<=dgr;j++)
LUdat[i][j]=dat[i][j];
}
det=1;
for(int i=1;i<=dgr;i++)
{
p=i;
for(int j=i+1;j<=dgr;j++)
if(ABS(LUdat[j][i])>ABS(LUdat[p][i]))
p=j;
if(LUdat[p][i]==0)
return false;
if(p!=i)
{
det=-det;
EXCHANGE(pi[i],pi[p]);
}
for(int j=1;j<=dgr;j++)
EXCHANGE(LUdat[i][j],LUdat[p][j]);
for(int j=i+1;j<=dgr;j++)
{
LUdat[j][i]/=LUdat[i][i];
for(int k=i+1;k<=dgr;k++)
LUdat[j][k]-=LUdat[j][i]*LUdat[i][k];
}
}
return true;
}
时间复杂度 O ( n 3 ) O(n^3) O(n3)。
矩阵求解线性方程组
对于线性方程组
{
a
11
x
1
+
a
12
x
2
+
⋯
+
a
1
n
x
n
=
b
1
a
21
x
1
+
a
22
x
2
+
⋯
+
a
2
n
x
n
=
b
2
⋮
⋮
⋮
⋮
a
n
1
x
1
+
a
n
2
x
2
+
⋯
+
a
n
n
x
n
=
b
n
\begin{cases} a_{11} x_1 + & a_{12} x_2 + \cdots + & a_{1n} x_n = & b_1 \\ a_{21} x_1 + & a_{22} x_2 + \cdots + & a_{2n} x_n = & b_2 \\ \vdots & \vdots & \vdots & \vdots \\ a_{n1} x_1 + & a_{n2} x_2 + \cdots + & a_{nn} x_n = & b_n \end{cases}
⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧a11x1+a21x1+⋮an1x1+a12x2+⋯+a22x2+⋯+⋮an2x2+⋯+a1nxn=a2nxn=⋮annxn=b1b2⋮bn
构建矩阵 A 和列向量 B,使得
A
=
(
a
11
a
12
⋯
a
1
n
a
21
a
22
⋯
a
2
n
⋮
⋮
⋱
⋮
a
n
1
a
n
2
⋯
a
n
n
)
,
B
=
(
b
1
b
2
⋮
b
n
)
\begin{matrix} A \end{matrix} = \begin{pmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} \end{pmatrix}, \begin{matrix} B \end{matrix} = \begin{pmatrix} b_1 \\ b_2 \\ \vdots \\ b_n \end{pmatrix}
A=⎝⎜⎜⎜⎛a11a21⋮an1a12a22⋮an2⋯⋯⋱⋯a1na2n⋮ann⎠⎟⎟⎟⎞,B=⎝⎜⎜⎜⎛b1b2⋮bn⎠⎟⎟⎟⎞
矩阵求解线性方程组即对于矩阵 A 和列向量 B,找到列向量 X,使得
A
X
=
B
\begin{matrix} A \end{matrix} \begin{matrix} X \end{matrix} = \begin{matrix} B \end{matrix}
AX=B
两边同乘 P,得
L
U
X
=
P
A
X
=
P
B
\begin{matrix} L \end{matrix} \begin{matrix} U \end{matrix} \begin{matrix} X \end{matrix} = \begin{matrix} P \end{matrix} \begin{matrix} A \end{matrix} \begin{matrix} X \end{matrix} = \begin{matrix} P \end{matrix} \begin{matrix} B \end{matrix}
LUX=PAX=PB
设列向量 Y = UX,得
L
Y
=
P
B
\begin{matrix} L \end{matrix} \begin{matrix} Y \end{matrix} = \begin{matrix} P \end{matrix} \begin{matrix} B \end{matrix}
LY=PB
设数组
π
1..
n
\pi_{1 .. n}
π1..n,其中
π
i
\pi_i
πi 表示矩阵 P 中第 i 行第
π
i
\pi_i
πi 列元素为 1,第 i 行其余元素均为 0,则 LY = PB 对应的方程组为
{
l
11
y
1
=
b
π
1
l
21
y
1
+
l
22
y
2
=
b
π
2
⋮
⋮
⋮
l
n
1
y
1
+
l
n
2
y
2
+
⋯
+
l
n
n
y
n
=
b
π
n
\begin{cases} l_{11} y_1 & & & = & b_{\pi_1} \\ l_{21} y_1 + & l_{22} y_2 & & = & b_{\pi_2} \\ \vdots & \vdots & & & \vdots \\ l_{n1} y_1 + & l_{n2} y_2 + \cdots + & l_{nn} y_n & = & b_{\pi_n} \end{cases}
⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧l11y1l21y1+⋮ln1y1+l22y2⋮ln2y2+⋯+lnnyn===bπ1bπ2⋮bπn
由定义,矩阵 L 中对角线元素均为 1。从上往下依次解出
y
1..
n
y_{1 .. n}
y1..n,这样解出 Y 后,利用方程 UX = Y,解出 X 即可。对应方程组为
{
u
11
x
1
+
u
12
x
2
+
⋯
+
u
1
n
x
n
=
y
1
u
22
x
2
+
⋯
+
u
2
n
x
n
=
y
2
⋮
⋮
u
n
n
x
n
=
y
n
\begin{cases} u_{11} x_1 + & u_{12} x_2 + \cdots + & u_{1n} x_n & = & y_1 \\ & u_{22} x_2 + \cdots + & u_{2n} x_n & = & y_2 \\ & & \vdots & & \vdots \\ & & u_{nn} x_n & = & y_n \end{cases}
⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧u11x1+u12x2+⋯+u22x2+⋯+u1nxnu2nxn⋮unnxn===y1y2⋮yn
从下往上依次解出 x n . . 1 x_{n .. 1} xn..1。在具体实现过程中,由于结果中 Y 无需保留,列向量 X 和 Y 可共用存储空间。
对于求解线性方程组的代码如下:
struct VectorColumn
{
double dat[MAX_N];
int raw;
};
VectorColumn Matrix::LUP_Solve(VectorColumn b)
{
VectorColumn x;
x.raw=dgr;
for(int i=1;i<=dgr;i++)
{
x.dat[i]=b.dat[pi[i]];
for(int j=1;j<i;j++)
x.dat[i]-=LUdat[i][j]*x.dat[j];
}
for(int i=dgr;i>=1;i--)
{
for(int j=dgr;j>i;j--)
x.dat[i]-=LUdat[i][j]*x.dat[j];
x.dat[i]/=LUdat[i][i];
}
return x;
}
时间复杂度 O ( n 2 ) O(n^2) O(n2)。(若算上矩阵 LUP 分解,总时间复杂度 O ( n 3 ) O(n^3) O(n3))
矩阵的行列式值求解
由行列式性质可知
∣
P
∣
∣
A
∣
=
∣
L
∣
∣
U
∣
\begin{vmatrix} P \end{vmatrix} \begin{vmatrix} A \end{vmatrix} = \begin{vmatrix} L \end{vmatrix} \begin{vmatrix} U \end{vmatrix}
∣∣P∣∣∣∣A∣∣=∣∣L∣∣∣∣U∣∣
其中,L 为对角线元素均为 1 的下三角矩阵,故其行列式值为 1;U 为上三角矩阵,其行列式值为对角线元素之积;P 初始行列式值为 1,每对 P 进行一次行交换,相当于行列式值取其相反数。
在具体实现过程中,行列式值的初值为 1 或 -1 在 LUP 分解过程中决定,此部分代码见上文。
对于矩阵的行列式值求解代码如下:
void Matrix::GetDeterminant()
{
for(int i=1;i<=dgr;i++)
det*=LUdat[i][i];
}
时间复杂度 O ( n ) O(n) O(n)。(若算上矩阵 LUP 分解,总时间复杂度 O ( n 3 ) O(n^3) O(n3))
矩阵求逆
求矩阵 A A A 的逆矩阵,即对于矩阵 A A A 和单位矩阵 E E E,找到矩阵 A − 1 A^{-1} A−1,使得 A A − 1 = E A A^{-1} = E AA−1=E。
对于矩阵
A
−
1
A^{-1}
A−1 和矩阵
E
E
E 的任一列,设列数为
i
i
i,构建列向量
A
c
o
l
(
i
)
−
1
=
(
a
1
i
a
2
i
⋮
a
n
i
)
,
E
c
o
l
(
i
)
=
(
e
1
i
e
2
i
⋮
e
n
i
)
\begin{matrix} A_{col(i)}^{-1} \end{matrix} = \begin{pmatrix} a_{1i} \\ a_{2i} \\ \vdots \\ a_{ni} \end{pmatrix}, \begin{matrix} E_{col(i)} \end{matrix} = \begin{pmatrix} e_{1i} \\ e_{2i} \\ \vdots \\ e_{ni} \end{pmatrix}
Acol(i)−1=⎝⎜⎜⎜⎛a1ia2i⋮ani⎠⎟⎟⎟⎞,Ecol(i)=⎝⎜⎜⎜⎛e1ie2i⋮eni⎠⎟⎟⎟⎞
则有
A
A
c
o
l
(
i
)
−
1
=
E
c
o
l
(
i
)
\begin{matrix} A \end{matrix} \begin{matrix} A_{col(i)}^{-1} \end{matrix} = \begin{matrix} E_{col(i)} \end{matrix}
AAcol(i)−1=Ecol(i)
套用矩阵求解线性方程组的方法即可。
对于矩阵求逆的代码如下:
void Matrix::GetReverse()
{
VectorColumn curcol,revcol;
curcol.raw=dgr;
for(int i=1;i<=dgr;i++)
curcol.dat[i]=0;
for(int i=1;i<=dgr;i++)
{
curcol.dat[i-1]=0,curcol.dat[i]=1;
revcol=LUP_Solve(curcol);
for(int j=1;j<=dgr;j++)
rev[j][i]=revcol.dat[j];
}
}
时间复杂度 O ( n 3 ) O(n^3) O(n3)。