一、思路轉變

A為對稱正定矩陣，$$Ax=b$$

求解向量$x$這個問題可以轉化為一個求$f(x)$極小值點的問題，為什么可以這樣：

$$f(x)=\frac{1}{2}{x}^{T}Ax?x^{T}b+c$$

可以發現：

$$?f=\mathrm{grad}f=Ax?b$$

由$A$的正定性可以保證$f(x)$的駐點一定是極小值點。而$Ax?b=0$得到的就是$f(x)$的駐點，即：

$$f(x^{?})=\min f(x)??f(x^{?})=A{x}^{?}?b=0$$

把解線性方程組的問題，轉化為求函數$f(x)$的極小值點。

二、最速下降法

怎么找到這個極小值點?

已知一個多元函數沿其負梯度方向函數值下降得最快。

一種較為形象的解釋：

想象自己在半山腰上，要到山腳處：

• 首先要找好下降方向：負梯度方向
• 之后沿著選定方向直走
• 走到不能再下降為止（也就是選定方向的最低點），停下來，再找新的下降方向
• 重復上面的過程，便能到達山腳

翻譯成數學語言

• 給定任意初值$x_0$，計算殘量$r_0=b?A{x}_0$。

• 選擇$P=r_0$為前進方向，計算：

  $$\alpha =\frac{\left(r_0,r_0\right)}{\left(A{r}_0,r_0\right)},x_1=x_0+\alpha r_0$$

• 重復上面的過程。

算法如下：

三、北太天元 or matlab實現

最速下降法解線性方程組

function [x,k,r] = Gradient_Descent(A,b,x0,e_tol,N)% 最速下降法 解線性方程組% Input: A, b(列向量), x0(初始值),e_tol: error tolerant, N: 限制迭代次數小于 N 次             % Output: x , k(迭代次數), r%   Version:            1.0%   last modified:      2024/05/19n = length(b); k = 0; R = zeros(1,N); % 記錄殘差r = b - A*x0;x = zeros(n,N); % 記錄每次迭代結果x(:,1) = x0;norm_r = norm(r,2); R(1) = norm_r;while norm_r > e_tol && k < Nalpha = r'*r/(r'*A*r);  % 計算步長x(:,k+2) = x(:,k+1) + alpha * r;r = b - A * x(:,k+2); % 殘量norm_r = norm(r,2);R(k+2)=norm_r;k = k+1;endx = x(:,1:k+1); % 返回每次的迭代結果r = R(1:k+1);   % 返回每次的殘差if k>Nfprintf('迭代超出最大迭代次數');elsefprintf('迭代次數=%i\n',k);end
end

四、數值算例

下面例子中統一 $ N=100,e_tol=10^{-8},x_0 = 0 $

例1

$$Ax=b$$
$$A=\left[\begin{array}{cccc}4& 1& 0& 0\\ 1& 4& 1& 0\\ 0& 1& 4& 1\\ 0& 0& 1& 4\end{array}\right]b=\left[\begin{array}{c}6\\ 25\\ ?11\\ 15\end{array}\right]$$

用最速下降法求$x$ ;

實現

clc;clear all,format long;
N = 100; e_tol = 1e-8;
A = [4, 1, 0, 0; 1, 4, 1, 0;  0, 1, 4, 1; 0, 0, 1, 4];
b = [6; 25; -11; 15];
x0 = [0; 0; 0; 0];
[x11, k1, r11] = Gradient_Descent(A, b, x0, e_tol, N);
x_exact = gsem_column(A, b);% 作圖查看誤差變化
n = length(b);
k1 = k1 + 1;% 數值解
figure(1);
plot(1:k1, x11(1,:), '-*r', 1:k1, x11(2,:), '-og', 1:k1, x11(3,:), '-+b', 1:k1, x11(4,:), '-dk');
legend('x_1', 'x_2', 'x_3', 'x_4');
title('每個數值解的變化');% 殘差變化
figure(2);
plot(1:k1, r11, '-*r');
legend('殘差');
title('殘差變化');

運行后得到

通過這個例子可以初步看到方法是可行的.

例2

下面這個例子我將形象展示最速下降法的實現特點

A = [3 1; 1 5];
b = [-1;1];
c = 0;

對應函數：$$f(x,y)=\frac{1}{2}\left(3x^2+2?1?xy+5y^2\right)?(?x+y)+0$$

三維表示一下

clc;clear all;format long;
A = [3  1; 1 5]; 
b = [-1;1];
c = 0; 
N = 100; e_tol = 1e-8; x0 =zeros(length(b),1);%x0 =[-0.1;-0.1]x = linspace(-1,1,100); 
y = linspace(-1,1,100);
% 網格化、方便作圖
[x, y] = meshgrid(x,y);
% 定義函數 f(x) = 0.5 * x' * A * x - x'*b + c
% 為了作圖方便,如下定義
f=@(x,y) 0.5 * (A(1,1) * x.^2 + 2 * A(1,2) * x .* y + A(2,2) * y.^2) - (b(1) * x + b(2) * y) + c;
z = f(x,y);mesh(x,y,z)
[x11,k1,r11] = Gradient_Descent(A,b,x0,e_tol,N);figure(1)
mesh(x,y,z)
hold on
% 繪制最速下降法的每次迭代點
%scatter3(x11(1, :), x11(2, :), f(x11(1,:),x11(2,:)),'r','filled');
plot3(x11(1, :), x11(2, :), f(x11(1,:),x11(2,:)),'r-o');xlabel('x');
ylabel('y');
zlabel('f(x, y)');
title('函數的三維表示');
hold off;

運行后得到

繪制等高線圖

figure(2)
hold on 
contour(x,y,z,200)
plot(x11(1, :), x11(2, :), 'r-o');
title('最速下降法特點');
colorbar;

運行后得到

為了展示更清晰,將 $ x_0 $設為 [-0.2;0] ,可以得到這樣的圖像

由圖形可以看出,最速下降法是如何下降的.

從某一點,選定最快的下降方向,下降到不能再下降為止,再重新找新的最快的下降方向.就這樣依次進行下去.

由此可以看出最速下降法的優點是容易理解和實現較為簡單.

當然也可以看出它還存在很大的改進空間,在每一次選方向時,明明有著更快更好的方向(三角形任意的第三邊都更快).

---

以上圖形均在北太天元軟件中繪制，matlab同樣可以正常運行。