MATLAB Function Reference

bicg

雙共軛梯度（Gradients）法

Syntax

x = bicg(A,b)
bicg(A,b,tol)
bicg(A,b,tol,maxit)
bicg(A,b,tol,maxit,M)
bicg(A,b,tol,maxit,M1,M2)
bicg(A,b,tol,maxit,M1,M2,x0)
bicg(afun,b,tol,maxit,mfun1,mfun2,x0,p1,p2,...)
[x,flag] = bicg(A,b,...)
[x,flag,relres] = bicg(A,b,...)
[x,flag,relres,iter] = bicg(A,b,...)
[x,flag,relres,iter,resvec] = bicg(A,b,...)

Description

x = bicg(A,b) 解線性方程式 A*x = b 的 x。A 為 n*n 的係數矩陣而行向量 b 的長度也必須為 n。A 可為函數 afun 且 afun(x) 回傳 A*x，afun(x,'transp') 回傳 A'*x。

若 bicg 收歛，相關訊息將顯示出來。若 bicg 經過最大次數的重覆或任何理由被停止以致於沒有收歛，將會列印出一警告訊息顯示相關剩餘 norm(b-A*x)/norm(b)，而該作法下的重覆數字也將停止。

bicg(A,b,tol 定義方法的誤差。若 tol 為 []，則 bicg 使用預設值 1e-6。

bicg(A,b,tol,maxit) 定義最大重覆次數。若 maxit 為 []，則 bicg 使用預設值 min(n,20)。

bicg(A,b,tol,maxit,M) 及 bicg(A,b,tol,maxit,M1,M2) 使用先決條件 M 或 M = M1*M2 並有效地解方程式 inv(M)*A*x = inv(M)*b 的 x。若 M 為 []，則 bicg 不使用任何先決控制條件。M 可為函式 mfun 其可回傳 M\x 且 mfun(x,'transp') 回傳 M'\x。

bicg(A,b,tol,maxit,M1,M2,x0) 定義初始的推測。若 x0 為 []，則 bicg 使用預設值。即全為零的向量。

bicg(afun,b,tol,maxit,m1fun,m2fun,x0,p1,p2,...) 傳送參數 p1,p2,... 至函式 afun(x,p1,p2,...) 及 afun(x,p1,p2,...,'transp')，亦與先決函式 m1fun 及 m2fun 相似。

[x,flag] = bicg(A,b,...) 亦回傳收歛旗標。

旗標	收歛性
0	bicg 在 maxit 次內收歛至期望的誤差 tol。
1	bicg 重覆 maxit 次但並沒有收歛。
2	先決條件 M 不夠完善。
3	bicg 沉滯（stagnated）。（兩次連續的重覆視為相同）。
4	經 bicg 計算的數值變得太小或太大以致於沒辦法繼續做計算。

當 flag 不為 0 時，解 x 回傳經由所有重覆次數的最小基準餘數（minimal norm residual）。若有定義 flag 的輸出，則不會顯示任何訊息。

[x,flag,relres] = bicg(A,b,...) 亦回傳相關剩餘 norm(b-A*x)/norm(b)。若 flag 為 0, relres <= tol。

[x,flag,relres,iter] = bicg(A,b,...) 亦回傳 x 計算時的重覆次數，即 0  iter maxit。

[x,flag,relres,iter,resvec] = bicg(A,b,...) 亦回傳在每一次重覆的剩餘基準估計向量，包括 norm(b-A*x0)。

Examples

範例 1

n = 100; 
on = ones(n,1); 
A = spdiags([-2*on 4*on -on],-1:1,n,n);
b = sum(A,2); 
tol = 1e-8; 
maxit = 15;
M1 = spdiags([on/(-2) on],-1:0,n,n); 
M2 = spdiags([4*on -on],0:1,n,n);

x = bicg(A,b,tol,maxit,M1,M2,[]);
bicg converged at iteration 9 to a solution with relative 
residual 5.3e-009

另外選擇使用矩陣向量乘法的函式

function y = afun(x,n,transp_flag)
if (nargin > 2) & strcmp(transp_flag,'transp')
    y = 4 * x;
    y(1:n-1) = y(1:n-1) - 2 * x(2:n);
    y(2:n) = y(2:n) - x(1:n-1);
else
    y = 4 * x;
    y(2:n) = y(2:n) - 2 * x(1:n-1);
    y(1:n-1) = y(1:n-1) - x(2:n);
end

如同 bicg 的輸入。

   x1 = bicg(@afun,b,tol,maxit,M1,M2,[],n);

範例 2 以 A = west0479 為開始。

load west0479;
A = west0479;
b = sum(A,2);

若 A 不是很大，則可精確地使用反斜線來解 A*x = b。

x = A \ b;
norm(b-A*x) / norm(b)

ans =
    1.2454e-017

現在試著用 bicg 來解 A*x = b。

[x,flag,relres,iter,resvec] = bicg(A,b)

flag =
          1
relres =
          1
iter =
          0

flag 的值表示 bicg 經預設的 20 次重覆後沒有收歛。iter 的值顯示初始化全零的推測勝於隨後的重覆。relres 的值支援：relres = norm(b-A*x)/norm(b) = norm(b)/norm(b) = 1。

semilogy(0:20,resvec/norm(b),'-o')
xlabel('iteration number')
ylabel('relative residual')

以下試著用先決條件的微誤差為 1e-5 的不完全 LU 分解。

[L1,U1] = luinc(A,1e-5);
Warning: Incomplete upper triangular factor has 1 zero diagonal.
         It cannot be used as a preconditioner for an iterative 
         method.

nnz(A) 
ans =
         1887
nnz(L1) 
ans =
         5562
nnz(U1) 
ans =
         4320

在上三角 U1 主對角線上的零表示 U1 為單一的。

[x,flag,relres,iter,resvec] = bicg(A,b,1e-6,20,L1,U1)

flag =
         2
relres =
         1
iter =
         0
resvec =
         7.0557e+005

此方法在第一次重覆去使用反斜線解方程式時會失敗。若無其他重覆會產生，bicg 將回傳初始估計值。

再試著用輕微的稀疏先決條件。

[L2,U2] = luinc(A,1e-6)

nnz(L2) 
ans =
         6231
nnz(U2) 
ans =
         4559

此時 U2 不為單一且可能是個合適的先決條件。

[x,flag,relres,iter,resvec] = bicg(A,b,1e-15,10,L2,U2)

flag =
         0
relres =
         2.0248e-16
iter =
         8

bicg 在第八次重覆後會收歛至期望的誤差值。漸減微誤差的值將增加不完全因子的填滿，但也會增加原先矩陣近似值的精確度。因此，先決系統變得更近似於 inv(U)*inv(L)*L*U*x = inv(U)*inv(L)*b，其中 L 和 U 為 LU 的因子。

下圖顯示 bicg 使用來種不同的不完全 LU 因子當做先決條件的過程。圖中的每一行用先決條件的微誤差來標示。

下圖繪製出不完全 LU 因子相對於時間的微誤差，重覆一次先決條件也將被計算一次，而他們的總和，即總時間對將複用來解該問題。針對產生因子的時間並不會於填滿時非常快地增加，反而會慢下來。west0479 為一相當小的矩陣，僅 139*139，而先決條件 bicg 仍然花比反斜線還多的時間。

See Also

bicgstab, cgs, gmres, lsqr, luinc, minres, pcg, qmr, symmlq

@ (function handle), \ (backslash)

References

[1]  Barrett, R., M. Berry, T. F. Chan, et al., Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, SIAM, Philadelphia, 1994.

beta, betainc, betaln

bicgstab