当线性方程组 $\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.