|||
当线性方程组 $\textbfsymbol{Ax} = \textbfsymbol{b}$ 的系数矩阵 $\textbfsymbol{A}$ 对称正定时,直接三角分解法也可以简化,得到所谓的“平方根法”。
例如,线性回归模型 $\boldsymbol{Y} = \boldsymbol{X\beta}+\boldsymbol{\varepsilon}$ 的最小二乘解可以通过 $\boldsymbol{X}^T\boldsymbol{X}\boldsymbol{\beta} = \boldsymbol{X}^T\boldsymbol{Y}$ 求得,这里 $\boldsymbol{A} = \boldsymbol{X}^T\boldsymbol{X}$ 正定, $\boldsymbol{b} = \boldsymbol{X}^T\boldsymbol{Y}$,求回归参数 $\boldsymbol{\beta}$ 就是求对称正定方程组 $\boldsymbol{A\beta} = \boldsymbol{b}$ 的解.
方法的理论基础是下面的Cholesky(乔累斯基)分解定理。
Cholesky(乔累斯基)分解定理
设 $\textbfsymbol{A}\in \mathbb{R}^{n \times n}$ 对称正定,则存在一个实的非奇异下三角形矩阵 $\textbfsymbol{L}$ ,使得
$\textbfsymbol{A}=\textbfsymbol{L}\textbfsymbol{L}^T$
且当限定 $\textbfsymbol{L}$ 的对角元素为正时,此分解式是唯一的。
平方根法就是用 $\textbfsymbol{A}$ 的 $\textbfsymbol{L}\textbfsymbol{L}^T$ 分解来求解对称正定方程组 $\textbfsymbol{Ax} = \textbfsymbol{b}$ ,下面给出其求解步骤。
(1)求 $\textbfsymbol{A}$ 的 $\textbfsymbol{L}\textbfsymbol{L}^T$ 分解。
$\begin{bmatrix}
a_{11}&a_{12}&\cdots&a_{1n}\\
a_{21}&a_{22}&\cdots&a_{2n}\\
\cdots&\cdots&\cdots &\cdots\\
a_{n1}&a_{n2}&\cdots &a_{nn}\\
\end{bmatrix}
=\begin{bmatrix}
1_{11}\\
l_{21}&l_{22}\\
\cdots&\cdots&\ddots \\
l_{n1}&l_{n2}&\cdots &l_{nn}\\
\end{bmatrix}
\begin{bmatrix}
l_{11}&l_{21}&\cdots&l_{n1}\\
&l_{22}&\cdots&l_{n2}\\
& &\ddots&\vdots\\
& & &l_{nn}\\
\end{bmatrix}$
其中, $l_{ii}>0(i=1,2,\cdots,n)$ .
逐列(每列自上而下)考查 $\textbfsymbol{A}$ 的下三角部分,比较两边对应元素(用矩阵乘法规则),便可逐列确定 $\textbfsymbol{L}$ 的元素。
由 $\textbfsymbol{A}$ 的第 1 列元素 $a_{i1}$ 计算 $\textbfsymbol{L}$ 的第 1 列元素 $l_{i1}(i=1,2,\cdots,n.)$
$\left \{ \begin{array}{c} a_{11} = l_{11}\cdot l_{11}=l_{11}^2, \\ a_{i1} = l_{i1}\cdot l_{11}(i=2,\cdots,n) \end{array} \right. \Rightarrow \left \{ \begin{array}{c} l_{11} = \sqrt{a_{11}},\\ l_{i1} = a_{i1}/l_{11},(i=2,\cdots,n) \end{array} \right.$
一般地,由 $\textbfsymbol{A}$ 的第 $k$ 列元素 $a_{ik}$ 确定 $\textbfsymbol{L}$ 的第 $k$ 列元素 $l_{ik}(i=k,k+1,\cdots,n.)$
$\begin{array}{c} \left \{ \begin{array}{c} a_{kk} = l_{k1}^2+\cdots+ l_{kk}^2=\sum\limits_{r=1}^{k-1}l_{kr}^2+l_{kk}^2, \\ a_{ik} = l_{i1}l_{k1}+\cdots+l_{ik}l_{kk} =\sum\limits_{r=1}^{k-1}l_{ir}l_{kr}+l_{ik}l_{kk}, (i>k) \end{array} \right. \\ \Rightarrow \left \{ \begin{array}{c} l_{11} = \sqrt{a_{kk}-\sum\limits_{r=1}^{k-1}l_{kr}^2},\\ l_{i1} = (a_{ik}-\sum\limits_{r=1}^{k-1}l_{ir}l_{kr})/l_{11},(i>k) \end{array} \right. \end{array}$
(2) 求解下三角形方程组 $\textbfsymbol{L{}y} = \textbfsymbol{b}$ ,
$\left \{ \begin{array}{c} y_1 = b_1/l_{11},\\ y_i = (b_i-\sum\limits_{k=1}^{i-1}l_{ik}y_k)/l_{ii},(i=2,3,\cdots,n) \end{array} \right.$
(2) 求解上三角形方程组 $\textbfsymbol{L}^T{}\textbfsymbol{x} = \textbfsymbol{y}$ ,
$\left \{
\begin{array}{c}
x_n = y_n/l_{nn},\\
x_i = (y_i-\sum\limits_{k=i+1}^{n}l_{ki}x_k)/l_{ii},(i=n-1,n-2,\cdots,1)
\end{array}
\right.$
编写对应的函数
function [x,y,L] = Chol_Factor_Solve(A,b)
%% 矩阵A的 Cholesky 分解解线性方程组程序, 调用方法为
%% x = Chol_Factor_Solve(A,b)
%% 其中A为正定对称矩阵, L 为下三角阵, b 为 Ax = b 中常数项 b.
%% 矩阵A的 Cholesky 分解 A = L*L'
n = length(b); L = zeros(n);
for k = 1 : n
delta = A(k,k);
for j = 1 : k-1
delta = delta - L(k,j)^2;
end
if delta < 1e-15
error('% 矩阵不正定, 无法分解!');
end
L(k,k) = sqrt(delta);
for i = k+1 : n
L(i,k) = A(i,k);
for j = 1 : k-1
L(i,k) = L(i,k) - L(i,j)*L(k,j);
end
L(i,k) = L(i,k)/L(k,k);
end
end
% 解 Ly = b
y = b; y(1) = b(1)/L(1,1);
for i = 2:n
tempy = 0;
for k = 1:i-1
tempy = tempy + L(i,k)*y(k);
end
y(i) = (b(i)-tempy)/L(i,i);
end
% 解 L'x = y;
x = y; x(n) = y(n)/L(n,n);
for i = n-1:-1:1
tempx = 0;
for k = i+1:n
tempx = tempx + L(k,i)*x(k);
end
x(i) = (y(i)-tempx)/L(i,i);
end
例 用平方根法解对称正定方程组
$\begin{bmatrix} 16 & 4 & 8 \\ 4 & 5 & 4 \\ 8 & 4 & 22 \\ \end{bmatrix} \begin{bmatrix} x_1\\x_2\\x_3\\ \end{bmatrix}= \begin{bmatrix} 4 \\ 3 \\ 10 \\ \end{bmatrix}.$
>> A = [16 4 8;4 5 -4;8 -4 22];b = [-4;3;10];
>> [x,y,L] = Chol_Factor_Solve(A,b)
x =
-2.2500
4.0000
2.0000
y =
-1
2
6
L =
4 0 0
1 2 0
2 -3 3
(1) 算法稳定,平方根法是求解对称正定线性方程组的一种行之有效的方法。
(2) 可以进一步改进,避免开方运算。
参考文献:
姚仰新, 王福昌, 罗家洪, 庄楚强. 高等工程数学(第三版)[M].华南理工大学出版社, 广州, 2016-08.
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-5-20 13:06
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社